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

Update to use unions.

parent 92c7c5d0
No related branches found
No related tags found
No related merge requests found
...@@ -108,15 +108,17 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ ...@@ -108,15 +108,17 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
double y,hi,lo,c,t; double y,hi,lo,c,t;
int k,xsb; int k,xsb;
unsigned hx; unsigned hx;
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 */
xsb = (hx>>31)&1; /* sign bit of x */ xsb = (hx>>31)&1; /* sign bit of x */
hx &= 0x7fffffff; /* high word of |x| */ hx &= 0x7fffffff; /* high word of |x| */
/* filter out non-finite argument */ /* filter out non-finite argument */
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:0.0; /* exp(+-inf)={inf,0} */ else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
} }
...@@ -147,10 +149,14 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ ...@@ -147,10 +149,14 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
if(k==0) return one-((x*c)/(c-2.0)-x); if(k==0) return one-((x*c)/(c-2.0)-x);
else y = one-((lo-(x*c)/(2.0-c))-hi); else y = one-((lo-(x*c)/(2.0-c))-hi);
if(k >= -1021) { if(k >= -1021) {
__HI(y) += (k<<20); /* add k to y's exponent */ ux.d = y;
ux.i[HIWORD] += (k<<20); /* add k to y's exponent */
y = ux.d;
return y; return y;
} else { } else {
__HI(y) += ((k+1000)<<20);/* add k to y's exponent */ ux.d = y;
ux.i[HIWORD] += ((k+1000)<<20);/* add k to y's exponent */
y = ux.d;
return y*twom1000; return y*twom1000;
} }
} }
...@@ -92,9 +92,11 @@ static double zero = 0.0; ...@@ -92,9 +92,11 @@ static double zero = 0.0;
double hfsq,f,s,z,R,w,t1,t2,dk; double hfsq,f,s,z,R,w,t1,t2,dk;
int k,hx,i,j; int k,hx,i,j;
unsigned lx; unsigned lx;
union { int i[2]; double d; } ux;
hx = __HI(x); /* high word of x */ ux.d = x;
lx = __LO(x); /* low word of x */ hx = ux.i[HIWORD]; /* high word of x */
lx = ux.i[LOWORD]; /* low word of x */
k=0; k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */ if (hx < 0x00100000) { /* x < 2**-1022 */
...@@ -102,13 +104,16 @@ static double zero = 0.0; ...@@ -102,13 +104,16 @@ static double zero = 0.0;
return -two54/zero; /* log(+-0)=-inf */ return -two54/zero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */ k -= 54; x *= two54; /* subnormal number, scale up x */
hx = __HI(x); /* high word of x */ ux.d = x;
hx = ux.i[HIWORD]; /* high word of x */
} }
if (hx >= 0x7ff00000) return x+x; if (hx >= 0x7ff00000) return x+x;
k += (hx>>20)-1023; k += (hx>>20)-1023;
hx &= 0x000fffff; hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000; i = (hx+0x95f64)&0x100000;
__HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */ ux.d = x;
ux.i[HIWORD] = hx|(i^0x3ff00000); /* normalize x or x/2 */
x = ux.d;
k += (i>>20); k += (i>>20);
f = x-1.0; f = x-1.0;
if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */ if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment