* libc/include/machine/ieeefp.h: Comment about new configuration

macros _FLT_LARGEST_EXPONENT_IS_NORMAL and _FLT_NO_DENORMALS.
	* libm/common/fdlib.h: Define new macros for testing floats.
	* libm/common/sf_*: Use them.
	* libm/math/ef_*: Likewise.
	* libm/math/sf_*: Likewise.
This commit is contained in:
Richard Sandiford 2001-04-04 13:33:01 +00:00
parent 51fc3813e9
commit 16740220a2
41 changed files with 306 additions and 166 deletions

View File

@ -1,3 +1,12 @@
2001-04-04 Richard Sandiford <rsandifo@redhat.com>
* libc/include/machine/ieeefp.h: Comment about new configuration
macros _FLT_LARGEST_EXPONENT_IS_NORMAL and _FLT_NO_DENORMALS.
* libm/common/fdlib.h: Define new macros for testing floats.
* libm/common/sf_*: Use them.
* libm/math/ef_*: Likewise.
* libm/math/sf_*: Likewise.
2001-03-29 Jeff Johnston <jjohnstn@redhat.com> 2001-03-29 Jeff Johnston <jjohnstn@redhat.com>
* libc/sys/arm/setjmp.S: Added .code 16 specifier for thumb-mode * libc/sys/arm/setjmp.S: Added .code 16 specifier for thumb-mode

View File

@ -1,6 +1,24 @@
#ifndef __IEEE_BIG_ENDIAN #ifndef __IEEE_BIG_ENDIAN
#ifndef __IEEE_LITTLE_ENDIAN #ifndef __IEEE_LITTLE_ENDIAN
/* This file can define macros to choose variations of the IEEE float
format:
_FLT_LARGEST_EXPONENT_IS_NORMAL
Defined if the float format uses the largest exponent for finite
numbers rather than NaN and infinity representations. Such a
format cannot represent NaNs or infinities at all, but it's FLT_MAX
is twice the IEEE value.
_FLT_NO_DENORMALS
Defined if the float format does not support IEEE denormals. Every
float with a zero exponent is taken to be a zero representation.
??? At the moment, there are no equivalent macros for doubles and
the macros are not fully supported by --enable-newlib-hw-fp. */
#if defined(__arm__) || defined(__thumb__) #if defined(__arm__) || defined(__thumb__)
/* ARM always has big-endian words. Within those words the byte ordering /* ARM always has big-endian words. Within those words the byte ordering
appears to be big or little endian. Newlib doesn't seem to care about appears to be big or little endian. Newlib doesn't seem to care about

View File

@ -18,14 +18,116 @@
/* CYGNUS LOCAL: Default to XOPEN_MODE. */ /* CYGNUS LOCAL: Default to XOPEN_MODE. */
#define _XOPEN_MODE #define _XOPEN_MODE
/* Most routines need to check whether a float is finite, infinite, or not a
number, and many need to know whether the result of an operation will
overflow. These conditions depend on whether the largest exponent is
used for NaNs & infinities, or whether it's used for finite numbers. The
macros below wrap up that kind of information:
FLT_UWORD_IS_FINITE(X)
True if a positive float with bitmask X is finite.
FLT_UWORD_IS_NAN(X)
True if a positive float with bitmask X is not a number.
FLT_UWORD_IS_INFINITE(X)
True if a positive float with bitmask X is +infinity.
FLT_UWORD_MAX
The bitmask of FLT_MAX.
FLT_UWORD_HALF_MAX
The bitmask of FLT_MAX/2.
FLT_UWORD_EXP_MAX
The bitmask of the largest finite exponent (129 if the largest
exponent is used for finite numbers, 128 otherwise).
FLT_UWORD_LOG_MAX
The bitmask of log(FLT_MAX), rounded down. This value is the largest
input that can be passed to exp() without producing overflow.
FLT_UWORD_LOG_2MAX
The bitmask of log(2*FLT_MAX), rounded down. This value is the
largest input than can be passed to cosh() without producing
overflow.
FLT_LARGEST_EXP
The largest biased exponent that can be used for finite numbers
(255 if the largest exponent is used for finite numbers, 254
otherwise) */
#ifdef _FLT_LARGEST_EXPONENT_IS_NORMAL
#define FLT_UWORD_IS_FINITE(x) 1
#define FLT_UWORD_IS_NAN(x) 0
#define FLT_UWORD_IS_INFINITE(x) 0
#define FLT_UWORD_MAX 0x7fffffff
#define FLT_UWORD_EXP_MAX 0x43010000
#define FLT_UWORD_LOG_MAX 0x42b2d4fc
#define FLT_UWORD_LOG_2MAX 0x42b437e0
#define HUGE ((float)0X1.FFFFFEP128)
#else
#define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L)
#define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L)
#define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L)
#define FLT_UWORD_MAX 0x7f7fffff
#define FLT_UWORD_EXP_MAX 0x43000000
#define FLT_UWORD_LOG_MAX 0x42b17217
#define FLT_UWORD_LOG_2MAX 0x42b2d4fc
#define HUGE ((float)3.40282346638528860e+38)
#endif
#define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1<<23))
#define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23)
/* Many routines check for zero and subnormal numbers. Such things depend
on whether the target supports denormals or not:
FLT_UWORD_IS_ZERO(X)
True if a positive float with bitmask X is +0. Without denormals,
any float with a zero exponent is a +0 representation. With
denormals, the only +0 representation is a 0 bitmask.
FLT_UWORD_IS_SUBNORMAL(X)
True if a non-zero positive float with bitmask X is subnormal.
(Routines should check for zeros first.)
FLT_UWORD_MIN
The bitmask of the smallest float above +0. Call this number
REAL_FLT_MIN...
FLT_UWORD_EXP_MIN
The bitmask of the float representation of REAL_FLT_MIN's exponent.
FLT_UWORD_LOG_MIN
The bitmask of |log(REAL_FLT_MIN)|, rounding down.
FLT_SMALLEST_EXP
REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported,
-22 if they are).
*/
#ifdef _FLT_NO_DENORMALS
#define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L)
#define FLT_UWORD_IS_SUBNORMAL(x) 0
#define FLT_UWORD_MIN 0x00800000
#define FLT_UWORD_EXP_MIN 0x42fc0000
#define FLT_UWORD_LOG_MIN 0x42aeac50
#define FLT_SMALLEST_EXP 1
#else
#define FLT_UWORD_IS_ZERO(x) ((x)==0)
#define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L)
#define FLT_UWORD_MIN 0x00000001
#define FLT_UWORD_EXP_MIN 0x43160000
#define FLT_UWORD_LOG_MIN 0x42cff1b5
#define FLT_SMALLEST_EXP -22
#endif
#ifdef __STDC__ #ifdef __STDC__
#define __P(p) p #define __P(p) p
#else #else
#define __P(p) () #define __P(p) ()
#endif #endif
#define HUGE ((float)3.40282346638528860e+38)
/* /*
* set X_TLOSS = pi*2**52, which is possibly defined in <values.h> * set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
* (one may replace the following line by "#include <values.h>") * (one may replace the following line by "#include <values.h>")

View File

@ -53,13 +53,14 @@ G = 3.5714286566e-01; /* 5/14 = 0x3eb6db6e */
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
sign=hx&0x80000000; /* sign= sign(x) */ sign=hx&0x80000000; /* sign= sign(x) */
hx ^=sign; hx ^=sign;
if(hx>=0x7f800000) return(x+x); /* cbrt(NaN,INF) is itself */ if(!FLT_UWORD_IS_FINITE(hx))
if(hx==0) return(x+x); /* cbrt(NaN,INF) is itself */
return(x); /* cbrt(0) is itself */ if(FLT_UWORD_IS_ZERO(hx))
return(x); /* cbrt(0) is itself */
SET_FLOAT_WORD(x,hx); /* x <- |x| */ SET_FLOAT_WORD(x,hx); /* x <- |x| */
/* rough cbrt to 5 bits */ /* rough cbrt to 5 bits */
if(hx<0x00800000) /* subnormal number */ if(FLT_UWORD_IS_SUBNORMAL(hx)) /* subnormal number */
{SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */ {SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */
t*=x; GET_FLOAT_WORD(high,t); SET_FLOAT_WORD(t,high/3+B2); t*=x; GET_FLOAT_WORD(high,t); SET_FLOAT_WORD(t,high/3+B2);
} }

