Skip to content
Snippets Groups Projects
Commit e19642f9 authored by Raymond Toy's avatar Raymond Toy
Browse files

Modify to use unions and change name to include fdlibm prefix.

parent ae8eb31e
No related branches found
No related tags found
No related merge requests found
...@@ -127,17 +127,20 @@ Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */ ...@@ -127,17 +127,20 @@ Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
#ifdef __STDC__ #ifdef __STDC__
double expm1(double x) double fdlibm_expm1(double x)
#else #else
double expm1(x) double fdlibm_expm1(x)
double x; double x;
#endif #endif
{ {
double y,hi,lo,c,t,e,hxs,hfx,r1; double y,hi,lo,c,t,e,hxs,hfx,r1;
int k,xsb; int k,xsb;
unsigned hx; unsigned hx;
union { int i[2]; double d; } ux;
union { int i[2]; double d; } utmp;
hx = __HI(x); /* high word of x */ ux.d = x;
hx = ux.i[HIWORD]; /* high word of x */
xsb = hx&0x80000000; /* sign bit of x */ xsb = hx&0x80000000; /* sign bit of x */
if(xsb==0) y=x; else y= -x; /* y = |x| */ if(xsb==0) y=x; else y= -x; /* y = |x| */
hx &= 0x7fffffff; /* high word of |x| */ hx &= 0x7fffffff; /* high word of |x| */
...@@ -146,7 +149,7 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */ ...@@ -146,7 +149,7 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */ if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */
if(hx >= 0x40862E42) { /* if |x|>=709.78... */ if(hx >= 0x40862E42) { /* if |x|>=709.78... */
if(hx>=0x7ff00000) { if(hx>=0x7ff00000) {
if(((hx&0xfffff)|__LO(x))!=0) if(((hx&0xfffff)|ux.i[LOWORD])!=0)
return x+x; /* NaN */ return x+x; /* NaN */
else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
} }
...@@ -196,19 +199,29 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */ ...@@ -196,19 +199,29 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
else return one+2.0*(x-e); else return one+2.0*(x-e);
if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
y = one-(e-x); y = one-(e-x);
__HI(y) += (k<<20); /* add k to y's exponent */ utmp.d = y;
utmp.i[HIWORD] += (k<<20); /* add k to y's exponent */
y = utmp.d;
return y-one; return y-one;
} }
t = one; t = one;
if(k<20) { if(k<20) {
__HI(t) = 0x3ff00000 - (0x200000>>k); /* t=1-2^-k */ utmp.d = t;
utmp.i[HIWORD] = 0x3ff00000 - (0x200000>>k); /* t=1-2^-k */
t = utmp.d;
y = t-(e-x); y = t-(e-x);
__HI(y) += (k<<20); /* add k to y's exponent */ utmp.d = y;
utmp.i[HIWORD] += (k<<20); /* add k to y's exponent */
y = utmp.d;
} else { } else {
__HI(t) = ((0x3ff-k)<<20); /* 2^-k */ utmp.d = t;
utmp.i[HIWORD] = ((0x3ff-k)<<20); /* 2^-k */
t = utmp.d;
y = x-(e+t); y = x-(e+t);
y += one; y += one;
__HI(y) += (k<<20); /* add k to y's exponent */ utmp.d = y;
utmp.i[HIWORD] += (k<<20); /* add k to y's exponent */
y = utmp.d;
} }
} }
return y; return y;
......
...@@ -97,16 +97,18 @@ Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ ...@@ -97,16 +97,18 @@ Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static double zero = 0.0; static double zero = 0.0;
#ifdef __STDC__ #ifdef __STDC__
double log1p(double x) double fdlibm_log1p(double x)
#else #else
double log1p(x) double fdlibm_log1p(x)
double x; double x;
#endif #endif
{ {
double hfsq,f,c,s,z,R,u; double hfsq,f,c,s,z,R,u;
int k,hx,hu,ax; int k,hx,hu,ax;
union { int i[2]; double d; } ux;
hx = __HI(x); /* high word of x */ ux.d = x;
hx = ux.i[HIWORD]; /* high word of x */
ax = hx&0x7fffffff; ax = hx&0x7fffffff;
k = 1; k = 1;
...@@ -128,23 +130,29 @@ static double zero = 0.0; ...@@ -128,23 +130,29 @@ static double zero = 0.0;
if (hx >= 0x7ff00000) return x+x; if (hx >= 0x7ff00000) return x+x;
if(k!=0) { if(k!=0) {
if(hx<0x43400000) { if(hx<0x43400000) {
u = 1.0+x; u = 1.0+x;
hu = __HI(u); /* high word of u */ ux.d = u;
hu = ux.i[HIWORD]; /* high word of u */
k = (hu>>20)-1023; k = (hu>>20)-1023;
c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */ c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */
c /= u; c /= u;
} else { } else {
u = x; u = x;
hu = __HI(u); /* high word of u */ ux.d = u;
hu = ux.i[HIWORD]; /* high word of u */
k = (hu>>20)-1023; k = (hu>>20)-1023;
c = 0; c = 0;
} }
hu &= 0x000fffff; hu &= 0x000fffff;
if(hu<0x6a09e) { if(hu<0x6a09e) {
__HI(u) = hu|0x3ff00000; /* normalize u */ ux.d = u;
ux.i[HIWORD] = hu|0x3ff00000; /* normalize u */
u = ux.d;
} else { } else {
k += 1; k += 1;
__HI(u) = hu|0x3fe00000; /* normalize u/2 */ ux.d = u;
ux.i[HIWORD] = hu|0x3fe00000; /* normalize u/2 */
u = ux.d;
hu = (0x00100000-hu)>>2; hu = (0x00100000-hu)>>2;
} }
f = u-1.0; f = u-1.0;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment