summaryrefslogtreecommitdiff
path: root/libm/ldouble/floorl.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/ldouble/floorl.c')
-rw-r--r--libm/ldouble/floorl.c432
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);
+}
+