From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/float/floorf.c | 526 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 526 insertions(+) create mode 100644 libm/float/floorf.c (limited to 'libm/float/floorf.c') diff --git a/libm/float/floorf.c b/libm/float/floorf.c new file mode 100644 index 000000000..7a2f3530d --- /dev/null +++ b/libm/float/floorf.c @@ -0,0 +1,526 @@ +/* ceilf() + * floorf() + * frexpf() + * ldexpf() + * signbitf() + * isnanf() + * isfinitef() + * + * Single precision floating point numeric utilities + * + * + * + * SYNOPSIS: + * + * float x, y; + * float ceilf(), floorf(), frexpf(), ldexpf(); + * int signbit(), isnan(), isfinite(); + * int expnt, n; + * + * y = floorf(x); + * y = ceilf(x); + * y = frexpf( x, &expnt ); + * y = ldexpf( x, n ); + * n = signbit(x); + * n = isnan(x); + * n = isfinite(x); + * + * + * + * DESCRIPTION: + * + * All four routines return a single precision floating point + * result. + * + * sfloor() returns the largest integer less than or equal to x. + * It truncates toward minus infinity. + * + * sceil() returns the smallest integer greater than or equal + * to x. It truncates toward plus infinity. + * + * sfrexp() 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. + * + * ldexpf() 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.1: December, 1988 +Copyright 1984, 1987, 1988 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + + +#include +#ifdef DEC +#undef DENORMAL +#define DENORMAL 0 +#endif + +#ifdef UNK +#undef UNK +#if BIGENDIAN +#define MIEEE 1 +#else +#define IBMPC 1 +#endif +/* +char *unkmsg = "ceil(), floor(), frexp(), ldexp() must be rewritten!\n"; +*/ +#endif + +#define EXPMSK 0x807f +#define MEXP 255 +#define NBITS 24 + + +extern float MAXNUMF; /* (2^24 - 1) * 2^103 */ +#ifdef ANSIC +float floorf(float); +#else +float floorf(); +#endif + +float ceilf( float x ) +{ +float y; + +#ifdef UNK +printf( "%s\n", unkmsg ); +return(0.0); +#endif + +y = floorf( (float )x ); +if( y < x ) + y += 1.0; +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, +}; + + + +float floorf( float x ) +{ +unsigned short *p; +union + { + float y; + unsigned short i[2]; + } u; +int e; + +#ifdef UNK +printf( "%s\n", unkmsg ); +return(0.0); +#endif + +u.y = x; +/* find the exponent (power of 2) */ +#ifdef DEC +p = &u.i[0]; +e = (( *p >> 7) & 0377) - 0201; +p += 3; +#endif + +#ifdef IBMPC +p = &u.i[1]; +e = (( *p >> 7) & 0xff) - 0x7f; +p -= 1; +#endif + +#ifdef MIEEE +p = &u.i[0]; +e = (( *p >> 7) & 0xff) - 0x7f; +p += 1; +#endif + +if( e < 0 ) + { + if( u.y < 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); +} + + + +float frexpf( float x, int *pw2 ) +{ +union + { + float y; + unsigned short i[2]; + } u; +int i, k; +short *q; + +u.y = x; + +#ifdef UNK +printf( "%s\n", unkmsg ); +return(0.0); +#endif + +#ifdef IBMPC +q = &u.i[1]; +#endif + +#ifdef DEC +q = &u.i[0]; +#endif + +#ifdef MIEEE +q = &u.i[0]; +#endif + +/* find the exponent (power of 2) */ + +i = ( *q >> 7) & 0xff; +if( i == 0 ) + { + if( u.y == 0.0 ) + { + *pw2 = 0; + return(0.0); + } +/* Number is denormal or zero */ +#if DENORMAL +/* Handle denormal number. */ + do + { + u.y *= 2.0; + i -= 1; + k = ( *q >> 7) & 0xff; + } + while( k == 0 ); + i = i + k; +#else + *pw2 = 0; + return( 0.0 ); +#endif /* DENORMAL */ + } +i -= 0x7e; +*pw2 = i; +*q &= 0x807f; /* strip all exponent bits */ +*q |= 0x3f00; /* mantissa between 0.5 and 1 */ +return( u.y ); +} + + + + + +float ldexpf( float x, int pw2 ) +{ +union + { + float y; + unsigned short i[2]; + } u; +short *q; +int e; + +#ifdef UNK +printf( "%s\n", unkmsg ); +return(0.0); +#endif + +u.y = x; +#ifdef DEC +q = &u.i[0]; +#endif + +#ifdef IBMPC +q = &u.i[1]; +#endif +#ifdef MIEEE +q = &u.i[0]; +#endif +while( (e = ( *q >> 7) & 0xff) == 0 ) + { + if( u.y == (float )0.0 ) + { + return( 0.0 ); + } +/* Input is denormal. */ + if( pw2 > 0 ) + { + u.y *= 2.0; + pw2 -= 1; + } + if( pw2 < 0 ) + { + if( pw2 < -24 ) + return( 0.0 ); + u.y *= 0.5; + pw2 += 1; + } + if( pw2 == 0 ) + return(u.y); + } + +e += pw2; + +/* Handle overflow */ +if( e > MEXP ) + { + return( MAXNUMF ); + } + +*q &= 0x807f; + +/* Handle denormalized results */ +if( e < 1 ) + { +#if DENORMAL + if( e < -24 ) + return( 0.0 ); + *q |= 0x80; /* Set LSB of exponent. */ + /* 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 *= ldexpf(1.0f, e - 1); + return(u.y); +#else + return( 0.0 ); +#endif + } +*q |= (e & 0xff) << 7; +return(u.y); +} + + +/* Return 1 if the sign bit of x is 1, else 0. */ + +int signbitf(x) +float x; +{ +union + { + float f; + short s[4]; + int i; + } u; + +u.f = x; + +if( sizeof(int) == 4 ) + { +#ifdef IBMPC + return( u.i < 0 ); +#endif +#ifdef DEC + return( u.s[1] < 0 ); +#endif +#ifdef MIEEE + return( u.i < 0 ); +#endif + } +else + { +#ifdef IBMPC + return( u.s[1] < 0 ); +#endif +#ifdef DEC + return( u.s[1] < 0 ); +#endif +#ifdef MIEEE + return( u.s[0] < 0 ); +#endif + } +} + + +/* Return 1 if x is a number that is Not a Number, else return 0. */ + +int isnanf(x) +float x; +{ +#ifdef NANS +union + { + float f; + unsigned short s[2]; + unsigned int i; + } u; + +u.f = x; + +if( sizeof(int) == 4 ) + { +#ifdef IBMPC + if( ((u.i & 0x7f800000) == 0x7f800000) + && ((u.i & 0x007fffff) != 0) ) + return 1; +#endif +#ifdef DEC + if( (u.s[1] & 0x7f80) == 0) + { + if( (u.s[1] | u.s[0]) != 0 ) + return(1); + } +#endif +#ifdef MIEEE + if( ((u.i & 0x7f800000) == 0x7f800000) + && ((u.i & 0x007fffff) != 0) ) + return 1; +#endif + return(0); + } +else + { /* size int not 4 */ +#ifdef IBMPC + if( (u.s[1] & 0x7f80) == 0x7f80) + { + if( ((u.s[1] & 0x007f) | u.s[0]) != 0 ) + return(1); + } +#endif +#ifdef DEC + if( (u.s[1] & 0x7f80) == 0) + { + if( (u.s[1] | u.s[0]) != 0 ) + return(1); + } +#endif +#ifdef MIEEE + if( (u.s[0] & 0x7f80) == 0x7f80) + { + if( ((u.s[0] & 0x000f) | u.s[1]) != 0 ) + return(1); + } +#endif + return(0); + } /* size int not 4 */ + +#else +/* No NANS. */ +return(0); +#endif +} + + +/* Return 1 if x is not infinite and is not a NaN. */ + +int isfinitef(x) +float x; +{ +#ifdef INFINITIES +union + { + float f; + unsigned short s[2]; + unsigned int i; + } u; + +u.f = x; + +if( sizeof(int) == 4 ) + { +#ifdef IBMPC + if( (u.i & 0x7f800000) != 0x7f800000) + return 1; +#endif +#ifdef DEC + if( (u.s[1] & 0x7f80) == 0) + { + if( (u.s[1] | u.s[0]) != 0 ) + return(1); + } +#endif +#ifdef MIEEE + if( (u.i & 0x7f800000) != 0x7f800000) + return 1; +#endif + return(0); + } +else + { +#ifdef IBMPC + if( (u.s[1] & 0x7f80) != 0x7f80) + return 1; +#endif +#ifdef DEC + if( (u.s[1] & 0x7f80) == 0) + { + if( (u.s[1] | u.s[0]) != 0 ) + return(1); + } +#endif +#ifdef MIEEE + if( (u.s[0] & 0x7f80) != 0x7f80) + return 1; +#endif + return(0); + } +#else +/* No INFINITY. */ +return(1); +#endif +} -- cgit v1.2.3