Totally rework the math library, this time based on the MacOs X

math library (which is itself based on the math lib from FreeBSD).
 -Erik

1 files changed, 0 insertions, 432 deletions

diff --git a/libm/ldouble/floorl.c b/libm/ldouble/floorl.c
deleted file mode 100644
index 1abdfb2cd..000000000
--- a/libm/ldouble/floorl.c
+++ /dev/null
@@ -1,432 +0,0 @@
-/*							ceill()
- *							floorl()
- *							frexpl()
- *							ldexpl()
- *							fabsl()
- *							signbitl()
- *							isnanl()
- *							isfinitel()
- *
- *	Floating point numeric utilities
- *
- *
- *
- * SYNOPSIS:
- *
- * long double ceill(), floorl(), frexpl(), ldexpl(), fabsl();
- * int signbitl(), isnanl(), isfinitel();
- * long double x, y;
- * int expnt, n;
- *
- * y = floorl(x);
- * y = ceill(x);
- * y = frexpl( x, &expnt );
- * y = ldexpl( x, n );
- * y = fabsl( x );
- * n = signbitl(x);
- * n = isnanl(x);
- * n = isfinitel(x);
- *
- *
- *
- * DESCRIPTION:
- *
- * The following routines return a long double precision floating point
- * result:
- *
- * floorl() returns the largest integer less than or equal to x.
- * It truncates toward minus infinity.
- *
- * ceill() returns the smallest integer greater than or equal
- * to x.  It truncates toward plus infinity.
- *
- * frexpl() extracts the exponent from x.  It returns an integer
- * power of two to expnt and the significand between 0.5 and 1
- * to y.  Thus  x = y * 2**expn.
- *
- * ldexpl() multiplies x by 2**n.
- *
- * fabsl() returns the absolute value of its argument.
- *
- * These functions are part of the standard C run time library
- * for some but not all C compilers.  The ones supplied are
- * written in C for IEEE arithmetic.  They should
- * be used only if your compiler library does not already have
- * them.
- *
- * The IEEE versions assume that denormal numbers are implemented
- * in the arithmetic.  Some modifications will be required if
- * the arithmetic has abrupt rather than gradual underflow.
- */
-
-
-/*
-Cephes Math Library Release 2.7:  May, 1998
-Copyright 1984, 1987, 1988, 1992, 1998 by Stephen L. Moshier
-*/
-
-
-#include <math.h>
-
-/* This is defined in mconf.h. */
-/* #define DENORMAL 1 */
-
-#ifdef UNK
-/* Change UNK into something else.  */
-#undef UNK
-#if BIGENDIAN
-#define MIEEE 1
-#else
-#define IBMPC 1
-#endif
-#endif
-
-#ifdef IBMPC
-#define EXPMSK 0x800f
-#define MEXP 0x7ff
-#define NBITS 64
-#endif
-
-#ifdef MIEEE
-#define EXPMSK 0x800f
-#define MEXP 0x7ff
-#define NBITS 64
-#endif
-
-extern double MAXNUML;
-
-#ifdef ANSIPROT
-extern long double fabsl ( long double );
-extern long double floorl ( long double );
-extern int isnanl ( long double );
-#else
-long double fabsl(), floorl();
-int isnanl();
-#endif
-#ifdef INFINITIES
-extern long double INFINITYL;
-#endif
-#ifdef NANS
-extern long double NANL;
-#endif
-
-long double fabsl(x)
-long double x;
-{
-union
-  {
-    long double d;
-    short i[6];
-  } u;
-
-u.d = x;
-#ifdef IBMPC
-    u.i[4] &= 0x7fff;
-#endif
-#ifdef MIEEE
-    u.i[0] &= 0x7fff;
-#endif
-return( u.d );
-}
-
-
-
-long double ceill(x)
-long double x;
-{
-long double y;
-
-#ifdef UNK
-mtherr( "ceill", DOMAIN );
-return(0.0L);
-#endif
-#ifdef INFINITIES
-if(x == -INFINITYL)
-	return(x);
-#endif
-#ifdef MINUSZERO
-if(x == 0.0L)
-	return(x);
-#endif
-y = floorl(x);
-if( y < x )
-	y += 1.0L;
-return(y);
-}
-
-
-
-
-/* Bit clearing masks: */
-
-static unsigned short bmask[] = {
-0xffff,
-0xfffe,
-0xfffc,
-0xfff8,
-0xfff0,
-0xffe0,
-0xffc0,
-0xff80,
-0xff00,
-0xfe00,
-0xfc00,
-0xf800,
-0xf000,
-0xe000,
-0xc000,
-0x8000,
-0x0000,
-};
-
-
-
-
-long double floorl(x)
-long double x;
-{
-unsigned short *p;
-union
-  {
-    long double y;
-    unsigned short sh[6];
-  } u;
-int e;
-
-#ifdef UNK
-mtherr( "floor", DOMAIN );
-return(0.0L);
-#endif
-#ifdef INFINITIES
-if( x == INFINITYL )
-	return(x);
-#endif
-#ifdef MINUSZERO
-if(x == 0.0L)
-	return(x);
-#endif
-u.y = x;
-/* find the exponent (power of 2) */
-#ifdef IBMPC
-p = (unsigned short *)&u.sh[4];
-e = (*p & 0x7fff) - 0x3fff;
-p -= 4;
-#endif
-
-#ifdef MIEEE
-p = (unsigned short *)&u.sh[0];
-e = (*p & 0x7fff) - 0x3fff;
-p += 5;
-#endif
-
-if( e < 0 )
-	{
-	if( u.y < 0.0L )
-		return( -1.0L );
-	else
-		return( 0.0L );
-	}
-
-e = (NBITS -1) - e;
-/* clean out 16 bits at a time */
-while( e >= 16 )
-	{
-#ifdef IBMPC
-	*p++ = 0;
-#endif
-
-#ifdef MIEEE
-	*p-- = 0;
-#endif
-	e -= 16;
-	}
-
-/* clear the remaining bits */
-if( e > 0 )
-	*p &= bmask[e];
-
-if( (x < 0) && (u.y != x) )
-	u.y -= 1.0L;
-
-return(u.y);
-}
-
-
-
-long double frexpl( x, pw2 )
-long double x;
-int *pw2;
-{
-union
-  {
-    long double y;
-    unsigned short sh[6];
-  } u;
-int i, k;
-short *q;
-
-u.y = x;
-
-#ifdef NANS
-if(isnanl(x))
-	{
-	*pw2 = 0;
-	return(x);
-	}
-#endif
-#ifdef INFINITIES
-if(x == -INFINITYL)
-	{
-	*pw2 = 0;
-	return(x);
-	}
-#endif
-#ifdef MINUSZERO
-if(x == 0.0L)
-	{
-	*pw2 = 0;
-	return(x);
-	}
-#endif
-
-#ifdef UNK
-mtherr( "frexpl", DOMAIN );
-return(0.0L);
-#endif
-
-/* find the exponent (power of 2) */
-#ifdef IBMPC
-q = (short *)&u.sh[4];
-i  = *q & 0x7fff;
-#endif
-
-#ifdef MIEEE
-q = (short *)&u.sh[0];
-i  = *q & 0x7fff;
-#endif
-
-if( i == 0 )
-	{
-	if( u.y == 0.0L )
-		{
-		*pw2 = 0;
-		return(0.0L);
-		}
-/* Number is denormal or zero */
-#ifdef DENORMAL
-/* Handle denormal number. */
-do
-	{
-	u.y *= 2.0L;
-	i -= 1;
-	k  = *q & 0x7fff;
-	}
-while( (k == 0) && (i > -66) );
-i = i + k;
-#else
-	*pw2 = 0;
-	return(0.0L);
-#endif /* DENORMAL */
-	}
-
-*pw2 = i - 0x3ffe;
-/* *q = 0x3ffe; */
-/* Preserve sign of argument.  */
-*q &= 0x8000;
-*q |= 0x3ffe;
-return( u.y );
-}
-
-
-
-
-
-
-long double ldexpl( x, pw2 )
-long double x;
-int pw2;
-{
-union
-  {
-    long double y;
-    unsigned short sh[6];
-  } u;
-unsigned short *q;
-long e;
-
-#ifdef UNK
-mtherr( "ldexp", DOMAIN );
-return(0.0L);
-#endif
-
-u.y = x;
-#ifdef IBMPC
-q = (unsigned short *)&u.sh[4];
-#endif
-#ifdef MIEEE
-q = (unsigned short *)&u.sh[0];
-#endif
-while( (e = (*q & 0x7fffL)) == 0 )
-	{
-#ifdef DENORMAL
-	if( u.y == 0.0L )
-		{
-		return( 0.0L );
-		}
-/* Input is denormal. */
-	if( pw2 > 0 )
-		{
-		u.y *= 2.0L;
-		pw2 -= 1;
-		}
-	if( pw2 < 0 )
-		{
-		if( pw2 < -64 )
-			return(0.0L);
-		u.y *= 0.5L;
-		pw2 += 1;
-		}
-	if( pw2 == 0 )
-		return(u.y);
-#else
-	return( 0.0L );
-#endif
-	}
-
-e = e + pw2;
-
-/* Handle overflow */
-if( e > 0x7fffL )
-	{
-	return( MAXNUML );
-	}
-*q &= 0x8000;
-/* Handle denormalized results */
-if( e < 1 )
-	{
-#ifdef DENORMAL
-	if( e < -64 )
-		return(0.0L);
-
-#ifdef IBMPC
-	*(q-1) |= 0x8000;
-#endif
-#ifdef MIEEE
-	*(q+2) |= 0x8000;
-#endif
-
-	while( e < 1 )
-		{
-		u.y *= 0.5L;
-		e += 1;
-		}
-	e = 0;
-#else
-	return(0.0L);
-#endif
-	}
-
-*q |= (unsigned short) e & 0x7fff;
-return(u.y);
-}
-