View File

@ -27,7 +27,6 @@ static float
one = 1.0, one = 1.0,
huge = 1.0e+30, huge = 1.0e+30,
tiny = 1.0e-30, tiny = 1.0e-30,
o_threshold = 8.8721679688e+01,/* 0x42b17180 */
ln2_hi = 6.9313812256e-01,/* 0x3f317180 */ ln2_hi = 6.9313812256e-01,/* 0x3f317180 */
ln2_lo = 9.0580006145e-06,/* 0x3717f7d1 */ ln2_lo = 9.0580006145e-06,/* 0x3717f7d1 */
invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */ invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */
@ -56,13 +55,12 @@ Q5 = -2.0109921195e-07; /* 0xb457edbb */
/* filter out huge and non-finite argument */ /* filter out huge and non-finite argument */
if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */ if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */
if(hx >= 0x42b17218) { /* if |x|>=88.721... */ if(FLT_UWORD_IS_NAN(hx))
if(hx>0x7f800000) return x+x;
return x+x; /* NaN */ if(FLT_UWORD_IS_INFINITE(hx))
if(hx==0x7f800000) return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ if(xsb == 0 && hx > FLT_UWORD_LOG_MAX) /* if x>=o_threshold */
if(x > o_threshold) return huge*huge; /* overflow */ return huge*huge; /* overflow */
}
if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */ if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */
if(x+tiny<(float)0.0) /* raise inexact */ if(x+tiny<(float)0.0) /* raise inexact */
return tiny-one; /* return -1 */ return tiny-one; /* return -1 */

View File

@ -29,7 +29,8 @@
{ {
__int32_t ix; __int32_t ix;
GET_FLOAT_WORD(ix,x); GET_FLOAT_WORD(ix,x);
return (int)((__uint32_t)((ix&0x7fffffff)-0x7f800000)>>31); ix &= 0x7fffffff;
return (FLT_UWORD_IS_FINITE(ix));
} }
#ifdef _DOUBLE_IS_32BITS #ifdef _DOUBLE_IS_32BITS

View File

@ -27,15 +27,14 @@
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
hx &= 0x7fffffff; hx &= 0x7fffffff;
if(hx<0x00800000) { if(FLT_UWORD_IS_ZERO(hx))
if(hx==0) return - INT_MAX; /* ilogb(0) = 0x80000001 */
return - INT_MAX; /* ilogb(0) = 0x80000001 */ if(FLT_UWORD_IS_SUBNORMAL(hx)) {
else /* subnormal x */ for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1;
for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1;
return ix; return ix;
} }
else if (hx<0x7f800000) return (hx>>23)-127; else if (!FLT_UWORD_IS_FINITE(hx)) return INT_MAX;
else return INT_MAX; else return (hx>>23)-127;
} }
#ifdef _DOUBLE_IS_32BITS #ifdef _DOUBLE_IS_32BITS

View File

