uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 432 insertions, 0 deletions

diff --git a/libm/ldouble/floorl.c b/libm/ldouble/floorl.c
new file mode 100644
index 000000000..1abdfb2cd
--- /dev/null
+++ b/libm/ldouble/floorl.c
@@ -0,0 +1,432 @@
+/*							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);
+}
+