diff options
Diffstat (limited to 'libm/double/floor.c')
-rw-r--r-- | libm/double/floor.c | 531 |
1 files changed, 0 insertions, 531 deletions
diff --git a/libm/double/floor.c b/libm/double/floor.c deleted file mode 100644 index affc7753e..000000000 --- a/libm/double/floor.c +++ /dev/null @@ -1,531 +0,0 @@ -/* ceil() - * floor() - * frexp() - * ldexp() - * signbit() - * isnan() - * isfinite() - * - * Floating point numeric utilities - * - * - * - * SYNOPSIS: - * - * double ceil(), floor(), frexp(), ldexp(); - * int signbit(), isnan(), isfinite(); - * double x, y; - * int expnt, n; - * - * y = floor(x); - * y = ceil(x); - * y = frexp( x, &expnt ); - * y = ldexp( x, n ); - * n = signbit(x); - * n = isnan(x); - * n = isfinite(x); - * - * - * - * DESCRIPTION: - * - * All four routines return a double precision floating point - * result. - * - * floor() returns the largest integer less than or equal to x. - * It truncates toward minus infinity. - * - * ceil() returns the smallest integer greater than or equal - * to x. It truncates toward plus infinity. - * - * frexp() 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. - * - * ldexp() multiplies x by 2**n. - * - * signbit(x) returns 1 if the sign bit of x is 1, else 0. - * - * These functions are part of the standard C run time library - * for many but not all C compilers. The ones supplied are - * written in C for either DEC or 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.8: June, 2000 -Copyright 1984, 1995, 2000 by Stephen L. Moshier -*/ - - -#include <math.h> - -#ifdef UNK -/* ceil(), floor(), frexp(), ldexp() may need to be rewritten. */ -#undef UNK -#if BIGENDIAN -#define MIEEE 1 -#else -#define IBMPC 1 -#endif -#endif - -#ifdef DEC -#define EXPMSK 0x807f -#define MEXP 255 -#define NBITS 56 -#endif - -#ifdef IBMPC -#define EXPMSK 0x800f -#define MEXP 0x7ff -#define NBITS 53 -#endif - -#ifdef MIEEE -#define EXPMSK 0x800f -#define MEXP 0x7ff -#define NBITS 53 -#endif - -extern double MAXNUM, NEGZERO; -#ifdef ANSIPROT -double floor ( double ); -int isnan ( double ); -int isfinite ( double ); -double ldexp ( double, int ); -#else -double floor(); -int isnan(), isfinite(); -double ldexp(); -#endif - -double ceil(x) -double x; -{ -double y; - -#ifdef UNK -mtherr( "ceil", DOMAIN ); -return(0.0); -#endif -#ifdef NANS -if( isnan(x) ) - return( x ); -#endif -#ifdef INFINITIES -if(!isfinite(x)) - return(x); -#endif - -y = floor(x); -if( y < x ) - y += 1.0; -#ifdef MINUSZERO -if( y == 0.0 && x < 0.0 ) - return( NEGZERO ); -#endif -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, -}; - - - - - -double floor(x) -double x; -{ -union - { - double y; - unsigned short sh[4]; - } u; -unsigned short *p; -int e; - -#ifdef UNK -mtherr( "floor", DOMAIN ); -return(0.0); -#endif -#ifdef NANS -if( isnan(x) ) - return( x ); -#endif -#ifdef INFINITIES -if(!isfinite(x)) - return(x); -#endif -#ifdef MINUSZERO -if(x == 0.0L) - return(x); -#endif -u.y = x; -/* find the exponent (power of 2) */ -#ifdef DEC -p = (unsigned short *)&u.sh[0]; -e = (( *p >> 7) & 0377) - 0201; -p += 3; -#endif - -#ifdef IBMPC -p = (unsigned short *)&u.sh[3]; -e = (( *p >> 4) & 0x7ff) - 0x3ff; -p -= 3; -#endif - -#ifdef MIEEE -p = (unsigned short *)&u.sh[0]; -e = (( *p >> 4) & 0x7ff) - 0x3ff; -p += 3; -#endif - -if( e < 0 ) - { - if( u.y < 0.0 ) - return( -1.0 ); - else - return( 0.0 ); - } - -e = (NBITS -1) - e; -/* clean out 16 bits at a time */ -while( e >= 16 ) - { -#ifdef IBMPC - *p++ = 0; -#endif - -#ifdef DEC - *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.0; - -return(u.y); -} - - - - -double frexp( x, pw2 ) -double x; -int *pw2; -{ -union - { - double y; - unsigned short sh[4]; - } u; -int i; -#ifdef DENORMAL -int k; -#endif -short *q; - -u.y = x; - -#ifdef UNK -mtherr( "frexp", DOMAIN ); -return(0.0); -#endif - -#ifdef IBMPC -q = (short *)&u.sh[3]; -#endif - -#ifdef DEC -q = (short *)&u.sh[0]; -#endif - -#ifdef MIEEE -q = (short *)&u.sh[0]; -#endif - -/* find the exponent (power of 2) */ -#ifdef DEC -i = ( *q >> 7) & 0377; -if( i == 0 ) - { - *pw2 = 0; - return(0.0); - } -i -= 0200; -*pw2 = i; -*q &= 0x807f; /* strip all exponent bits */ -*q |= 040000; /* mantissa between 0.5 and 1 */ -return(u.y); -#endif - -#ifdef IBMPC -i = ( *q >> 4) & 0x7ff; -if( i != 0 ) - goto ieeedon; -#endif - -#ifdef MIEEE -i = *q >> 4; -i &= 0x7ff; -if( i != 0 ) - goto ieeedon; -#ifdef DENORMAL - -#else -*pw2 = 0; -return(0.0); -#endif - -#endif - - -#ifndef DEC -/* Number is denormal or zero */ -#ifdef DENORMAL -if( u.y == 0.0 ) - { - *pw2 = 0; - return( 0.0 ); - } - - -/* Handle denormal number. */ -do - { - u.y *= 2.0; - i -= 1; - k = ( *q >> 4) & 0x7ff; - } -while( k == 0 ); -i = i + k; -#endif /* DENORMAL */ - -ieeedon: - -i -= 0x3fe; -*pw2 = i; -*q &= 0x800f; -*q |= 0x3fe0; -return( u.y ); -#endif -} - - - - - - - -double ldexp( x, pw2 ) -double x; -int pw2; -{ -union - { - double y; - unsigned short sh[4]; - } u; -short *q; -int e; - -#ifdef UNK -mtherr( "ldexp", DOMAIN ); -return(0.0); -#endif - -u.y = x; -#ifdef DEC -q = (short *)&u.sh[0]; -e = ( *q >> 7) & 0377; -if( e == 0 ) - return(0.0); -#else - -#ifdef IBMPC -q = (short *)&u.sh[3]; -#endif -#ifdef MIEEE -q = (short *)&u.sh[0]; -#endif -while( (e = (*q & 0x7ff0) >> 4) == 0 ) - { - if( u.y == 0.0 ) - { - return( 0.0 ); - } -/* Input is denormal. */ - if( pw2 > 0 ) - { - u.y *= 2.0; - pw2 -= 1; - } - if( pw2 < 0 ) - { - if( pw2 < -53 ) - return(0.0); - u.y /= 2.0; - pw2 += 1; - } - if( pw2 == 0 ) - return(u.y); - } -#endif /* not DEC */ - -e += pw2; - -/* Handle overflow */ -#ifdef DEC -if( e > MEXP ) - return( MAXNUM ); -#else -if( e >= MEXP ) - return( 2.0*MAXNUM ); -#endif - -/* Handle denormalized results */ -if( e < 1 ) - { -#ifdef DENORMAL - if( e < -53 ) - return(0.0); - *q &= 0x800f; - *q |= 0x10; - /* For denormals, significant bits may be lost even - when dividing by 2. Construct 2^-(1-e) so the result - is obtained with only one multiplication. */ - u.y *= ldexp(1.0, e-1); - return(u.y); -#else - return(0.0); -#endif - } -else - { -#ifdef DEC - *q &= 0x807f; /* strip all exponent bits */ - *q |= (e & 0xff) << 7; -#else - *q &= 0x800f; - *q |= (e & 0x7ff) << 4; -#endif - return(u.y); - } -} - -/**********************************************************************/ -/* - * trunc is just a slightly modified version of floor above. - */ - -double trunc(double x) -{ - union { - double y; - unsigned short sh[4]; - } u; - unsigned short *p; - int e; - -#ifdef UNK - mtherr( "trunc", DOMAIN ); - return(0.0); -#endif -#ifdef NANS - if( isnan(x) ) - return( x ); -#endif -#ifdef INFINITIES - if(!isfinite(x)) - return(x); -#endif -#ifdef MINUSZERO - if(x == 0.0L) - return(x); -#endif - u.y = x; - /* find the exponent (power of 2) */ -#ifdef DEC - p = (unsigned short *)&u.sh[0]; - e = (( *p >> 7) & 0377) - 0201; - p += 3; -#endif - -#ifdef IBMPC - p = (unsigned short *)&u.sh[3]; - e = (( *p >> 4) & 0x7ff) - 0x3ff; - p -= 3; -#endif - -#ifdef MIEEE - p = (unsigned short *)&u.sh[0]; - e = (( *p >> 4) & 0x7ff) - 0x3ff; - p += 3; -#endif - - if( e < 0 ) - return( 0.0 ); - - e = (NBITS -1) - e; - /* clean out 16 bits at a time */ - while( e >= 16 ) - { -#ifdef IBMPC - *p++ = 0; -#endif - -#ifdef DEC - *p-- = 0; -#endif - -#ifdef MIEEE - *p-- = 0; -#endif - e -= 16; - } - - /* clear the remaining bits */ - if( e > 0 ) - *p &= bmask[e]; - - return(u.y); -} |