@ -51,6 +51,7 @@ static float zero = 0.0;
ax = hx&0x7fffffff; ax = hx&0x7fffffff;
k = 1; k = 1;
if (!FLT_UWORD_IS_FINITE(hx)) return x+x;
if (hx < 0x3ed413d7) { /* x < 0.41422 */ if (hx < 0x3ed413d7) { /* x < 0.41422 */
if(ax>=0x3f800000) { /* x <= -1.0 */ if(ax>=0x3f800000) { /* x <= -1.0 */
if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=+inf */ if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=+inf */
@ -65,8 +66,7 @@ static float zero = 0.0;
} }
if(hx>0||hx<=((__int32_t)0xbe95f61f)) { if(hx>0||hx<=((__int32_t)0xbe95f61f)) {
k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */ k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */
} }
if (hx >= 0x7f800000) return x+x;
if(k!=0) { if(k!=0) {
if(hx<0x5a000000) { if(hx<0x5a000000) {
u = (float)1.0+x; u = (float)1.0+x;

View File

@ -25,8 +25,8 @@
__int32_t ix; __int32_t ix;
GET_FLOAT_WORD(ix,x); GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff; /* high |x| */ ix &= 0x7fffffff; /* high |x| */
if(ix==0) return (float)-1.0/fabsf(x); if(FLT_UWORD_IS_ZERO(ix)) return (float)-1.0/fabsf(x);
if(ix>=0x7f800000) return x*x; if(!FLT_UWORD_IS_FINITE(ix)) return x*x;
if((ix>>=23)==0) /* IEEE 754 logb */ if((ix>>=23)==0) /* IEEE 754 logb */
return -126.0; return -126.0;
else else

View File

@ -21,3 +21,4 @@
} }
#endif /* defined(_DOUBLE_IS_32BITS) */ #endif /* defined(_DOUBLE_IS_32BITS) */

View File

@ -29,15 +29,15 @@
ix = hx&0x7fffffff; /* |x| */ ix = hx&0x7fffffff; /* |x| */
iy = hy&0x7fffffff; /* |y| */ iy = hy&0x7fffffff; /* |y| */
if((ix>0x7f800000) || /* x is nan */ if(FLT_UWORD_IS_NAN(ix) ||
(iy>0x7f800000)) /* y is nan */ FLT_UWORD_IS_NAN(iy))
return x+y; return x+y;
if(x==y) return x; /* x=y, return x */ if(x==y) return x; /* x=y, return x */
if(ix==0) { /* x == 0 */ if(FLT_UWORD_IS_ZERO(ix)) { /* x == 0 */
SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */ SET_FLOAT_WORD(x,(hy&0x80000000)|FLT_UWORD_MIN);
y = x*x; y = x*x;
if(y==x) return y; else return x; /* raise underflow flag */ if(y==x) return y; else return x; /* raise underflow flag */
} }
if(hx>=0) { /* x > 0 */ if(hx>=0) { /* x > 0 */
if(hx>hy) { /* x > y, x -= ulp */ if(hx>hy) { /* x > y, x -= ulp */
hx -= 1; hx -= 1;
@ -52,7 +52,7 @@
} }
} }
hy = hx&0x7f800000; hy = hx&0x7f800000;
if(hy>=0x7f800000) return x+x; /* overflow */ if(hy>FLT_UWORD_MAX) return x+x; /* overflow */
if(hy<0x00800000) { /* underflow */ if(hy<0x00800000) { /* underflow */
y = x*x; y = x*x;
if(y!=x) { /* raise underflow flag */ if(y!=x) { /* raise underflow flag */

View File

@ -33,15 +33,17 @@ TWO23[2]={
#endif #endif
{ {
__int32_t i0,j0,sx; __int32_t i0,j0,sx;
__uint32_t i,i1; __uint32_t i,i1,ix;
float t; float t;
volatile float w; volatile float w;
GET_FLOAT_WORD(i0,x); GET_FLOAT_WORD(i0,x);
sx = (i0>>31)&1; sx = (i0>>31)&1;
j0 = ((i0>>23)&0xff)-0x7f; ix = (i0&0x7fffffff);
j0 = (ix>>23)-0x7f;
if(j0<23) { if(j0<23) {
if(j0<0) { if(FLT_UWORD_IS_ZERO(ix))
if((i0&0x7fffffff)==0) return x; return x;
if(j0<0) {
i1 = (i0&0x07fffff); i1 = (i0&0x07fffff);
i0 &= 0xfff00000; i0 &= 0xfff00000;
i0 |= ((i1|-i1)>>9)&0x400000; i0 |= ((i1|-i1)>>9)&0x400000;
@ -58,8 +60,9 @@ TWO23[2]={
if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0); if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0);
} }
} else { } else {
if(j0==0x80) return x+x; /* inf or NaN */ if(!FLT_UWORD_IS_FINITE(ix)) return x+x; /* inf or NaN */
else return x; /* x is integral */ else
return x; /* x is integral */
} }
SET_FLOAT_WORD(x,i0); SET_FLOAT_WORD(x,i0);
w = TWO23[sx]+x; w = TWO23[sx]+x;

View File

@ -15,6 +15,7 @@
#include "fdlibm.h" #include "fdlibm.h"
#include <limits.h> #include <limits.h>
#include <float.h>
#if INT_MAX > 50000 #if INT_MAX > 50000
#define OVERFLOW_INT 50000 #define OVERFLOW_INT 50000
@ -40,25 +41,30 @@ tiny = 1.0e-30;
#endif #endif
{ {
__int32_t k,ix; __int32_t k,ix;
__uint32_t hx;
GET_FLOAT_WORD(ix,x); GET_FLOAT_WORD(ix,x);
k = (ix&0x7f800000)>>23; /* extract exponent */ hx = ix&0x7fffffff;
if (k==0) { /* 0 or subnormal x */ k = hx>>23; /* extract exponent */
if ((ix&0x7fffffff)==0) return x; /* +-0 */ if (FLT_UWORD_IS_ZERO(hx))
return x;
if (!FLT_UWORD_IS_FINITE(hx))
return x+x; /* NaN or Inf */
if (FLT_UWORD_IS_SUBNORMAL(hx)) {
x *= two25; x *= two25;
GET_FLOAT_WORD(ix,x); GET_FLOAT_WORD(ix,x);
k = ((ix&0x7f800000)>>23) - 25; k = ((ix&0x7f800000)>>23) - 25;
if (n< -50000) return tiny*x; /*underflow*/ if (n< -50000) return tiny*x; /*underflow*/
} }
if (k==0xff) return x+x; /* NaN or Inf */
k = k+n; k = k+n;
if (k > 0xfe) return huge*copysignf(huge,x); /* overflow */ if (k > FLT_LARGEST_EXP) return huge*copysignf(huge,x); /* overflow */
if (k > 0) /* normal result */ if (k > 0) /* normal result */
{SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;} {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
if (k <= -25) { if (k < FLT_SMALLEST_EXP) {
if (n > OVERFLOW_INT) /* in case integer overflow in n+k */ if (n > OVERFLOW_INT) /* in case integer overflow in n+k */
return huge*copysignf(huge,x); /*overflow*/ return huge*copysignf(huge,x); /*overflow*/
else return tiny*copysignf(tiny,x); /*underflow*/ else return tiny*copysignf(tiny,x); /*underflow*/
} }
k += 25; /* subnormal result */ k += 25; /* subnormal result */
SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
return x*twom25; return x*twom25;

View File

@ -37,7 +37,7 @@ ln2 = 6.9314718246e-01; /* 0x3f317218 */
if(hx<0x3f800000) { /* x < 1 */ if(hx<0x3f800000) { /* x < 1 */
return (x-x)/(x-x); return (x-x)/(x-x);
} else if(hx >=0x4d800000) { /* x > 2**28 */ } else if(hx >=0x4d800000) { /* x > 2**28 */
if(hx >=0x7f800000) { /* x is inf of NaN */ if(!FLT_UWORD_IS_FINITE(hx)) { /* x is inf of NaN */
return x+x; return x+x;
} else } else
return __ieee754_logf(x)+ln2; /* acosh(huge)=log(2x) */ return __ieee754_logf(x)+ln2; /* acosh(huge)=log(2x) */

View File

@ -42,14 +42,14 @@ pi_lo = 1.5099578832e-07; /* 0x34222168 */
ix = hx&0x7fffffff; ix = hx&0x7fffffff;
GET_FLOAT_WORD(hy,y); GET_FLOAT_WORD(hy,y);
iy = hy&0x7fffffff; iy = hy&0x7fffffff;
if((ix>0x7f800000)|| if(FLT_UWORD_IS_NAN(ix)||
(iy>0x7f800000)) /* x or y is NaN */ FLT_UWORD_IS_NAN(iy)) /* x or y is NaN */
return x+y; return x+y;
if(hx==0x3f800000) return atanf(y); /* x=1.0 */ if(hx==0x3f800000) return atanf(y); /* x=1.0 */
m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */ m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */
/* when y = 0 */ /* when y = 0 */
if(iy==0) { if(FLT_UWORD_IS_ZERO(iy)) {
switch(m) { switch(m) {
case 0: case 0:
case 1: return y; /* atan(+-0,+anything)=+-0 */ case 1: return y; /* atan(+-0,+anything)=+-0 */
@ -58,11 +58,11 @@ pi_lo = 1.5099578832e-07; /* 0x34222168 */
} }
} }
/* when x = 0 */ /* when x = 0 */
if(ix==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; if(FLT_UWORD_IS_ZERO(ix)) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
/* when x is INF */ /* when x is INF */
if(ix==0x7f800000) { if(FLT_UWORD_IS_INFINITE(ix)) {
if(iy==0x7f800000) { if(FLT_UWORD_IS_INFINITE(iy)) {
switch(m) { switch(m) {
case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */ case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */
case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */ case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
@ -79,7 +79,7 @@ pi_lo = 1.5099578832e-07; /* 0x34222168 */
} }
} }
/* when y is INF */ /* when y is INF */
if(iy==0x7f800000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; if(FLT_UWORD_IS_INFINITE(iy)) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
/* compute y/x */ /* compute y/x */
k = (iy-ix)>>23; k = (iy-ix)>>23;

View File

@ -39,7 +39,7 @@ static float one = 1.0, half=0.5, huge = 1.0e30;
ix &= 0x7fffffff; ix &= 0x7fffffff;
/* x is INF or NaN */ /* x is INF or NaN */
if(ix>=0x7f800000) return x*x; if(!FLT_UWORD_IS_FINITE(ix)) return x*x;
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if(ix<0x3eb17218) { if(ix<0x3eb17218) {
@ -56,10 +56,11 @@ static float one = 1.0, half=0.5, huge = 1.0e30;
} }
/* |x| in [22, log(maxdouble)] return half*exp(|x|) */ /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
if (ix < 0x42b17180) return half*__ieee754_expf(fabsf(x)); if (ix <= FLT_UWORD_LOG_MAX)
return half*__ieee754_expf(fabsf(x));
/* |x| in [log(maxdouble), overflowthresold] */ /* |x| in [log(maxdouble), overflowthresold] */
if (ix<=0x42b2d4fc) { if (ix <= FLT_UWORD_LOG_2MAX) {
w = __ieee754_expf(half*fabsf(x)); w = __ieee754_expf(half*fabsf(x));
t = half*w; t = half*w;
return t*w; return t*w;

View File

@ -28,8 +28,6 @@ one = 1.0,
halF[2] = {0.5,-0.5,}, halF[2] = {0.5,-0.5,},
huge = 1.0e+30, huge = 1.0e+30,
twom100 = 7.8886090522e-31, /* 2**-100=0x0d800000 */ twom100 = 7.8886090522e-31, /* 2**-100=0x0d800000 */
o_threshold= 8.8721679688e+01, /* 0x42b17180 */
u_threshold= -1.0397208405e+02, /* 0xc2cff1b5 */
ln2HI[2] ={ 6.9313812256e-01, /* 0x3f317180 */ ln2HI[2] ={ 6.9313812256e-01, /* 0x3f317180 */
-6.9313812256e-01,}, /* 0xbf317180 */ -6.9313812256e-01,}, /* 0xbf317180 */
ln2LO[2] ={ 9.0580006145e-06, /* 0x3717f7d1 */ ln2LO[2] ={ 9.0580006145e-06, /* 0x3717f7d1 */
@ -49,23 +47,23 @@ P5 = 4.1381369442e-08; /* 0x3331bb4c */
#endif #endif
{ {
float y,hi,lo,c,t; float y,hi,lo,c,t;
__int32_t k,xsb; __int32_t k,xsb,sx;
__uint32_t hx; __uint32_t hx;
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(sx,x);
xsb = (hx>>31)&1; /* sign bit of x */ xsb = (sx>>31)&1; /* sign bit of x */
hx &= 0x7fffffff; /* high word of |x| */ hx = sx & 0x7fffffff; /* high word of |x| */
/* filter out non-finite argument */ /* filter out non-finite argument */
if(hx >= 0x42b17218) { /* if |x|>=88.721... */ if(FLT_UWORD_IS_NAN(hx))
if(hx>0x7f800000) return x+x; /* NaN */
return x+x; /* NaN */ if(FLT_UWORD_IS_INFINITE(hx))
if(hx==0x7f800000) return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ if(sx > FLT_UWORD_LOG_MAX)
if(x > o_threshold) return huge*huge; /* overflow */ return huge*huge; /* overflow */
if(x < u_threshold) return twom100*twom100; /* underflow */ if(sx < 0 && hx > FLT_UWORD_LOG_MIN)
} return twom100*twom100; /* underflow */
/* argument reduction */ /* argument reduction */
if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */

View File

@ -43,20 +43,23 @@ static float one = 1.0, Zero[] = {0.0, -0.0,};
hy &= 0x7fffffff; /* |y| */ hy &= 0x7fffffff; /* |y| */
/* purge off exception values */ /* purge off exception values */
if(hy==0||(hx>=0x7f800000)|| /* y=0,or x not finite */ if(FLT_UWORD_IS_ZERO(hy)||
(hy>0x7f800000)) /* or y is NaN */ !FLT_UWORD_IS_FINITE(hx)||
FLT_UWORD_IS_NAN(hy))
return (x*y)/(x*y); return (x*y)/(x*y);
if(hx<hy) return x; /* |x|<|y| return x */ if(hx<hy) return x; /* |x|<|y| return x */
if(hx==hy) if(hx==hy)
return Zero[(__uint32_t)sx>>31]; /* |x|=|y| return x*0*/ return Zero[(__uint32_t)sx>>31]; /* |x|=|y| return x*0*/
/* Note: y cannot be zero if we reach here. */
/* determine ix = ilogb(x) */ /* determine ix = ilogb(x) */
if(hx<0x00800000) { /* subnormal x */ if(FLT_UWORD_IS_SUBNORMAL(hx)) { /* subnormal x */
for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1; for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
} else ix = (hx>>23)-127; } else ix = (hx>>23)-127;
/* determine iy = ilogb(y) */ /* determine iy = ilogb(y) */
if(hy<0x00800000) { /* subnormal y */ if(FLT_UWORD_IS_SUBNORMAL(hy)) { /* subnormal y */
for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1; for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1;
} else iy = (hy>>23)-127; } else iy = (hy>>23)-127;
@ -99,6 +102,8 @@ static float one = 1.0, Zero[] = {0.0, -0.0,};
hx = ((hx-0x00800000)|((iy+127)<<23)); hx = ((hx-0x00800000)|((iy+127)<<23));
SET_FLOAT_WORD(x,hx|sx); SET_FLOAT_WORD(x,hx|sx);
} else { /* subnormal output */ } else { /* subnormal output */
/* If denormals are not supported, this code will generate a
zero representation. */
n = -126 - iy; n = -126 - iy;
hx >>= n; hx >>= n;
SET_FLOAT_WORD(x,hx|sx); SET_FLOAT_WORD(x,hx|sx);

View File

@ -35,10 +35,10 @@
if((ha-hb)>0xf000000L) {return a+b;} /* x/y > 2**30 */ if((ha-hb)>0xf000000L) {return a+b;} /* x/y > 2**30 */
k=0; k=0;
if(ha > 0x58800000L) { /* a>2**50 */ if(ha > 0x58800000L) { /* a>2**50 */
if(ha >= 0x7f800000L) { /* Inf or NaN */ if(!FLT_UWORD_IS_FINITE(ha)) { /* Inf or NaN */
w = a+b; /* for sNaN */ w = a+b; /* for sNaN */
if(ha == 0x7f800000L) w = a; if(FLT_UWORD_IS_INFINITE(ha)) w = a;
if(hb == 0x7f800000L) w = b; if(FLT_UWORD_IS_INFINITE(hb)) w = b;
return w; return w;
} }
/* scale a and b by 2**-60 */ /* scale a and b by 2**-60 */
@ -47,8 +47,9 @@
SET_FLOAT_WORD(b,hb); SET_FLOAT_WORD(b,hb);
} }
if(hb < 0x26800000L) { /* b < 2**-50 */ if(hb < 0x26800000L) { /* b < 2**-50 */
if(hb <= 0x007fffffL) { /* subnormal b or 0 */ if(FLT_UWORD_IS_ZERO(hb)) {
if(hb==0) return a; return a;
} else if(FLT_UWORD_IS_SUBNORMAL(hb)) {
SET_FLOAT_WORD(t1,0x3f000000L); /* t1=2^126 */ SET_FLOAT_WORD(t1,0x3f000000L); /* t1=2^126 */
b *= t1; b *= t1;
a *= t1; a *= t1;

View File

@ -58,14 +58,14 @@ static float zero = 0.0;
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff; ix = hx&0x7fffffff;
if(ix>=0x7f800000) return one/(x*x); if(!FLT_UWORD_IS_FINITE(ix)) return one/(x*x);
x = fabsf(x); x = fabsf(x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */ if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x); s = sinf(x);
c = cosf(x); c = cosf(x);
ss = s-c; ss = s-c;
cc = s+c; cc = s+c;
if(ix<0x7f000000) { /* make sure x+x not overflow */ if(ix<=FLT_UWORD_HALF_MAX) { /* make sure x+x not overflow */
z = -cosf(x+x); z = -cosf(x+x);
if ((s*c)<zero) cc = z/ss; if ((s*c)<zero) cc = z/ss;
else ss = z/cc; else ss = z/cc;
@ -128,8 +128,8 @@ v04 = 4.4111031494e-10; /* 0x2ff280c2 */
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx; ix = 0x7fffffff&hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if(ix>=0x7f800000) return one/(x+x*x); if(!FLT_UWORD_IS_FINITE(ix)) return one/(x+x*x);
if(ix==0) return -one/zero; if(FLT_UWORD_IS_ZERO(ix)) return -one/zero;
if(hx<0) return zero/zero; if(hx<0) return zero/zero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */ if(ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
@ -151,7 +151,7 @@ v04 = 4.4111031494e-10; /* 0x2ff280c2 */
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/ */
if(ix<0x7f000000) { /* make sure x+x not overflow */ if(ix<=FLT_UWORD_HALF_MAX) { /* make sure x+x not overflow */
z = -cosf(x+x); z = -cosf(x+x);
if ((s*c)<zero) cc = z/ss; if ((s*c)<zero) cc = z/ss;
else ss = z/cc; else ss = z/cc;

View File

@ -59,14 +59,14 @@ static float zero = 0.0;
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff; ix = hx&0x7fffffff;
if(ix>=0x7f800000) return one/x; if(!FLT_UWORD_IS_FINITE(ix)) return one/x;
y = fabsf(x); y = fabsf(x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */ if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(y); s = sinf(y);
c = cosf(y); c = cosf(y);
ss = -s-c; ss = -s-c;
cc = s-c; cc = s-c;
if(ix<0x7f000000) { /* make sure y+y not overflow */ if(ix<=FLT_UWORD_HALF_MAX) { /* make sure y+y not overflow */
z = cosf(y+y); z = cosf(y+y);
if ((s*c)>zero) cc = z/ss; if ((s*c)>zero) cc = z/ss;
else ss = z/cc; else ss = z/cc;
@ -129,15 +129,15 @@ static float V0[5] = {
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx; ix = 0x7fffffff&hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if(ix>=0x7f800000) return one/(x+x*x); if(!FLT_UWORD_IS_FINITE(ix)) return one/(x+x*x);
if(ix==0) return -one/zero; if(FLT_UWORD_IS_ZERO(ix)) return -one/zero;
if(hx<0) return zero/zero; if(hx<0) return zero/zero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */ if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x); s = sinf(x);
c = cosf(x); c = cosf(x);
ss = -s-c; ss = -s-c;
cc = s-c; cc = s-c;
if(ix<0x7f000000) { /* make sure x+x not overflow */ if(ix<=FLT_UWORD_HALF_MAX) { /* make sure x+x not overflow */
z = cosf(x+x); z = cosf(x+x);
if ((s*c)>zero) cc = z/ss; if ((s*c)>zero) cc = z/ss;
else ss = z/cc; else ss = z/cc;

View File

@ -47,7 +47,7 @@ static float zero = 0.0000000000e+00;
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx; ix = 0x7fffffff&hx;
/* if J(n,NaN) is NaN */ /* if J(n,NaN) is NaN */
if(ix>0x7f800000) return x+x; if(FLT_UWORD_IS_NAN(ix)) return x+x;
if(n<0){ if(n<0){
n = -n; n = -n;
x = -x; x = -x;
@ -57,7 +57,7 @@ static float zero = 0.0000000000e+00;
if(n==1) return(__ieee754_j1f(x)); if(n==1) return(__ieee754_j1f(x));
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabsf(x); x = fabsf(x);
if(ix==0||ix>=0x7f800000) /* if x is 0 or inf */ if(FLT_UWORD_IS_ZERO(ix)||FLT_UWORD_IS_INFINITE(ix))
b = zero; b = zero;
else if((float)n<=x) { else if((float)n<=x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
@ -181,8 +181,8 @@ static float zero = 0.0000000000e+00;
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx; ix = 0x7fffffff&hx;
/* if Y(n,NaN) is NaN */ /* if Y(n,NaN) is NaN */
if(ix>0x7f800000) return x+x; if(FLT_UWORD_IS_NAN(ix)) return x+x;
if(ix==0) return -one/zero; if(FLT_UWORD_IS_ZERO(ix)) return -one/zero;
if(hx<0) return zero/zero; if(hx<0) return zero/zero;
sign = 1; sign = 1;
if(n<0){ if(n<0){
@ -191,7 +191,7 @@ static float zero = 0.0000000000e+00;
} }
if(n==0) return(__ieee754_y0f(x)); if(n==0) return(__ieee754_y0f(x));
if(n==1) return(sign*__ieee754_y1f(x)); if(n==1) return(sign*__ieee754_y1f(x));
if(ix==0x7f800000) return zero; if(FLT_UWORD_IS_INFINITE(ix)) return zero;
a = __ieee754_y0f(x); a = __ieee754_y0f(x);
b = __ieee754_y1f(x); b = __ieee754_y1f(x);

View File

@ -50,14 +50,14 @@ static float zero = 0.0;
GET_FLOAT_WORD(ix,x); GET_FLOAT_WORD(ix,x);
k=0; k=0;
if (ix < 0x00800000) { /* x < 2**-126 */ if (FLT_UWORD_IS_ZERO(ix&0x7fffffff))
if ((ix&0x7fffffff)==0) return -two25/zero; /* log(+-0)=-inf */
return -two25/zero; /* log(+-0)=-inf */ if (ix<0) return (x-x)/zero; /* log(-#) = NaN */
if (ix<0) return (x-x)/zero; /* log(-#) = NaN */ if (!FLT_UWORD_IS_FINITE(ix)) return x+x;
if (FLT_UWORD_IS_SUBNORMAL(ix)) {
k -= 25; x *= two25; /* subnormal number, scale up x */ k -= 25; x *= two25; /* subnormal number, scale up x */
GET_FLOAT_WORD(ix,x); GET_FLOAT_WORD(ix,x);
} }
if (ix >= 0x7f800000) return x+x;
k += (ix>>23)-127; k += (ix>>23)-127;
ix &= 0x007fffff; ix &= 0x007fffff;
i = (ix+(0x95f64<<3))&0x800000; i = (ix+(0x95f64<<3))&0x800000;

View File

@ -44,14 +44,14 @@ static float zero = 0.0;
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
k=0; k=0;
if (hx < 0x00800000) { /* x < 2**-126 */ if (FLT_UWORD_IS_ZERO(hx&0x7fffffff))
if ((hx&0x7fffffff)==0) return -two25/zero; /* log(+-0)=-inf */
return -two25/zero; /* log(+-0)=-inf */ if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ if (!FLT_UWORD_IS_FINITE(hx)) return x+x;
if (FLT_UWORD_IS_SUBNORMAL(hx)) {
k -= 25; x *= two25; /* subnormal number, scale up x */ k -= 25; x *= two25; /* subnormal number, scale up x */
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
} }
if (hx >= 0x7f800000) return x+x;
k += (hx>>23)-127; k += (hx>>23)-127;
i = ((__uint32_t)k&0x80000000)>>31; i = ((__uint32_t)k&0x80000000)>>31;
hx = (hx&0x007fffff)|((0x7f-i)<<23); hx = (hx&0x007fffff)|((0x7f-i)<<23);

View File

@ -73,11 +73,11 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
ix = hx&0x7fffffff; iy = hy&0x7fffffff; ix = hx&0x7fffffff; iy = hy&0x7fffffff;
/* y==zero: x**0 = 1 */ /* y==zero: x**0 = 1 */
if(iy==0) return one; if(FLT_UWORD_IS_ZERO(iy)) return one;
/* +-NaN return x+y */ /* +-NaN return x+y */
if(ix > 0x7f800000 || if(FLT_UWORD_IS_NAN(ix) ||
iy > 0x7f800000) FLT_UWORD_IS_NAN(iy))
return x+y; return x+y;
/* determine if y is an odd int when x < 0 /* determine if y is an odd int when x < 0
@ -96,7 +96,7 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
} }
/* special value of y */ /* special value of y */
if (iy==0x7f800000) { /* y is +-inf */ if (FLT_UWORD_IS_INFINITE(iy)) { /* y is +-inf */
if (ix==0x3f800000) if (ix==0x3f800000)
return y - y; /* inf**+-1 is NaN */ return y - y; /* inf**+-1 is NaN */
else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */ else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
@ -115,7 +115,7 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
ax = fabsf(x); ax = fabsf(x);
/* special value of x */ /* special value of x */
if(ix==0x7f800000||ix==0||ix==0x3f800000){ if(FLT_UWORD_IS_INFINITE(ix)||FLT_UWORD_IS_ZERO(ix)||ix==0x3f800000){
z = ax; /*x is +-0,+-inf,+-1*/ z = ax; /*x is +-0,+-inf,+-1*/
if(hy<0) z = one/z; /* z = (1/|x|) */ if(hy<0) z = one/z; /* z = (1/|x|) */
if(hx<0) { if(hx<0) {
@ -149,7 +149,7 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
float s2,s_h,s_l,t_h,t_l; float s2,s_h,s_l,t_h,t_l;
n = 0; n = 0;
/* take care subnormal number */ /* take care subnormal number */
if(ix<0x00800000) if(FLT_UWORD_IS_SUBNORMAL(ix))
{ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); } {ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); }
n += ((ix)>>23)-0x7f; n += ((ix)>>23)-0x7f;
j = ix&0x007fffff; j = ix&0x007fffff;
@ -209,20 +209,21 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
p_h = y1*t1; p_h = y1*t1;
z = p_l+p_h; z = p_l+p_h;
GET_FLOAT_WORD(j,z); GET_FLOAT_WORD(j,z);
if (j>0x43000000) /* if z > 128 */ i = j&0x7fffffff;
return s*huge*huge; /* overflow */ if (j>0) {
else if (j==0x43000000) { /* if z == 128 */ if (i>FLT_UWORD_EXP_MAX)
if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ return s*huge*huge; /* overflow */
} else if (i==FLT_UWORD_EXP_MAX)
else if ((j&0x7fffffff)>0x43160000) /* z <= -150 */ if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
return s*tiny*tiny; /* underflow */ } else {
else if (j==0xc3160000){ /* z == -150 */ if (i>FLT_UWORD_EXP_MIN)
if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ return s*tiny*tiny; /* underflow */
else if (i==FLT_UWORD_EXP_MIN)
if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
} }
/* /*
* compute 2**(p_h+p_l) * compute 2**(p_h+p_l)
*/ */
i = j&0x7fffffff;
k = (i>>23)-0x7f; k = (i>>23)-0x7f;
n = 0; n = 0;
if(i>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */ if(i>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */

View File

@ -174,7 +174,7 @@ pio2_3t = 6.1232342629e-17; /* 0x248d3132 */
/* /*
* all other (large) arguments * all other (large) arguments
*/ */
if(ix>=0x7f800000) { /* x is inf or NaN */ if(!FLT_UWORD_IS_FINITE(ix)) {
y[0]=y[1]=x-x; return 0; y[0]=y[1]=x-x; return 0;
} }
/* set z = scalbn(|x|,ilogb(x)-7) */ /* set z = scalbn(|x|,ilogb(x)-7) */

View File

@ -40,13 +40,13 @@ static float zero = 0.0;
hx &= 0x7fffffff; hx &= 0x7fffffff;
/* purge off exception values */ /* purge off exception values */
if(hp==0) return (x*p)/(x*p); /* p = 0 */ if(FLT_UWORD_IS_ZERO(hp)||
if((hx>=0x7f800000)|| /* x not finite */ !FLT_UWORD_IS_FINITE(hx)||
((hp>0x7f800000))) /* p is NaN */ FLT_UWORD_IS_NAN(hp))
return (x*p)/(x*p); return (x*p)/(x*p);
if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */ if (hp<=FLT_UWORD_HALF_MAX) x = __ieee754_fmodf(x,p+p); /* now x < 2p */
if ((hx-hp)==0) return zero*x; if ((hx-hp)==0) return zero*x;
x = fabsf(x); x = fabsf(x);
p = fabsf(p); p = fabsf(p);

View File

@ -35,7 +35,7 @@ static float one = 1.0, shuge = 1.0e37;
ix = jx&0x7fffffff; ix = jx&0x7fffffff;
/* x is INF or NaN */ /* x is INF or NaN */
if(ix>=0x7f800000) return x+x; if(!FLT_UWORD_IS_FINITE(ix)) return x+x;
h = 0.5; h = 0.5;
if (jx<0) h = -h; if (jx<0) h = -h;
@ -49,10 +49,10 @@ static float one = 1.0, shuge = 1.0e37;
} }
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */ /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x42b17180) return h*__ieee754_expf(fabsf(x)); if (ix<=FLT_UWORD_LOG_MAX) return h*__ieee754_expf(fabsf(x));
/* |x| in [log(maxdouble), overflowthresold] */ /* |x| in [log(maxdouble), overflowthresold] */
if (ix<=0x42b2d4fc) { if (ix<=FLT_UWORD_LOG_2MAX) {
w = __ieee754_expf((float)0.5*fabsf(x)); w = __ieee754_expf((float)0.5*fabsf(x));
t = h*w; t = h*w;
return t*w; return t*w;

View File

@ -30,25 +30,23 @@ static float one = 1.0, tiny=1.0e-30;
{ {
float z; float z;
__int32_t sign = (__int32_t)0x80000000; __int32_t sign = (__int32_t)0x80000000;
__uint32_t r; __uint32_t r,hx;
__int32_t ix,s,q,m,t,i; __int32_t ix,s,q,m,t,i;
GET_FLOAT_WORD(ix,x); GET_FLOAT_WORD(ix,x);
hx = ix&0x7fffffff;
/* take care of Inf and NaN */ /* take care of Inf and NaN */
if((ix&0x7f800000L)==0x7f800000L) { if(!FLT_UWORD_IS_FINITE(hx))
return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
sqrt(-inf)=sNaN */ sqrt(-inf)=sNaN */
} /* take care of zero and -ves */
/* take care of zero */ if(FLT_UWORD_IS_ZERO(hx)) return x;/* sqrt(+-0) = +-0 */
if(ix<=0) { if(ix<0) return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */
else if(ix<0)
return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
}
/* normalize x */ /* normalize x */
m = (ix>>23); m = (ix>>23);
if(m==0) { /* subnormal x */ if(FLT_UWORD_IS_SUBNORMAL(hx)) { /* subnormal x */
for(i=0;(ix&0x00800000L)==0;i++) ix<<=1; for(i=0;(ix&0x00800000L)==0;i++) ix<<=1;
m -= i-1; m -= i-1;
} }

View File

@ -35,7 +35,7 @@ huge= 1.0000000000e+30;
__int32_t hx,ix; __int32_t hx,ix;
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff; ix = hx&0x7fffffff;
if(ix>=0x7f800000) return x+x; /* x is inf or NaN */ if(!FLT_UWORD_IS_FINITE(ix)) return x+x; /* x is inf or NaN */
if(ix< 0x31800000) { /* |x|<2**-28 */ if(ix< 0x31800000) { /* |x|<2**-28 */
if(huge+x>one) return x; /* return x inexact except 0 */ if(huge+x>one) return x; /* return x inexact except 0 */
} }

View File

@ -77,7 +77,7 @@ huge = 1.0e30;
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff; ix = hx&0x7fffffff;
if(ix>=0x50800000) { /* if |x| >= 2^34 */ if(ix>=0x50800000) { /* if |x| >= 2^34 */
if(ix>0x7f800000) if(FLT_UWORD_IS_NAN(ix))
return x+x; /* NaN */ return x+x; /* NaN */
if(hx>0) return atanhi[3]+atanlo[3]; if(hx>0) return atanhi[3]+atanlo[3];
else return -atanhi[3]-atanlo[3]; else return -atanhi[3]-atanlo[3];

View File

@ -29,14 +29,15 @@ static float huge = 1.0e30;
#endif #endif
{ {
__int32_t i0,j0; __int32_t i0,j0;
__uint32_t i; __uint32_t i,ix;
GET_FLOAT_WORD(i0,x); GET_FLOAT_WORD(i0,x);
j0 = ((i0>>23)&0xff)-0x7f; ix = (i0&0x7fffffff);
j0 = (ix>>23)-0x7f;
if(j0<23) { if(j0<23) {
if(j0<0) { /* raise inexact if x != 0 */ if(j0<0) { /* raise inexact if x != 0 */
if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */ if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */
if(i0<0) {i0=0x80000000;} if(i0<0) {i0=0x80000000;}
else if(i0!=0) { i0=0x3f800000;} else if(!FLT_UWORD_IS_ZERO(ix)) { i0=0x3f800000;}
} }
} else { } else {
i = (0x007fffff)>>j0; i = (0x007fffff)>>j0;
@ -47,7 +48,7 @@ static float huge = 1.0e30;
} }
} }
} else { } else {
if(j0==0x80) return x+x; /* inf or NaN */ if(!FLT_UWORD_IS_FINITE(ix)) return x+x; /* inf or NaN */
else return x; /* x is integral */ else return x; /* x is integral */
} }
SET_FLOAT_WORD(x,i0); SET_FLOAT_WORD(x,i0);

View File

@ -38,7 +38,7 @@ static float one=1.0;
if(ix <= 0x3f490fd8) return __kernel_cosf(x,z); if(ix <= 0x3f490fd8) return __kernel_cosf(x,z);
/* cos(Inf or NaN) is NaN */ /* cos(Inf or NaN) is NaN */
else if (ix>=0x7f800000) return x-x; else if (!FLT_UWORD_IS_FINITE(ix)) return x-x;
/* argument reduction needed */ /* argument reduction needed */
else { else {

View File

@ -109,7 +109,7 @@ sb7 = -2.2440952301e+01; /* 0xc1b38712 */
float R,S,P,Q,s,y,z,r; float R,S,P,Q,s,y,z,r;
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff; ix = hx&0x7fffffff;
if(ix>=0x7f800000) { /* erf(nan)=nan */ if(!FLT_UWORD_IS_FINITE(ix)) { /* erf(nan)=nan */
i = ((__uint32_t)hx>>31)<<1; i = ((__uint32_t)hx>>31)<<1;
return (float)(1-i)+one/x; /* erf(+-inf)=+-1 */ return (float)(1-i)+one/x; /* erf(+-inf)=+-1 */
} }
@ -166,7 +166,7 @@ sb7 = -2.2440952301e+01; /* 0xc1b38712 */
float R,S,P,Q,s,y,z,r; float R,S,P,Q,s,y,z,r;
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff; ix = hx&0x7fffffff;
if(ix>=0x7f800000) { /* erfc(nan)=nan */ if(!FLT_UWORD_IS_FINITE(ix)) { /* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */ /* erfc(+-inf)=0,2 */
return (float)(((__uint32_t)hx>>31)<<1)+one/x; return (float)(((__uint32_t)hx>>31)<<1)+one/x;
} }

View File

@ -38,14 +38,15 @@ static float huge = 1.0e30;
#endif #endif
{ {
__int32_t i0,j0; __int32_t i0,j0;
__uint32_t i; __uint32_t i,ix;
GET_FLOAT_WORD(i0,x); GET_FLOAT_WORD(i0,x);
j0 = ((i0>>23)&0xff)-0x7f; ix = (i0&0x7fffffff);
j0 = (ix>>23)-0x7f;
if(j0<23) { if(j0<23) {
if(j0<0) { /* raise inexact if x != 0 */ if(j0<0) { /* raise inexact if x != 0 */
if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */ if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */
if(i0>=0) {i0=0;} if(i0>=0) {i0=0;}
else if((i0&0x7fffffff)!=0) else if(!FLT_UWORD_IS_ZERO(ix))
{ i0=0xbf800000;} { i0=0xbf800000;}
} }
} else { } else {
@ -57,7 +58,7 @@ static float huge = 1.0e30;
} }
} }
} else { } else {
if(j0==0x80) return x+x; /* inf or NaN */ if(!FLT_UWORD_IS_FINITE(ix)) return x+x; /* inf or NaN */
else return x; /* x is integral */ else return x; /* x is integral */
} }
SET_FLOAT_WORD(x,i0); SET_FLOAT_WORD(x,i0);

View File

@ -33,8 +33,8 @@ two25 = 3.3554432000e+07; /* 0x4c000000 */
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx; ix = 0x7fffffff&hx;
*eptr = 0; *eptr = 0;
if(ix>=0x7f800000||(ix==0)) return x; /* 0,inf,nan */ if(!FLT_UWORD_IS_FINITE(ix)||FLT_UWORD_IS_ZERO(ix)) return x; /* 0,inf,nan */
if (ix<0x00800000) { /* subnormal */ if (FLT_UWORD_IS_SUBNORMAL(ix)) { /* subnormal */
x *= two25; x *= two25;
GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff; ix = hx&0x7fffffff;

View File

@ -1,6 +1,5 @@
/* /*
* isinff(x) returns 1 if x is infinity, else 0; * isinff(x) returns 1 if x is +-infinity, else 0;
* no branching!
* Added by Cygnus Support. * Added by Cygnus Support.
*/ */
@ -16,8 +15,7 @@
__int32_t ix; __int32_t ix;
GET_FLOAT_WORD(ix,x); GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff; ix &= 0x7fffffff;
ix = 0x7f800000 - ix; return FLT_UWORD_IS_INFINITE(ix);
return 1 - (int)((__uint32_t)(ix|(-ix))>>31);
} }
#ifdef _DOUBLE_IS_32BITS #ifdef _DOUBLE_IS_32BITS

View File

@ -15,7 +15,6 @@
/* /*
* isnanf(x) returns 1 is x is nan, else 0; * isnanf(x) returns 1 is x is nan, else 0;
* no branching!
*/ */
#include "fdlibm.h" #include "fdlibm.h"
@ -30,8 +29,7 @@
__int32_t ix; __int32_t ix;
GET_FLOAT_WORD(ix,x); GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff; ix &= 0x7fffffff;
ix = 0x7f800000 - ix; return FLT_UWORD_IS_NAN(ix);
return (int)(((__uint32_t)(ix))>>31);
} }
#ifdef _DOUBLE_IS_32BITS #ifdef _DOUBLE_IS_32BITS

View File

@ -32,7 +32,7 @@
if(ix <= 0x3f490fd8) return __kernel_sinf(x,z,0); if(ix <= 0x3f490fd8) return __kernel_sinf(x,z,0);
/* sin(Inf or NaN) is NaN */ /* sin(Inf or NaN) is NaN */
else if (ix>=0x7f800000) return x-x; else if (!FLT_UWORD_IS_FINITE(ix)) return x-x;
/* argument reduction needed */ /* argument reduction needed */
else { else {

View File

@ -32,7 +32,7 @@
if(ix <= 0x3f490fda) return __kernel_tanf(x,z,1); if(ix <= 0x3f490fda) return __kernel_tanf(x,z,1);
/* tan(Inf or NaN) is NaN */ /* tan(Inf or NaN) is NaN */
else if (ix>=0x7f800000) return x-x; /* NaN */ else if (!FLT_UWORD_IS_FINITE(ix)) return x-x; /* NaN */
/* argument reduction needed */ /* argument reduction needed */
else { else {

View File

@ -35,7 +35,7 @@ static float one=1.0, two=2.0, tiny = 1.0e-30;
ix = jx&0x7fffffff; ix = jx&0x7fffffff;
/* x is INF or NaN */ /* x is INF or NaN */
if(ix>=0x7f800000) { if(!FLT_UWORD_IS_FINITE(ix)) {
if (jx>=0) return one/x+one; /* tanh(+-inf)=+-1 */ if (jx>=0) return one/x+one; /* tanh(+-inf)=+-1 */
else return one/x-one; /* tanh(NaN) = NaN */ else return one/x-one; /* tanh(NaN) = NaN */
} }