Index: libc/include/math.h =================================================================== RCS file: /cvs/src/src/newlib/libc/include/math.h,v retrieving revision 1.34 diff -p -u -r1.34 math.h --- libc/include/math.h 18 Mar 2009 18:11:18 -0000 1.34 +++ libc/include/math.h 20 Mar 2009 21:27:02 -0000 @@ -229,6 +229,7 @@ extern long int lrint _PARAMS((double)); extern _LONG_LONG_TYPE int llrint _PARAMS((double)); extern double round _PARAMS((double)); extern long int lround _PARAMS((double)); +extern long long int llround _PARAMS((double)); extern double trunc _PARAMS((double)); extern double remquo _PARAMS((double, double, int *)); extern double copysign _PARAMS((double, double)); @@ -297,6 +298,7 @@ extern long int lrintf _PARAMS((float)); extern _LONG_LONG_TYPE llrintf _PARAMS((float)); extern float roundf _PARAMS((float)); extern long int lroundf _PARAMS((float)); +extern long long int llroundf _PARAMS((float)); extern float truncf _PARAMS((float)); extern float remquof _PARAMS((float, float, int *)); extern float copysignf _PARAMS((float, float)); Index: libm/common/sf_lround.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/sf_lround.c,v retrieving revision 1.3 diff -p -u -r1.3 sf_lround.c --- libm/common/sf_lround.c 24 Aug 2007 20:49:59 -0000 1.3 +++ libm/common/sf_lround.c 20 Mar 2009 21:27:02 -0000 @@ -56,7 +56,7 @@ double x; #endif { - return (double) lroundf((float) x); + return lroundf((float) x); } #endif /* defined(_DOUBLE_IS_32BITS) */ Index: libm/common/sf_lrint.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/sf_lrint.c,v retrieving revision 1.3 diff -p -u -r1.3 sf_lrint.c --- libm/common/sf_lrint.c 8 Sep 2005 21:03:54 -0000 1.3 +++ libm/common/sf_lrint.c 20 Mar 2009 21:27:02 -0000 @@ -95,7 +95,7 @@ TWO23[2]={ double x; #endif { - return (double) lrintf((float) x); + return lrintf((float) x); } #endif /* defined(_DOUBLE_IS_32BITS) */ Index: libm/common/sf_fdim.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/sf_fdim.c,v retrieving revision 1.1 diff -p -u -r1.1 sf_fdim.c --- libm/common/sf_fdim.c 7 Jun 2002 21:59:56 -0000 1.1 +++ libm/common/sf_fdim.c 20 Mar 2009 21:27:02 -0000 @@ -15,8 +15,10 @@ #endif { int c = __fpclassifyf(x); - if (c == FP_NAN || c == FP_INFINITE) - return HUGE_VAL; + if (c == FP_NAN) return(x); + if (__fpclassifyf(y) == FP_NAN) return(y); + if (c == FP_INFINITE) + return HUGE_VALF; return x > y ? x - y : 0.0; } Index: libm/common/sf_logb.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/sf_logb.c,v retrieving revision 1.2 diff -p -u -r1.2 sf_logb.c --- libm/common/sf_logb.c 4 Apr 2001 13:30:59 -0000 1.2 +++ libm/common/sf_logb.c 20 Mar 2009 21:27:02 -0000 @@ -1,7 +1,4 @@ -/* sf_logb.c -- float version of s_logb.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - +/* 2009 for Newlib: Sun's sf_ilogb.c converted to be sf_logb.c. */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. @@ -13,24 +10,41 @@ * ==================================================== */ +/* float logb(float x) + * return the binary exponent of non-zero x + * logbf(0) = -inf, raise divide-by-zero floating point exception + * logbf(+inf|-inf) = +inf (no signal is raised) + * logbf(NaN) = NaN (no signal is raised) + * Per C99 recommendation, a NaN argument is returned unchanged. + */ + #include "fdlibm.h" +float #ifdef __STDC__ - float logbf(float x) +logbf(float x) #else - float logbf(x) - float x; +logbf(x) +float x; #endif { - __int32_t ix; - GET_FLOAT_WORD(ix,x); - ix &= 0x7fffffff; /* high |x| */ - if(FLT_UWORD_IS_ZERO(ix)) return (float)-1.0/fabsf(x); - if(!FLT_UWORD_IS_FINITE(ix)) return x*x; - if((ix>>=23)==0) /* IEEE 754 logb */ - return -126.0; - else - return (float) (ix-127); + __int32_t hx,ix; + + GET_FLOAT_WORD(hx,x); + hx &= 0x7fffffff; + if(FLT_UWORD_IS_ZERO(hx)) { + float xx; + /* arg==0: return -inf and raise divide-by-zero exception */ + SET_FLOAT_WORD(xx,hx); /* +0.0 */ + return -1./xx; /* logbf(0) = -inf */ + } + if(FLT_UWORD_IS_SUBNORMAL(hx)) { + for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1; + return (float) ix; + } + else if (FLT_UWORD_IS_INFINITE(hx)) return HUGE_VALF; /* x==+|-inf */ + else if (FLT_UWORD_IS_NAN(hx)) return x; + else return (float) ((hx>>23)-127); } #ifdef _DOUBLE_IS_32BITS Index: libm/common/s_signbit.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_signbit.c,v retrieving revision 1.1 diff -p -u -r1.1 s_signbit.c --- libm/common/s_signbit.c 7 Jun 2002 21:59:56 -0000 1.1 +++ libm/common/s_signbit.c 20 Mar 2009 21:27:02 -0000 @@ -3,6 +3,35 @@ * Permission to use, copy, modify, and distribute this software * is freely granted, provided that this notice is preserved. */ +/* +FUNCTION +<>--Does floating-point number have negative sign? + +INDEX + signbit + +ANSI_SYNOPSIS + #include + int signbit(real-floating <[x]>); + +DESCRIPTION +The <> macro determines whether the sign of its argument value is +negative. The macro reports the sign of all values, including infinities, +zeros, and NaNs. If zero is unsigned, it is treated as positive. As shown in +the synopsis, the argument is "real-floating," meaning that any of the real +floating-point types (float, double, etc.) may be given to it. + +Note that because of the possibilities of signed 0 and NaNs, the expression +"<[x]> < 0.0" does not give the same result as <> in all cases. + +RETURNS +The <> macro returns a nonzero value if and only if the sign of its +argument value is negative. + +PORTABILITY +C99, POSIX. + +*/ #include "fdlibm.h" Index: libm/common/s_log2.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_log2.c,v retrieving revision 1.1 diff -p -u -r1.1 s_log2.c --- libm/common/s_log2.c 18 Mar 2009 18:11:18 -0000 1.1 +++ libm/common/s_log2.c 20 Mar 2009 21:27:02 -0000 @@ -14,7 +14,7 @@ /* FUNCTION - <>, <>---base 2 logarithm + <>, <>--base 2 logarithm INDEX log2 INDEX @@ -25,28 +25,40 @@ ANSI_SYNOPSIS double log2(double <[x]>); float log2f(float <[x]>); -TRAD_SYNOPSIS - #include - double log2(<[x]>); - double <[x]>; - - float log2f(<[x]>); - float <[x]>; - DESCRIPTION - <> returns the base 2 logarithm of <[x]>. - - <> is identical, save that it takes and returns <> values. +The <> functions compute the base-2 logarithm of <[x]>. A domain error +occurs if the argument is less than zero. A range error occurs if the +argument is zero. + +The Newlib implementations are not full, intrinisic calculations, but +rather are derivatives based on <>. (Accuracy might be slightly off from +a direct calculation.) In addition to functions, they are also implemnted as +macros defined in math.h: +. #define log2(x) (log (x) / _M_LOG2_E) +. #define log2f(x) (logf (x) / (float) _M_LOG2_E) +To use the functions instead, just undefine the macros first. + +You can use the (non-ANSI) function <> to specify error +handling for these functions, indirectly through the respective <> +function. RETURNS - On success, <> and <> return the calculated value. - If the result underflows, the returned value is <<0>>. If the - result overflows, the returned value is <>. In - either case, <> is set to <>. +The <> functions return +@ifnottex +<)>> +@end ifnottex +@tex +$log_2(x)$ +@end tex +on success. +When <[x]> is zero, the +returned value is <<-HUGE_VAL>> and <> is set to <>. +When <[x]> is negative, the returned value is NaN (not a number) and +<> is set to <>. You can control the error behavior via +<>. PORTABILITY - <> and <> are required by ISO/IEC 9899:1999 and the - System V Interface Definition (Issue 6). +C99, POSIX, System V Interface Definition (Issue 6). */ /* Index: libm/common/sf_fma.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/sf_fma.c,v retrieving revision 1.1 diff -p -u -r1.1 sf_fma.c --- libm/common/sf_fma.c 7 Jun 2002 21:59:56 -0000 1.1 +++ libm/common/sf_fma.c 20 Mar 2009 21:27:02 -0000 @@ -15,8 +15,14 @@ float z; #endif { - /* Let the implementation handle this. */ - return (x * y) + z; + /* NOTE: The floating-point exception behavior of this is not as + * required. But since the basic function is not really done properly, + * it is not worth bothering to get the exceptions right, either. */ + /* Let the implementation handle this. */ /* <= NONSENSE! */ + /* In floating-point implementations in which double is larger than float, + * computing as double should provide the desired function. Otherwise, + * the behavior will not be as specified in the standards. */ + return (float) (((double) x * (double) y) + (double) z); } #ifdef _DOUBLE_IS_32BITS Index: libm/common/sf_remquo.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/sf_remquo.c,v retrieving revision 1.2 diff -p -u -r1.2 sf_remquo.c --- libm/common/sf_remquo.c 19 Jun 2002 17:15:47 -0000 1.2 +++ libm/common/sf_remquo.c 20 Mar 2009 21:27:02 -0000 @@ -1,50 +1,130 @@ -/* Copyright (C) 2002 by Red Hat, Incorporated. All rights reserved. +/* Adapted for Newlib, 2009. (Allow for int < 32 bits; return *quo=0 during + * errors to make test scripts easier.) */ +/* @(#)e_fmod.c 1.3 95/01/18 */ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * - * Permission to use, copy, modify, and distribute this software - * is freely granted, provided that this notice is preserved. + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== */ +#include #include "fdlibm.h" -#ifdef __STDC__ - float remquof(float x, float y, int *quo) /* wrapper remquof */ -#else - float remquof(x,y,quo) /* wrapper remquof */ - float x,y; - int *quo; +/* For quotient, return either all 31 bits that can from calculation (using + * int32_t), or as many as can fit into an int that is smaller than 32 bits. */ +#if INT_MAX > 0x7FFFFFFFL + #define QUO_MASK 0x7FFFFFFF +# else + #define QUO_MASK INT_MAX #endif -{ - int signx, signy, signres; - int wx; - int wy; - float x_over_y; - - GET_FLOAT_WORD(wx, x); - GET_FLOAT_WORD(wy, y); - - signx = (wx & 0x80000000) >> 31; - signy = (wy & 0x80000000) >> 31; - - signres = (signx ^ signy) ? -1 : 1; - - x_over_y = fabsf(x / y); - - *quo = signres * (lrintf(x_over_y) & 0x7f); - return remainderf(x,y); -} - -#ifdef _DOUBLE_IS_32BITS +static const float Zero[] = {0.0, -0.0,}; -#ifdef __STDC__ - double remquo(double x, double y, int *quo) /* wrapper remquof */ -#else - double remquo(x,y,quo) /* wrapper remquof */ - double x,y; - int *quo; -#endif +/* + * Return the IEEE remainder and set *quo to the last n bits of the + * quotient, rounded to the nearest integer. We choose n=31--if that many fit-- + * we wind up computing all the integer bits of the quotient anyway as + * a side-effect of computing the remainder by the shift and subtract + * method. In practice, this is far more bits than are needed to use + * remquo in reduction algorithms. + */ +float +remquof(float x, float y, int *quo) { - return (double) remquof((float) x, (float) y, quo); -} + int32_t n,hx,hy,hz,ix,iy,sx,i; + u_int32_t q,sxy; -#endif /* defined(_DOUBLE_IS_32BITS) */ + GET_FLOAT_WORD(hx,x); + GET_FLOAT_WORD(hy,y); + sxy = (hx ^ hy) & 0x80000000; + sx = hx&0x80000000; /* sign of x */ + hx ^=sx; /* |x| */ + hy &= 0x7fffffff; /* |y| */ + + /* purge off exception values */ + if(hy==0||hx>=0x7f800000||hy>0x7f800000) { /* y=0,NaN;or x not finite */ + *quo = 0; /* Not necessary, but return consistent value */ + return (x*y)/(x*y); + } + if(hx>31]; /* |x|=|y| return x*0*/ + } + + /* determine ix = ilogb(x) */ + if(hx<0x00800000) { /* subnormal x */ + for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1; + } else ix = (hx>>23)-127; + + /* determine iy = ilogb(y) */ + if(hy<0x00800000) { /* subnormal y */ + for (iy = -126,i=(hy<<8); i>0; i<<=1) iy -=1; + } else iy = (hy>>23)-127; + + /* set up {hx,lx}, {hy,ly} and align y to x */ + if(ix >= -126) + hx = 0x00800000|(0x007fffff&hx); + else { /* subnormal x, shift x to normal */ + n = -126-ix; + hx <<= n; + } + if(iy >= -126) + hy = 0x00800000|(0x007fffff&hy); + else { /* subnormal y, shift y to normal */ + n = -126-iy; + hy <<= n; + } + + /* fix point fmod */ + n = ix - iy; + q = 0; + while(n--) { + hz=hx-hy; + if(hz<0) hx = hx << 1; + else {hx = hz << 1; q++;} + q <<= 1; + } + hz=hx-hy; + if(hz>=0) {hx=hz;q++;} + + /* convert back to floating value and restore the sign */ + if(hx==0) { /* return sign(x)*0 */ + *quo = (sxy ? -q : q); + return Zero[(u_int32_t)sx>>31]; + } + while(hx<0x00800000) { /* normalize x */ + hx <<= 1; + iy -= 1; + } + if(iy>= -126) { /* normalize output */ + hx = ((hx-0x00800000)|((iy+127)<<23)); + } else { /* subnormal output */ + n = -126 - iy; + hx >>= n; + } +fixup: + SET_FLOAT_WORD(x,hx); + y = fabsf(y); + if (y < 0x1p-125f) { + if (x+x>y || (x+x==y && (q & 1))) { + q++; + x-=y; + } + } else if (x>0.5f*y || (x==0.5f*y && (q & 1))) { + q++; + x-=y; + } + GET_FLOAT_WORD(hx,x); + SET_FLOAT_WORD(x,hx^sx); + q &= 0x7fffffff; + *quo = (sxy ? -q : q); + return x; +} Index: libm/common/s_remquo.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_remquo.c,v retrieving revision 1.1 diff -p -u -r1.1 s_remquo.c --- libm/common/s_remquo.c 7 Jun 2002 21:59:56 -0000 1.1 +++ libm/common/s_remquo.c 20 Mar 2009 21:27:02 -0000 @@ -1,39 +1,208 @@ -/* Copyright (C) 2002 by Red Hat, Incorporated. All rights reserved. +/* Adapted for Newlib, 2009. (Allow for int < 32 bits; return *quo=0 during + * errors to make test scripts easier.) */ +/* @(#)e_fmod.c 1.3 95/01/18 */ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * - * Permission to use, copy, modify, and distribute this software - * is freely granted, provided that this notice is preserved. + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== */ +/* +FUNCTION +<>, <>--remainder and part of quotient +INDEX + remquo +INDEX + remquof + +ANSI_SYNOPSIS + #include + double remquo(double <[x]>, double <[y]>, int *<[quo]>); + float remquof(float <[x]>, float <[y]>, int *<[quo]>); + +DESCRIPTION +The <> functions compute the same remainder as the <> +functions; this value is in the range -<[y]>/2 ... +<[y]>/2. In the object +pointed to by <> they store a value whose sign is the sign of <>/<> +and whose magnitude is congruent modulo 2**n to the magnitude of the integral +quotient of <>/<>. (That is, <> is given the n lsbs of the +quotient, not counting the sign.) This implementation uses n=31 if int is 32 +bits or more, otherwise, n is 1 less than the width of int. + +For example: +. remquo(-29.0, 3.0, &<[quo]>) +returns -1.0 and sets <[quo]>=10, and +. remquo(-98307.0, 3.0, &<[quo]>) +returns -0.0 and sets <[quo]>=-32769, although for 16-bit int, <[quo]>=-1. In +the latter case, the actual quotient of -(32769=0x8001) is reduced to -1 +because of the 15-bit limitation for the quotient. + +RETURNS +When either argument is NaN, NaN is returned. If <[y]> is 0 or <[x]> is +infinite (and neither is NaN), a domain error occurs (i.e. the "invalid" +floating point exception is raised or errno is set to EDOM), and NaN is +returned. +Otherwise, the <> functions return <[x]> REM <[y]>. + +BUGS +IEEE754-2008 calls for <>(subnormal, inf) to cause the "underflow" +floating-point exception. This implementation does not. -#include "fdlibm.h" - -#ifndef _DOUBLE_IS_32BITS - -#ifdef __STDC__ - double remquo(double x, double y, int *quo) /* wrapper remquo */ -#else - double remquo(x,y,quo) /* wrapper remquo */ - double x,y; - int *quo; -#endif -{ - int signx, signy, signres; - int mswx; - int mswy; - double x_over_y; +PORTABILITY +C99, POSIX. - GET_HIGH_WORD(mswx, x); - GET_HIGH_WORD(mswy, y); +*/ - signx = (mswx & 0x80000000) >> 31; - signy = (mswy & 0x80000000) >> 31; +#include +#include +#include "fdlibm.h" - signres = (signx ^ signy) ? -1 : 1; +/* For quotient, return either all 31 bits that can from calculation (using + * int32_t), or as many as can fit into an int that is smaller than 32 bits. */ +#if INT_MAX > 0x7FFFFFFFL + #define QUO_MASK 0x7FFFFFFF +# else + #define QUO_MASK INT_MAX +#endif - x_over_y = fabs(x / y); +static const double Zero[] = {0.0, -0.0,}; - *quo = signres * (lrint(x_over_y) & 0x7f); +/* + * Return the IEEE remainder and set *quo to the last n bits of the + * quotient, rounded to the nearest integer. We choose n=31--if that many fit-- + * because we wind up computing all the integer bits of the quotient anyway as + * a side-effect of computing the remainder by the shift and subtract + * method. In practice, this is far more bits than are needed to use + * remquo in reduction algorithms. + */ +double +remquo(double x, double y, int *quo) +{ + int32_t n,hx,hy,hz,ix,iy,sx,i; + u_int32_t lx,ly,lz,q,sxy; - return remainder(x,y); + EXTRACT_WORDS(hx,lx,x); + EXTRACT_WORDS(hy,ly,y); + sxy = (hx ^ hy) & 0x80000000; + sx = hx&0x80000000; /* sign of x */ + hx ^=sx; /* |x| */ + hy &= 0x7fffffff; /* |y| */ + + /* purge off exception values */ + if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ + ((hy|((ly|-ly)>>31))>0x7ff00000)) { /* or y is NaN */ + *quo = 0; /* Not necessary, but return consistent value */ + return (x*y)/(x*y); + } + if(hx<=hy) { + if((hx>31]; /* |x|=|y| return x*0*/ + } + } + + /* determine ix = ilogb(x) */ + if(hx<0x00100000) { /* subnormal x */ + if(hx==0) { + for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; + } else { + for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; + } + } else ix = (hx>>20)-1023; + + /* determine iy = ilogb(y) */ + if(hy<0x00100000) { /* subnormal y */ + if(hy==0) { + for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; + } else { + for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; + } + } else iy = (hy>>20)-1023; + + /* set up {hx,lx}, {hy,ly} and align y to x */ + if(ix >= -1022) + hx = 0x00100000|(0x000fffff&hx); + else { /* subnormal x, shift x to normal */ + n = -1022-ix; + if(n<=31) { + hx = (hx<>(32-n)); + lx <<= n; + } else { + hx = lx<<(n-32); + lx = 0; + } + } + if(iy >= -1022) + hy = 0x00100000|(0x000fffff&hy); + else { /* subnormal y, shift y to normal */ + n = -1022-iy; + if(n<=31) { + hy = (hy<>(32-n)); + ly <<= n; + } else { + hy = ly<<(n-32); + ly = 0; + } + } + + /* fix point fmod */ + n = ix - iy; + q = 0; + while(n--) { + hz=hx-hy;lz=lx-ly; if(lx>31); lx = lx+lx;} + else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;} + q <<= 1; + } + hz=hx-hy;lz=lx-ly; if(lx=0) {hx=hz;lx=lz;q++;} + + /* convert back to floating value and restore the sign */ + if((hx|lx)==0) { /* return sign(x)*0 */ + q &= QUO_MASK; + *quo = (sxy ? -q : q); + return Zero[(u_int32_t)sx>>31]; + } + while(hx<0x00100000) { /* normalize x */ + hx = hx+hx+(lx>>31); lx = lx+lx; + iy -= 1; + } + if(iy>= -1022) { /* normalize output */ + hx = ((hx-0x00100000)|((iy+1023)<<20)); + } else { /* subnormal output */ + n = -1022 - iy; + if(n<=20) { + lx = (lx>>n)|((u_int32_t)hx<<(32-n)); + hx >>= n; + } else if (n<=31) { + lx = (hx<<(32-n))|(lx>>n); hx = sx; + } else { + lx = hx>>(n-32); hx = sx; + } + } +fixup: + INSERT_WORDS(x,hx,lx); + y = fabs(y); + if (y < 0x1p-1021) { + if (x+x>y || (x+x==y && (q & 1))) { + q++; + x-=y; + } + } else if (x>0.5*y || (x==0.5*y && (q & 1))) { + q++; + x-=y; + } + GET_HIGH_WORD(hx,x); + SET_HIGH_WORD(x,hx^sx); + q &= QUO_MASK; + *quo = (sxy ? -q : q); + return x; } - -#endif /* defined(_DOUBLE_IS_32BITS) */ Index: libm/common/s_fmax.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_fmax.c,v retrieving revision 1.2 diff -p -u -r1.2 s_fmax.c --- libm/common/s_fmax.c 10 Aug 2005 21:02:28 -0000 1.2 +++ libm/common/s_fmax.c 20 Mar 2009 21:27:02 -0000 @@ -3,6 +3,31 @@ * Permission to use, copy, modify, and distribute this software * is freely granted, provided that this notice is preserved. */ +/* +FUNCTION +<>, <>--maximum +INDEX + fmax +INDEX + fmaxf + +ANSI_SYNOPSIS + #include + double fmax(double <[x]>, double <[y]>); + float fmaxf(float <[x]>, float <[y]>); + +DESCRIPTION +The <> functions determine the maximum numeric value of their arguments. +NaN arguments are treated as missing data: if one argument is a NaN and the +other numeric, then the <> functions choose the numeric value. + +RETURNS +The <> functions return the maximum numeric value of their arguments. + +PORTABILITY +ANSI C, POSIX. + +*/ #include "fdlibm.h" Index: libm/common/s_fmin.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_fmin.c,v retrieving revision 1.2 diff -p -u -r1.2 s_fmin.c --- libm/common/s_fmin.c 10 Aug 2005 21:02:28 -0000 1.2 +++ libm/common/s_fmin.c 20 Mar 2009 21:27:02 -0000 @@ -3,6 +3,31 @@ * Permission to use, copy, modify, and distribute this software * is freely granted, provided that this notice is preserved. */ +/* +FUNCTION +<>, <>--minimum +INDEX + fmin +INDEX + fminf + +ANSI_SYNOPSIS + #include + double fmin(double <[x]>, double <[y]>); + float fminf(float <[x]>, float <[y]>); + +DESCRIPTION +The <> functions determine the minimum numeric value of their arguments. +NaN arguments are treated as missing data: if one argument is a NaN and the +other numeric, then the <> functions choose the numeric value. + +RETURNS +The <> functions return the minimum numeric value of their arguments. + +PORTABILITY +ANSI C, POSIX. + +*/ #include "fdlibm.h" Index: libm/common/s_fdim.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_fdim.c,v retrieving revision 1.1 diff -p -u -r1.1 s_fdim.c --- libm/common/s_fdim.c 7 Jun 2002 21:59:56 -0000 1.1 +++ libm/common/s_fdim.c 20 Mar 2009 21:27:02 -0000 @@ -3,6 +3,39 @@ * Permission to use, copy, modify, and distribute this software * is freely granted, provided that this notice is preserved. */ +/* +FUNCTION +<>, <>--positive difference +INDEX + fdim +INDEX + fdimf + +ANSI_SYNOPSIS + #include + double fdim(double <[x]>, double <[y]>); + float fdimf(float <[x]>, float <[y]>); + +DESCRIPTION +The <> functions determine the positive difference between their +arguments, returning: +. <[x]> - <[y]> if <[x]> > <[y]>, or + @ifnottex +. +0 if <[x]> <= <[y]>, or + @end ifnottex + @tex +. +0 if <[x]> $\leq$ <[y]>, or + @end tex +. NAN if either argument is NAN. +A range error may occur. + +RETURNS +The <> functions return the positive difference value. + +PORTABILITY +ANSI C, POSIX. + +*/ #include "fdlibm.h" @@ -17,7 +50,9 @@ #endif { int c = __fpclassifyd(x); - if (c == FP_NAN || c == FP_INFINITE) + if (c == FP_NAN) return(x); + if (__fpclassifyd(y) == FP_NAN) return(y); + if (c == FP_INFINITE) return HUGE_VAL; return x > y ? x - y : 0.0; Index: libm/common/s_trunc.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_trunc.c,v retrieving revision 1.1 diff -p -u -r1.1 s_trunc.c --- libm/common/s_trunc.c 7 Jun 2002 21:59:56 -0000 1.1 +++ libm/common/s_trunc.c 20 Mar 2009 21:27:02 -0000 @@ -8,6 +8,35 @@ * is preserved. * ==================================================== */ +/* +FUNCTION +<>, <>--round to integer, towards zero +INDEX + trunc +INDEX + truncf + +ANSI_SYNOPSIS + #include + double trunc(double <[x]>); + float truncf(float <[x]>); + +DESCRIPTION + The <> functions round their argument to the integer value, in + floating format, nearest to but no larger in magnitude than the + argument, regardless of the current rounding direction. (While the + "inexact" floating-point exception behavior is unspecified by the C + standard, the <> functions are written so that "inexact" is not + raised if the result does not equal the argument, which behavior is as + recommended by IEEE 754 for its related functions.) + +RETURNS +<[x]> truncated to an integral value. + +PORTABILITY +ANSI C, POSIX + +*/ #include "fdlibm.h" Index: libm/common/s_rint.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_rint.c,v retrieving revision 1.1.1.1 diff -p -u -r1.1.1.1 s_rint.c --- libm/common/s_rint.c 17 Feb 2000 19:39:51 -0000 1.1.1.1 +++ libm/common/s_rint.c 20 Mar 2009 21:27:02 -0000 @@ -10,6 +10,40 @@ * is preserved. * ==================================================== */ +/* +FUNCTION +<>, <>--round to integer +INDEX + rint +INDEX + rintf + +ANSI_SYNOPSIS + #include + double rint(double <[x]>); + float rintf(float <[x]>); + +DESCRIPTION + The <> functions round their argument to an integer value in + floating-point format, using the current rounding direction. They + raise the "inexact" floating-point exception if the result differs + in value from the argument. See the <> functions for the + same function with the "inexact" floating-point exception never being + raised. Newlib does not directly support floating-point exceptions. + The <> functions are written so that the "inexact" exception is + raised in hardware implementations that support it, even though Newlib + does not provide access. + +RETURNS +<[x]> rounded to an integral value, using the current rounding direction. + +PORTABILITY +ANSI C, POSIX + +SEEALSO +<>, <> + +*/ /* * rint(x) Index: libm/common/s_round.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_round.c,v retrieving revision 1.1 diff -p -u -r1.1 s_round.c --- libm/common/s_round.c 7 Jun 2002 21:59:56 -0000 1.1 +++ libm/common/s_round.c 20 Mar 2009 21:27:02 -0000 @@ -8,6 +8,38 @@ * is preserved. * ==================================================== */ +/* +FUNCTION +<>, <>--round to integer, to nearest +INDEX + round +INDEX + roundf + +ANSI_SYNOPSIS + #include + double round(double <[x]>); + float roundf(float <[x]>); + +DESCRIPTION + The <> functions round their argument to the nearest integer + value in floating-point format, rounding halfway cases away from zero, + regardless of the current rounding direction. (While the "inexact" + floating-point exception behavior is unspecified by the C standard, the + <> functions are written so that "inexact" is not raised if the + result does not equal the argument, which behavior is as recommended by + IEEE 754 for its related functions.) + +RETURNS +<[x]> rounded to an integral value. + +PORTABILITY +ANSI C, POSIX + +SEEALSO +<>, <> + +*/ #include "fdlibm.h" Index: libm/common/s_scalbn.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_scalbn.c,v retrieving revision 1.3 diff -p -u -r1.3 s_scalbn.c --- libm/common/s_scalbn.c 20 Oct 2003 18:46:37 -0000 1.3 +++ libm/common/s_scalbn.c 20 Mar 2009 21:27:02 -0000 @@ -13,38 +13,44 @@ /* FUNCTION -<>, <>---scale by power of two +<>, <>, <>, <>--scale by power of FLT_RADIX (=2) INDEX scalbn INDEX scalbnf +INDEX + scalbln +INDEX + scalblnf ANSI_SYNOPSIS #include - double scalbn(double <[x]>, int <[y]>); - float scalbnf(float <[x]>, int <[y]>); - -TRAD_SYNOPSIS - #include - double scalbn(<[x]>,<[y]>) - double <[x]>; - int <[y]>; - float scalbnf(<[x]>,<[y]>) - float <[x]>; - int <[y]>; + double scalbn(double <[x]>, int <[n]>); + float scalbnf(float <[x]>, int <[n]>); + double scalbln(double <[x]>, long int <[n]>); + float scalblnf(float <[x]>, long int <[n]>); DESCRIPTION -<> and <> scale <[x]> by <[n]>, returning <[x]> times -2 to the power <[n]>. The result is computed by manipulating the -exponent, rather than by actually performing an exponentiation or -multiplication. +The <> and <> functions compute + @ifnottex + <[x]> times FLT_RADIX to the power <[n]>. + @end ifnottex + @tex + $x \cdot FLT\_RADIX^n$. + @end tex +efficiently. The result is computed by manipulating the exponent, rather than +by actually performing an exponentiation or multiplication. In this +floating-point implementation FLT_RADIX=2, which makes the <> +functions equivalent to the <> functions. RETURNS -<[x]> times 2 to the power <[n]>. +<[x]> times 2 to the power <[n]>. A range error may occur. PORTABILITY -Neither <> nor <> is required by ANSI C or by the System V -Interface Definition (Issue 2). +ANSI C, POSIX + +SEEALSO +<> */ Index: libm/common/s_infinity.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_infinity.c,v retrieving revision 1.2 diff -p -u -r1.2 s_infinity.c --- libm/common/s_infinity.c 20 Oct 2003 18:46:37 -0000 1.2 +++ libm/common/s_infinity.c 20 Mar 2009 21:27:02 -0000 @@ -5,7 +5,7 @@ /* FUNCTION - <>, <>---representation of infinity + <>, <>--representation of infinity INDEX infinity @@ -17,17 +17,18 @@ ANSI_SYNOPSIS double infinity(void); float infinityf(void); -TRAD_SYNOPSIS - #include - double infinity(); - float infinityf(); - - DESCRIPTION <> and <> return the special number IEEE infinity in double- and single-precision arithmetic respectively. +PORTABILITY +<> and <> are neither standard C nor POSIX. C and +POSIX require macros HUGE_VAL and HUGE_VALF to be defined in math.h, which +Newlib defines to be infinities corresponding to these archaic infinity() +and infinityf() functions in floating-point implementations which do have +infinities. + QUICKREF infinity - pure Index: libm/common/s_lround.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_lround.c,v retrieving revision 1.1 diff -p -u -r1.1 s_lround.c --- libm/common/s_lround.c 7 Jun 2002 21:59:56 -0000 1.1 +++ libm/common/s_lround.c 20 Mar 2009 21:27:02 -0000 @@ -8,6 +8,44 @@ * is preserved. * ==================================================== */ +/* +FUNCTION +<>, <>, <>, <>--round to integer, to nearest +INDEX + lround +INDEX + lroundf +INDEX + llround +INDEX + llroundf + +ANSI_SYNOPSIS + #include + long int lround(double <[x]>); + long int lroundf(float <[x]>); + long long int llround(double <[x]>); + long long int llroundf(float <[x]>); + +DESCRIPTION + The <> and <> functions round their argument to the + nearest integer value, rounding halfway cases away from zero, regardless + of the current rounding direction. If the rounded value is outside the + range of the return type, the numeric result is unspecified (depending + upon the floating-point implementation, not the library). A range + error may occur if the magnitude of x is too large. + +RETURNS +<[x]> rounded to an integral value as an integer. + +SEEALSO +See the <> functions for the return being the same floating-point type +as the argument. <>, <>. + +PORTABILITY +ANSI C, POSIX + +*/ #include "fdlibm.h" Index: libm/common/s_lrint.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_lrint.c,v retrieving revision 1.3 diff -p -u -r1.3 s_lrint.c --- libm/common/s_lrint.c 8 Sep 2005 21:03:54 -0000 1.3 +++ libm/common/s_lrint.c 20 Mar 2009 21:27:02 -0000 @@ -10,6 +10,44 @@ * is preserved. * ==================================================== */ +/* +FUNCTION +<>, <>, <>, <>--round to integer +INDEX + lrint +INDEX + lrintf +INDEX + llrint +INDEX + llrintf + +ANSI_SYNOPSIS + #include + long int lrint(double <[x]>); + long int lrintf(float <[x]>); + long long int llrint(double <[x]>); + long long int llrintf(float <[x]>); + +DESCRIPTION +The <> and <> functions round their argument to the nearest +integer value, using the current rounding direction. If the rounded value is +outside the range of the return type, the numeric result is unspecified. A +range error may occur if the magnitude of <[x]> is too large. +The "inexact" floating-point exception is raised in implementations that +support it when the result differs in value from the argument (i.e., when +a fraction actually has been truncated). + +RETURNS +<[x]> rounded to an integral value, using the current rounding direction. + +SEEALSO +<> + +PORTABILITY +ANSI C, POSIX + +*/ /* * lrint(x) Index: libm/common/s_fma.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_fma.c,v retrieving revision 1.1 diff -p -u -r1.1 s_fma.c --- libm/common/s_fma.c 7 Jun 2002 21:59:56 -0000 1.1 +++ libm/common/s_fma.c 20 Mar 2009 21:27:02 -0000 @@ -1,3 +1,41 @@ +/* +FUNCTION +<>, <>--floating multiply add +INDEX + fma +INDEX + fmaf + +ANSI_SYNOPSIS + #include + double fma(double <[x]>, double <[y]>, double <[z]>); + float fmaf(float <[x]>, float <[y]>, float <[z]>); + +DESCRIPTION +The <> functions compute (<[x]> * <[y]>) + <[z]>, rounded as one ternary +operation: they compute the value (as if) to infinite precision and round once +to the result format, according to the rounding mode characterized by the value +of FLT_ROUNDS. That is, they are supposed to do this: see below. + +RETURNS +The <> functions return (<[x]> * <[y]>) + <[z]>, rounded as one ternary +operation. + +BUGS +This implementation does not provide the function that it should, purely +returning "(<[x]> * <[y]>) + <[z]>;" with no attempt at all to provide the +simulated infinite precision intermediates which are required. DO NOT USE THEM. + +If float has enough more precision than double, then <> should provide +the expected numeric results, as it does use double for the calculation. But +since this is not the case for all platforms, this manual cannot determine +if it is so for your case. + +PORTABILITY +ANSI C, POSIX. + +*/ + #include "fdlibm.h" #ifndef _DOUBLE_IS_32BITS Index: libm/common/s_logb.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_logb.c,v retrieving revision 1.1.1.1 diff -p -u -r1.1.1.1 s_logb.c --- libm/common/s_logb.c 17 Feb 2000 19:39:51 -0000 1.1.1.1 +++ libm/common/s_logb.c 20 Mar 2009 21:27:02 -0000 @@ -1,5 +1,5 @@ - -/* @(#)s_logb.c 5.1 93/09/24 */ +/* 2009 for Newlib: Sun's s_ilogb.c converted to be s_logb.c. */ +/* @(#)s_ilogb.c 5.1 93/09/24 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. @@ -10,33 +10,101 @@ * is preserved. * ==================================================== */ - /* - * double logb(x) - * IEEE 754 logb. Included to pass IEEE test suite. Not recommend. - * Use ilogb instead. +FUNCTION + <>, <>--get exponent of floating-point number +INDEX + logb +INDEX + logbf + +ANSI_SYNOPSIS + #include + double logb(double <[x]>); + float logbf(float <[x]>); + +DESCRIPTION +The <> functions extract the exponent of <[x]>, as a signed integer value +in floating-point format. If <[x]> is subnormal it is treated as though it were +normalized; thus, for positive finite <[x]>, +@ifnottex +1 <= (<[x]> * FLT_RADIX to the power (-logb(<[x]>))) < FLT_RADIX. +@end ifnottex +@tex +$1 \leq ( x \cdot FLT\_RADIX ^ {-logb(x)} ) < FLT\_RADIX$. +@end tex +A domain error may occur if the argument is zero. +In this floating-point implementation, FLT_RADIX is 2. Which also means +that for finite <[x]>, <>(<[x]>) = <>(<>(<>(<[x]>))). + +All nonzero, normal numbers can be described as +@ifnottex +<[m]> * 2**<[p]>, where 1.0 <= <[m]> < 2.0. +@end ifnottex +@tex +$m \cdot 2^p$, where $1.0 \leq m < 2.0$. +@end tex +The <> functions examine the argument <[x]>, and return <[p]>. +The <> functions are similar to the <> functions, but +returning <[m]> adjusted to the interval [.5, 1) or 0, and <[p]>+1. + +RETURNS +@comment Formatting note: "$@" forces a new line +When <[x]> is:@* ++inf or -inf, +inf is returned;@* +NaN, NaN is returned;@* +0, -inf is returned, and the divide-by-zero exception is raised;@* +otherwise, the <> functions return the signed exponent of <[x]>. + +PORTABILITY +ANSI C, POSIX + +SEEALSO +frexp, ilogb +*/ + +/* double logb(double x) + * return the binary exponent of non-zero x + * logb(0) = -inf, raise divide-by-zero floating point exception + * logb(+inf|-inf) = +inf (no signal is raised) + * logb(NaN) = NaN (no signal is raised) + * Per C99 recommendation, a NaN argument is returned unchanged. */ #include "fdlibm.h" #ifndef _DOUBLE_IS_32BITS +double #ifdef __STDC__ - double logb(double x) +logb(double x) #else - double logb(x) - double x; +logb(x) +double x; #endif { - __int32_t lx,ix; - EXTRACT_WORDS(ix,lx,x); - ix &= 0x7fffffff; /* high |x| */ - if((ix|lx)==0) return -1.0/fabs(x); - if(ix>=0x7ff00000) return x*x; - if((ix>>=20)==0) /* IEEE 754 logb */ - return -1022.0; - else - return (double) (ix-1023); + __int32_t hx,lx,ix; + + EXTRACT_WORDS(hx,lx,x); + hx &= 0x7fffffff; /* high |x| */ + if(hx<0x00100000) { /* 0 or subnormal */ + if((hx|lx)==0) { + double xx; + /* arg==0: return -inf and raise divide-by-zero exception */ + INSERT_WORDS(xx,hx,lx); /* +0.0 */ + return -1./xx; /* logb(0) = -inf */ + } + else /* subnormal x */ + if(hx==0) { + for (ix = -1043; lx>0; lx<<=1) ix -=1; + } else { + for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1; + } + return (double) ix; + } + else if (hx<0x7ff00000) return (hx>>20)-1023; /* normal # */ + else if (hx>0x7ff00000 || lx) return x; /* x==NaN */ + else return HUGE_VAL; /* x==inf (+ or -) */ } #endif /* _DOUBLE_IS_32BITS */ Index: libm/common/s_isnan.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/s_isnan.c,v retrieving revision 1.1 diff -p -u -r1.1 s_isnan.c --- libm/common/s_isnan.c 27 Feb 2006 23:51:28 -0000 1.1 +++ libm/common/s_isnan.c 20 Mar 2009 21:27:02 -0000 @@ -13,8 +13,21 @@ /* FUNCTION - <>, <>, <>, <>, <>, <>---test for exceptional numbers +<>, <>, <>, <>, and <>--floating-point classification macros; <>, <>, <>, <>, <>, <>--test for exceptional numbers +@c C99 (start +INDEX + fpclassify +INDEX + isfinite +INDEX + isinf +INDEX + isnan +INDEX + isnormal +@c C99 end) +@c SUSv2 (start INDEX isnan INDEX @@ -28,8 +41,18 @@ INDEX isinff INDEX finitef +@c SUSv2 end) ANSI_SYNOPSIS + [C99 standard macros:] + #include + int fpclassify(real-floating <[x]>); + int isfinite(real-floating <[x]>); + int isinf(real-floating <[x]>); + int isnan(real-floating <[x]>); + int isnormal(real-floating <[x]>); + + [Archaic SUSv2 functions:] #include int isnan(double <[arg]>); int isinf(double <[arg]>); @@ -38,30 +61,64 @@ ANSI_SYNOPSIS int isinff(float <[arg]>); int finitef(float <[arg]>); -TRAD_SYNOPSIS - #include - int isnan(<[arg]>) - double <[arg]>; - int isinf(<[arg]>) - double <[arg]>; - int finite(<[arg]>); - double <[arg]>; - int isnanf(<[arg]>); - float <[arg]>; - int isinff(<[arg]>); - float <[arg]>; - int finitef(<[arg]>); - float <[arg]>; - - DESCRIPTION - These functions provide information on the floating-point +<>, <>, <>, <>, and <> are macros +defined for use in classifying floating-point numbers. This is a help because +of special "values" like NaN and infinities. In the synopses shown, +"real-floating" indicates that the argument is an expression of real floating +type. These function-like macros are C99 and POSIX-compliant, and should be +used instead of the now-archaic SUSv2 functions. + +The <> macro classifies its argument value as NaN, infinite, normal, +subnormal, zero, or into another implementation-defined category. First, an +argument represented in a format wider than its semantic type is converted to +its semantic type. Then classification is based on the type of the argument. +The <> macro returns the value of the number classification macro +appropriate to the value of its argument: + +o+ +o FP_INFINITE + <[x]> is either plus or minus infinity; +o FP_NAN + <[x]> is "Not A Number" (plus or minus); +o FP_NORMAL + <[x]> is a "normal" number (i.e. is none of the other special forms); +o FP_SUBNORMAL + <[x]> is too small be stored as a regular normalized number (i.e. loss of precision is likely); or +o FP_ZERO + <[x]> is 0 (either plus or minus). +o- + +The "<>" set of macros provide a useful set of shorthand ways for +classifying floating-point numbers, providing the following equivalent +relations: + +o+ +o <>(<[x]>) +returns non-zero if <[x]> is finite. (It is equivalent to +(<>(<[x]>) != FP_INFINITE && <>(<[x]>) != FP_NAN).) + +o <>(<[x]>) +returns non-zero if <[x]> is infinite. (It is equivalent to +(<>(<[x]>) == FP_INFINITE).) + +o <>(<[x]>) +returns non-zero if <[x]> is NaN. (It is equivalent to +(<>(<[x]>) == FP_NAN).) + +o <>(<[x]>) +returns non-zero if <[x]> is normal. (It is equivalent to +(<>(<[x]>) == FP_NORMAL).) +o- + + The archaic SUSv2 functions provide information on the floating-point argument supplied. - There are five major number formats: + There are five major number formats ("exponent" referring to the + biased exponent in the binary-encoded number): o+ o zero - A number which contains all zero bits. + A number which contains all zero bits, excluding the sign bit. o subnormal A number with a zero exponent but a nonzero fraction. o normal @@ -85,7 +142,21 @@ DESCRIPTION and <> are macros that operate on multiple types of floating-point. The SUSv2 standard declares <> as a function taking double. Newlib has decided to declare - them both as macros in math.h and as functions in ieeefp.h. + them both as macros in math.h and as functions in ieeefp.h to + maintain backward compatibility. + +RETURNS +@comment Formatting note: "$@" forces a new line +The fpclassify macro returns the value corresponding to the appropriate FP_ macro.@* +The isfinite macro returns nonzero if <[x]> is finite, else 0.@* +The isinf macro returns nonzero if <[x]> is infinite, else 0.@* +The isnan macro returns nonzero if <[x]> is an NaN, else 0.@* +The isnormal macro returns nonzero if <[x]> has a normal value, else 0. + +PORTABILITY +math.h macros are C99, POSIX. + +ieeefp.h funtions are outdated and should be avoided. QUICKREF isnan - pure