reverting old volatile fix for 24892 and replacing with cls' union patch.

spider and xpshell now pass the related testcases.
author=cls, r=me
This commit is contained in:
rginda%netscape.com 2000-09-11 20:56:33 +00:00
Родитель d50223fd65
Коммит 84d6c25c59
49 изменённых файлов: 489 добавлений и 247 удалений

Просмотреть файл

@ -89,7 +89,7 @@ include $(topsrcdir)/config/rules.mk
#
# Default IEEE libm
#
CFLAGS += -D_IEEE_LIBM -DGCC_OPT_BUG
CFLAGS += -D_IEEE_LIBM
ifeq ($(OS_ARCH),Linux)
LDFLAGS += -ldl

Просмотреть файл

@ -97,17 +97,15 @@ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
double x;
#endif
{
#ifdef GCC_OPT_BUG
volatile double df;
#else
fd_twoints u;
double df;
#endif
double z,p,q,r,w,s,c;
int hx,ix;
hx = __HI(x);
u.d = x;
hx = __HI(u);
ix = hx&0x7fffffff;
if(ix>=0x3ff00000) { /* |x| >= 1 */
if(((ix-0x3ff00000)|__LO(x))==0) { /* |x|==1 */
if(((ix-0x3ff00000)|__LO(u))==0) { /* |x|==1 */
if(hx>0) return 0.0; /* acos(1) = 0 */
else return pi+2.0*pio2_lo; /* acos(-1)= pi */
}
@ -131,8 +129,9 @@ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
} else { /* x > 0.5 */
z = (one-x)*0.5;
s = fd_sqrt(z);
df = s;
__LO(df) = 0;
u.d = s;
__LO(u) = 0;
df = u.d;
c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));

Просмотреть файл

@ -76,9 +76,11 @@ ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
double x;
#endif
{
fd_twoints u;
double t;
int hx;
hx = __HI(x);
u.d = x;
hx = __HI(u);
if(hx<0x3ff00000) { /* x < 1 */
return (x-x)/(x-x);
} else if(hx >=0x41b00000) { /* x > 2**28 */
@ -86,7 +88,7 @@ ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
return x+x;
} else
return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */
} else if(((hx-0x3ff00000)|__LO(x))==0) {
} else if(((hx-0x3ff00000)|__LO(u))==0) {
return 0.0; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t=x*x;

Просмотреть файл

@ -106,17 +106,15 @@ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
double x;
#endif
{
#ifdef GCC_OPT_BUG
volatile double w;
#else
double w;
#endif
double t,p,q,c,r,s;
fd_twoints u;
double w,t,p,q,c,r,s;
int hx,ix;
hx = __HI(x);
u.d = x;
hx = __HI(u);
x = u.d;
ix = hx&0x7fffffff;
if(ix>= 0x3ff00000) { /* |x|>= 1 */
if(((ix-0x3ff00000)|__LO(x))==0)
if(((ix-0x3ff00000)|__LO(u))==0)
/* asin(1)=+-pi/2 with inexact */
return x*pio2_hi+x*pio2_lo;
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
@ -140,8 +138,9 @@ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
w = p/q;
t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
} else {
w = s;
__LO(w) = 0;
u.d = s;
__LO(u) = 0;
w = u.d;
c = (t-w*w)/(s+w);
r = p/q;
p = 2.0*s*r-(pio2_lo-2.0*c);

Просмотреть файл

@ -93,14 +93,16 @@ pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
double y,x;
#endif
{
fd_twoints ux, uy, uz;
double z;
int k,m,hx,hy,ix,iy;
unsigned lx,ly;
hx = __HI(x); ix = hx&0x7fffffff;
lx = __LO(x);
hy = __HI(y); iy = hy&0x7fffffff;
ly = __LO(y);
ux.d = x; uy.d = y;
hx = __HI(ux); ix = hx&0x7fffffff;
lx = __LO(ux);
hy = __HI(uy); iy = hy&0x7fffffff;
ly = __LO(uy);
if(((ix|((lx|-(int)lx)>>31))>0x7ff00000)||
((iy|((ly|-(int)ly)>>31))>0x7ff00000)) /* x or y is NaN */
return x+y;
@ -147,7 +149,9 @@ pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
else z=fd_atan(fd_fabs(y/x)); /* safe to do y/x */
switch (m) {
case 0: return z ; /* atan(+,+) */
case 1: __HI(z) ^= 0x80000000;
case 1: uz.d = z;
__HI(uz) ^= 0x80000000;
z = uz.d;
return z ; /* atan(-,+) */
case 2: return pi-(z-pi_lo);/* atan(+,-) */
default: /* case 3 */

Просмотреть файл

@ -83,15 +83,19 @@ static double zero = 0.0;
double t;
int hx,ix;
unsigned lx;
hx = __HI(x); /* high word */
lx = __LO(x); /* low word */
fd_twoints u;
u.d = x;
hx = __HI(u); /* high word */
lx = __LO(u); /* low word */
ix = hx&0x7fffffff;
if ((ix|((lx|(-(int)lx))>>31))>0x3ff00000) /* |x|>1 */
return (x-x)/(x-x);
if(ix==0x3ff00000)
return x/zero;
if(ix<0x3e300000&&(really_big+x)>zero) return x; /* x<2**-28 */
__HI(x) = ix; /* x <- |x| */
u.d = x;
__HI(u) = ix; /* x <- |x| */
x = u.d;
if(ix<0x3fe00000) { /* x < 0.5 */
t = x+x;
t = 0.5*fd_log1p(t+t*x/(one-x));

Просмотреть файл

@ -84,12 +84,14 @@ static double one = 1.0, half=0.5, really_big = 1.0e300;
double x;
#endif
{
fd_twoints u;
double t,w;
int ix;
unsigned lx;
/* High word of |x|. */
ix = __HI(x);
u.d = x;
ix = __HI(u);
ix &= 0x7fffffff;
/* x is INF or NaN */

Просмотреть файл

@ -139,23 +139,21 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
double x;
#endif
{
#ifdef GCC_OPT_BUG
volatile int xsb;
#else
int xsb;
#endif
fd_twoints u;
double y,hi,lo,c,t;
int k;
int k, xsb;
unsigned hx;
hx = __HI(x); /* high word of x */
u.d = x;
hx = __HI(u); /* high word of x */
xsb = (hx>>31)&1; /* sign bit of x */
hx &= 0x7fffffff; /* high word of |x| */
/* filter out non-finite argument */
if(hx >= 0x40862E42) { /* if |x|>=709.78... */
if(hx>=0x7ff00000) {
if(((hx&0xfffff)|__LO(x))!=0)
u.d = x;
if(((hx&0xfffff)|__LO(u))!=0)
return x+x; /* NaN */
else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
}
@ -186,10 +184,14 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
if(k==0) return one-((x*c)/(c-2.0)-x);
else y = one-((lo-(x*c)/(2.0-c))-hi);
if(k >= -1021) {
__HI(y) += (k<<20); /* add k to y's exponent */
u.d = y;
__HI(u) += (k<<20); /* add k to y's exponent */
y = u.d;
return y;
} else {
__HI(y) += ((k+1000)<<20);/* add k to y's exponent */
u.d = y;
__HI(u) += ((k+1000)<<20);/* add k to y's exponent */
y = u.d;
return y*twom1000;
}
}

Просмотреть файл

@ -65,13 +65,15 @@ static double one = 1.0, Zero[] = {0.0, -0.0,};
double x,y ;
#endif
{
fd_twoints ux, uy;
int n,hx,hy,hz,ix,iy,sx,i;
unsigned lx,ly,lz;
hx = __HI(x); /* high word of x */
lx = __LO(x); /* low word of x */
hy = __HI(y); /* high word of y */
ly = __LO(y); /* low word of y */
ux.d = x; uy.d = y;
hx = __HI(ux); /* high word of x */
lx = __LO(ux); /* low word of x */
hy = __HI(uy); /* high word of y */
ly = __LO(uy); /* low word of y */
sx = hx&0x80000000; /* sign of x */
hx ^=sx; /* |x| */
hy &= 0x7fffffff; /* |y| */
@ -153,8 +155,10 @@ static double one = 1.0, Zero[] = {0.0, -0.0,};
}
if(iy>= -1022) { /* normalize output */
hx = ((hx-0x00100000)|((iy+1023)<<20));
__HI(x) = hx|sx;
__LO(x) = lx;
ux.d = x;
__HI(ux) = hx|sx;
__LO(ux) = lx;
x = ux.d;
} else { /* subnormal output */
n = -1022 - iy;
if(n<=20) {
@ -165,8 +169,10 @@ static double one = 1.0, Zero[] = {0.0, -0.0,};
} else {
lx = hx>>(n-32); hx = sx;
}
__HI(x) = hx|sx;
__LO(x) = lx;
ux.d = x;
__HI(ux) = hx|sx;
__LO(ux) = lx;
x = ux.d;
x *= one; /* create necessary signal */
}
return x; /* exact output */

Просмотреть файл

@ -85,33 +85,43 @@
double x, y;
#endif
{
fd_twoints ux, uy;
double a=x,b=y,t1,t2,y1,y2,w;
int j,k,ha,hb;
ha = __HI(x)&0x7fffffff; /* high word of x */
hb = __HI(y)&0x7fffffff; /* high word of y */
ux.d = x; uy.d = y;
ha = __HI(ux)&0x7fffffff; /* high word of x */
hb = __HI(uy)&0x7fffffff; /* high word of y */
if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
__HI(a) = ha; /* a <- |a| */
__HI(b) = hb; /* b <- |b| */
ux.d = a; uy.d = b;
__HI(ux) = ha; /* a <- |a| */
__HI(uy) = hb; /* b <- |b| */
a = ux.d; b = uy.d;
if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
k=0;
if(ha > 0x5f300000) { /* a>2**500 */
if(ha >= 0x7ff00000) { /* Inf or NaN */
w = a+b; /* for sNaN */
if(((ha&0xfffff)|__LO(a))==0) w = a;
if(((hb^0x7ff00000)|__LO(b))==0) w = b;
ux.d = a; uy.d = b;
if(((ha&0xfffff)|__LO(ux))==0) w = a;
if(((hb^0x7ff00000)|__LO(uy))==0) w = b;
return w;
}
/* scale a and b by 2**-600 */
ha -= 0x25800000; hb -= 0x25800000; k += 600;
__HI(a) = ha;
__HI(b) = hb;
ux.d = a; uy.d = b;
__HI(ux) = ha;
__HI(uy) = hb;
a = ux.d; b = uy.d;
}
if(hb < 0x20b00000) { /* b < 2**-500 */
if(hb <= 0x000fffff) { /* subnormal b or 0 */
if((hb|(__LO(b)))==0) return a;
uy.d = b;
if((hb|(__LO(uy)))==0) return a;
t1=0;
__HI(t1) = 0x7fd00000; /* t1=2^1022 */
ux.d = t1;
__HI(ux) = 0x7fd00000; /* t1=2^1022 */
t1 = ux.d;
b *= t1;
a *= t1;
k -= 1022;
@ -119,30 +129,40 @@
ha += 0x25800000; /* a *= 2^600 */
hb += 0x25800000; /* b *= 2^600 */
k -= 600;
__HI(a) = ha;
__HI(b) = hb;
ux.d = a; uy.d = b;
__HI(ux) = ha;
__HI(uy) = hb;
a = ux.d; b = uy.d;
}
}
/* medium size a and b */
w = a-b;
if (w>b) {
t1 = 0;
__HI(t1) = ha;
ux.d = t1;
__HI(ux) = ha;
t1 = ux.d;
t2 = a-t1;
w = fd_sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
} else {
a = a+a;
y1 = 0;
__HI(y1) = hb;
ux.d = y1;
__HI(ux) = hb;
y1 = ux.d;
y2 = b - y1;
t1 = 0;
__HI(t1) = ha+0x00100000;
ux.d = t1;
__HI(ux) = ha+0x00100000;
t1 = ux.d;
t2 = a - t1;
w = fd_sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
}
if(k!=0) {
t1 = 1.0;
__HI(t1) += (k<<20);
ux.d = t1;
__HI(ux) += (k<<20);
t1 = ux.d;
return t1*w;
} else return w;
}

Просмотреть файл

@ -125,10 +125,12 @@ static double zero = 0.0;
double x;
#endif
{
fd_twoints un;
double z, s,c,ss,cc,r,u,v;
int hx,ix;
hx = __HI(x);
un.d = x;
hx = __HI(un);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) return one/(x*x);
x = fd_fabs(x);
@ -194,12 +196,14 @@ v04 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
double x;
#endif
{
fd_twoints un;
double z, s,c,ss,cc,u,v;
int hx,ix,lx;
hx = __HI(x);
un.d = x;
hx = __HI(un);
ix = 0x7fffffff&hx;
lx = __LO(x);
lx = __LO(un);
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if(ix>=0x7ff00000) return one/(x+x*x);
if((ix|lx)==0) return -one/zero;
@ -362,9 +366,11 @@ static double pS2[5] = {
#else
double *p,*q;
#endif
fd_twoints u;
double z,r,s;
int ix;
ix = 0x7fffffff&__HI(x);
u.d = x;
ix = 0x7fffffff&__HI(u);
if(ix>=0x40200000) {p = pR8; q= pS8;}
else if(ix>=0x40122E8B){p = pR5; q= pS5;}
else if(ix>=0x4006DB6D){p = pR3; q= pS3;}
@ -497,9 +503,11 @@ static double qS2[6] = {
#else
double *p,*q;
#endif
fd_twoints u;
double s,r,z;
int ix;
ix = 0x7fffffff&__HI(x);
u.d = x;
ix = 0x7fffffff&__HI(u);
if(ix>=0x40200000) {p = qR8; q= qS8;}
else if(ix>=0x40122E8B){p = qR5; q= qS5;}
else if(ix>=0x4006DB6D){p = qR3; q= qS3;}

Просмотреть файл

@ -126,10 +126,12 @@ static double zero = 0.0;
double x;
#endif
{
fd_twoints un;
double z, s,c,ss,cc,r,u,v,y;
int hx,ix;
hx = __HI(x);
un.d = x;
hx = __HI(un);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) return one/x;
y = fd_fabs(x);
@ -195,12 +197,14 @@ static double V0[5] = {
double x;
#endif
{
fd_twoints un;
double z, s,c,ss,cc,u,v;
int hx,ix,lx;
hx = __HI(x);
un.d = x;
hx = __HI(un);
ix = 0x7fffffff&hx;
lx = __LO(x);
lx = __LO(un);
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if(ix>=0x7ff00000) return one/(x+x*x);
if((ix|lx)==0) return -one/zero;
@ -360,9 +364,11 @@ static double ps2[5] = {
#else
double *p,*q;
#endif
fd_twoints un;
double z,r,s;
int ix;
ix = 0x7fffffff&__HI(x);
un.d = x;
ix = 0x7fffffff&__HI(un);
if(ix>=0x40200000) {p = pr8; q= ps8;}
else if(ix>=0x40122E8B){p = pr5; q= ps5;}
else if(ix>=0x4006DB6D){p = pr3; q= ps3;}
@ -496,9 +502,11 @@ static double qs2[6] = {
#else
double *p,*q;
#endif
fd_twoints un;
double s,r,z;
int ix;
ix = 0x7fffffff&__HI(x);
un.d = x;
ix = 0x7fffffff&__HI(un);
if(ix>=0x40200000) {p = qr8; q= qs8;}
else if(ix>=0x40122E8B){p = qr5; q= qs5;}
else if(ix>=0x4006DB6D){p = qr3; q= qs3;}

Просмотреть файл

@ -90,6 +90,7 @@ static double zero = 0.00000000000000000000e+00;
int n; double x;
#endif
{
fd_twoints u;
int i,hx,ix,lx, sgn;
double a, b, temp, di;
double z, w;
@ -97,9 +98,10 @@ static double zero = 0.00000000000000000000e+00;
/* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
* Thus, J(-n,x) = J(n,-x)
*/
hx = __HI(x);
u.d = x;
hx = __HI(u);
ix = 0x7fffffff&hx;
lx = __LO(x);
lx = __LO(u);
/* if J(n,NaN) is NaN */
if((ix|((unsigned)(lx|-lx))>>31)>0x7ff00000) return x+x;
if(n<0){
@ -251,13 +253,15 @@ static double zero = 0.00000000000000000000e+00;
int n; double x;
#endif
{
fd_twoints u;
int i,hx,ix,lx;
int sign;
double a, b, temp;
hx = __HI(x);
u.d = x;
hx = __HI(u);
ix = 0x7fffffff&hx;
lx = __LO(x);
lx = __LO(u);
/* if Y(n,NaN) is NaN */
if((ix|((unsigned)(lx|-lx))>>31)>0x7ff00000) return x+x;
if((ix|lx)==0) return -one/zero;
@ -295,7 +299,8 @@ static double zero = 0.00000000000000000000e+00;
a = __ieee754_y0(x);
b = __ieee754_y1(x);
/* quit if b is -inf */
for(i=1;i<n&&(__HI(b) != 0xfff00000);i++){
u.d = b;
for(i=1;i<n&&(__HI(u) != 0xfff00000);i++){
temp = b;
b = ((double)(i+i)/x)*b - a;
a = temp;

Просмотреть файл

@ -195,10 +195,12 @@ static double zero= 0.00000000000000000000e+00;
double x;
#endif
{
fd_twoints u;
double y,z;
int n,ix;
ix = 0x7fffffff&__HI(x);
u.d = x;
ix = 0x7fffffff&__HI(u);
if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0);
y = -x; /* x is assume negative */
@ -217,7 +219,8 @@ static double zero= 0.00000000000000000000e+00;
y = zero; n = 0; /* y must be even */
} else {
if(ix<0x43300000) z = y+two52; /* exact */
n = __LO(z)&1; /* lower word of z */
u.d = z;
n = __LO(u)&1; /* lower word of z */
y = n;
n<<= 2;
}
@ -243,11 +246,13 @@ static double zero= 0.00000000000000000000e+00;
double x; int *signgamp;
#endif
{
fd_twoints u;
double t,y,z,nadj,p,p1,p2,p3,q,r,w;
int i,hx,lx,ix;
hx = __HI(x);
lx = __LO(x);
u.d = x;
hx = __HI(u);
lx = __LO(u);
/* purge off +-inf, NaN, +-0, and negative arguments */
*signgamp = 1;

Просмотреть файл

@ -122,12 +122,14 @@ static double zero = 0.0;
double x;
#endif
{
fd_twoints u;
double hfsq,f,s,z,R,w,t1,t2,dk;
int k,hx,i,j;
unsigned lx;
hx = __HI(x); /* high word of x */
lx = __LO(x); /* low word of x */
u.d = x;
hx = __HI(u); /* high word of x */
lx = __LO(u); /* low word of x */
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
@ -135,13 +137,16 @@ static double zero = 0.0;
return -two54/zero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
hx = __HI(x); /* high word of x */
u.d = x;
hx = __HI(u); /* high word of x */
}
if (hx >= 0x7ff00000) return x+x;
k += (hx>>20)-1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;
__HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */
u.d = x;
__HI(u) = hx|(i^0x3ff00000); /* normalize x or x/2 */
x = u.d;
k += (i>>20);
f = x-1.0;
if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */

Просмотреть файл

@ -98,12 +98,14 @@ static double zero = 0.0;
double x;
#endif
{
fd_twoints u;
double y,z;
int i,k,hx;
unsigned lx;
hx = __HI(x); /* high word of x */
lx = __LO(x); /* low word of x */
u.d = x;
hx = __HI(u); /* high word of x */
lx = __LO(u); /* low word of x */
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
@ -111,14 +113,17 @@ static double zero = 0.0;
return -two54/zero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
hx = __HI(x); /* high word of x */
u.d = x;
hx = __HI(u); /* high word of x */
}
if (hx >= 0x7ff00000) return x+x;
k += (hx>>20)-1023;
i = ((unsigned)k&0x80000000)>>31;
hx = (hx&0x000fffff)|((0x3ff-i)<<20);
y = (double)(k+i);
__HI(x) = hx;
u.d = x;
__HI(u) = hx;
x = u.d;
z = y*log10_2lo + ivln10*__ieee754_log(x);
return z+y*log10_2hi;
}

Просмотреть файл

@ -135,19 +135,17 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
double x, y;
#endif
{
#ifdef GCC_OPT_BUG
volatile double y1,t1,p_h,t,z,ax;
#else
fd_twoints ux, uy, uz;
double y1,t1,p_h,t,z,ax;
#endif
double z_h,z_l,p_l;
double t2,r,s,u,v,w;
int i,j,k,yisint,n;
int hx,hy,ix,iy;
unsigned lx,ly;
hx = __HI(x); lx = __LO(x);
hy = __HI(y); ly = __LO(y);
ux.d = x; uy.d = y;
hx = __HI(ux); lx = __LO(ux);
hy = __HI(uy); ly = __LO(uy);
ix = hx&0x7fffffff; iy = hy&0x7fffffff;
/* y==zero: x**0 = 1 */
@ -214,7 +212,9 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
} else if(yisint==1) {
#ifdef HPUX
__HI(z) ^= 1<<31; /* some HPUXes cannot negate 0.. */
uz.d = z;
__HI(uz) ^= 1<<31; /* some HPUXes cannot negate 0.. */
z = uz.d;
#else
z = -z; /* (x<0)**odd = -(|x|**odd) */
#endif
@ -243,19 +243,17 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
v = t*ivln2_l-w*ivln2;
t1 = u+v;
__LO(t1) = 0;
uz.d = t1;
__LO(uz) = 0;
t1 = uz.d;
t2 = v-(t1-u);
} else {
#ifdef GCC_OPT_BUG
volatile double s_h,t_h;
#else
double s_h,t_h;
#endif
double s2,s_l,t_l;
n = 0;
/* take care subnormal number */
if(ix<0x00100000)
{ax *= two53; n -= 53; ix = __HI(ax); }
if(ix<0x00100000)
{ax *= two53; n -= 53; uz.d = ax; ix = __HI(uz); }
n += ((ix)>>20)-0x3ff;
j = ix&0x000fffff;
/* determine interval */
@ -263,17 +261,23 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */
else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */
else {k=0;n+=1;ix -= 0x00100000;}
__HI(ax) = ix;
uz.d = ax;
__HI(uz) = ix;
ax = uz.d;
/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
v = one/(ax+bp[k]);
s = u*v;
s_h = s;
__LO(s_h) = 0;
uz.d = s_h;
__LO(uz) = 0;
s_h = uz.d;
/* t_h=ax+bp[k] High */
t_h = zero;
__HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18);
uz.d = t_h;
__HI(uz)=((ix>>1)|0x20000000)+0x00080000+(k<<18);
t_h = uz.d;
t_l = ax - (t_h-bp[k]);
s_l = v*((u-s_h*t_h)-s_h*t_l);
/* compute log(ax) */
@ -282,21 +286,27 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
r += s_l*(s_h+s);
s2 = s_h*s_h;
t_h = 3.0+s2+r;
__LO(t_h) = 0;
uz.d = t_h;
__LO(uz) = 0;
t_h = uz.d;
t_l = r-((t_h-3.0)-s2);
/* u+v = s*(1+...) */
u = s_h*t_h;
v = s_l*t_h+t_l*s;
/* 2/(3log2)*(s+...) */
p_h = u+v;
__LO(p_h) = 0;
uz.d = p_h;
__LO(uz) = 0;
p_h = uz.d;
p_l = v-(p_h-u);
z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
z_l = cp_l*p_h+p_l*cp+dp_l[k];
/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
t = (double)n;
t1 = (((z_h+z_l)+dp_h[k])+t);
__LO(t1) = 0;
uz.d = t1;
__LO(uz) = 0;
t1 = uz.d;
t2 = z_l-(((t1-t)-dp_h[k])-z_h);
}
@ -305,12 +315,15 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
/* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
y1 = y;
__LO(y1) = 0;
uy.d = y1;
__LO(uy) = 0;
y1 = uy.d;
p_l = (y-y1)*t1+y*t2;
p_h = y1*t1;
z = p_l+p_h;
j = __HI(z);
i = __LO(z);
uz.d = z;
j = __HI(uz);
i = __LO(uz);
if (j>=0x40900000) { /* z >= 1024 */
if(((j-0x40900000)|i)!=0) /* if z > 1024 */
@ -335,13 +348,17 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
n = j+(0x00100000>>(k+1));
k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
t = zero;
__HI(t) = (n&~(0x000fffff>>k));
uz.d = t;
__HI(uz) = (n&~(0x000fffff>>k));
t = uz.d;
n = ((n&0x000fffff)|0x00100000)>>(20-k);
if(j<0) n = -n;
p_h -= t;
}
t = p_l+p_h;
__LO(t) = 0;
uz.d = t;
__LO(uz) = 0;
t = uz.d;
u = t*lg2_h;
v = (p_l-(t-p_h))*lg2+t*lg2_l;
z = u+v;
@ -350,9 +367,10 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
r = (z*t1)/(t1-two)-(w+z*w);
z = one-(r-z);
j = __HI(z);
uz.d = z;
j = __HI(uz);
j += (n<<20);
if((j>>20)<=0) z = fd_scalbn(z,n); /* subnormal output */
else __HI(z) += (n<<20);
else { uz.d = z; __HI(uz) += (n<<20); z = uz.d; }
return s*z;
}

Просмотреть файл

@ -120,11 +120,13 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
double x,y[];
#endif
{
fd_twoints u, ux, uz;
double z,w,t,r,fn;
double tx[3];
int e0,i,j,nx,n,ix,hx;
hx = __HI(x); /* high word of x */
u.d = x;
hx = __HI(u); /* high word of x */
ix = hx&0x7fffffff;
if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */
{y[0] = x; y[1] = 0; return 0;}
@ -163,15 +165,17 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
y[0] = r-w; /* quick check no cancellation */
} else {
j = ix>>20;
y[0] = r-w;
i = j-(((__HI(y[0]))>>20)&0x7ff);
y[0] = r-w;
u.d = y[0];
i = j-(((__HI(u))>>20)&0x7ff);
if(i>16) { /* 2nd iteration needed, good to 118 */
t = r;
w = fn*pio2_2;
r = t-w;
w = fn*pio2_2t-((t-r)-w);
y[0] = r-w;
i = j-(((__HI(y[0]))>>20)&0x7ff);
u.d = y[0];
i = j-(((__HI(u))>>20)&0x7ff);
if(i>49) { /* 3rd iteration need, 151 bits acc */
t = r; /* will cover all possible cases */
w = fn*pio2_3;
@ -192,9 +196,13 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
y[0]=y[1]=x-x; return 0;
}
/* set z = scalbn(|x|,ilogb(x)-23) */
__LO(z) = __LO(x);
ux.d = x; uz.d = z;
__LO(uz) = __LO(ux);
z = uz.d;
e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */
__HI(z) = ix - (e0<<20);
uz.d = z;
__HI(uz) = ix - (e0<<20);
z = uz.d;
for(i=0;i<2;i++) {
tx[i] = (double)((int)(z));
z = (z-tx[i])*two24;

Просмотреть файл

@ -69,14 +69,17 @@ static double zero = 0.0;
double x,p;
#endif
{
fd_twoints u;
int hx,hp;
unsigned sx,lx,lp;
double p_half;
hx = __HI(x); /* high word of x */
lx = __LO(x); /* low word of x */
hp = __HI(p); /* high word of p */
lp = __LO(p); /* low word of p */
u.d = x;
hx = __HI(u); /* high word of x */
lx = __LO(u); /* low word of x */
u.d = p;
hp = __HI(u); /* high word of p */
lp = __LO(u); /* low word of p */
sx = hx&0x80000000;
hp &= 0x7fffffff;
hx &= 0x7fffffff;
@ -105,6 +108,8 @@ static double zero = 0.0;
if(x>=p_half) x -= p;
}
}
__HI(x) ^= sx;
u.d = x;
__HI(u) ^= sx;
x = u.d;
return x;
}

Просмотреть файл

@ -77,12 +77,14 @@ static double one = 1.0, shuge = 1.0e307;
double x;
#endif
{
fd_twoints u;
double t,w,h;
int ix,jx;
unsigned lx;
/* High word of |x|. */
jx = __HI(x);
u.d = x;
jx = __HI(u);
ix = jx&0x7fffffff;
/* x is INF or NaN */

Просмотреть файл

@ -128,13 +128,15 @@ static double one = 1.0, tiny=1.0e-300;
double x;
#endif
{
fd_twoints u;
double z;
int sign = (int)0x80000000;
unsigned r,t1,s1,ix1,q1;
int ix0,s0,q,m,t,i;
ix0 = __HI(x); /* high word of x */
ix1 = __LO(x); /* low word of x */
u.d = x;
ix0 = __HI(u); /* high word of x */
ix1 = __LO(u); /* low word of x */
/* take care of Inf and NaN */
if((ix0&0x7ff00000)==0x7ff00000) {
@ -219,8 +221,10 @@ static double one = 1.0, tiny=1.0e-300;
ix1 = q1>>1;
if ((q&1)==1) ix1 |= sign;
ix0 += (m <<20);
__HI(z) = ix0;
__LO(z) = ix1;
u.d = z;
__HI(u) = ix0;
__LO(u) = ix1;
z = u.d;
return z;
}

Просмотреть файл

@ -45,12 +45,14 @@
*/
/* Modified defines start here.. */
#undef __LITTLE_ENDIAN
#ifdef _WIN32
#define huge myhuge
#define __LITTLE_ENDIAN
#endif
#ifdef X86_LINUX
#if defined(linux) && defined(__i386__)
#define __LITTLE_ENDIAN
#endif
@ -63,19 +65,19 @@
#endif
#endif
typedef union {
#ifdef __LITTLE_ENDIAN
#define __HI(x) *(1+(int*)&x)
#define __LO(x) *(int*)&x
#define __HIp(x) *(1+(int*)x)
#define __LOp(x) *(int*)x
struct { int lo, hi; } ints;
#else
#define __HI(x) *(int*)&x
#define __LO(x) *(1+(int*)&x)
#define __HIp(x) *(int*)x
#define __LOp(x) *(1+(int*)x)
struct { int hi, lo; } ints;
#endif
double d;
} fd_twoints;
#define __HI(x) x.ints.hi
#define __LO(x) x.ints.lo
#undef __P
#ifdef __STDC__
#define __P(p) p
#else

Просмотреть файл

@ -101,9 +101,11 @@ C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
double x,y;
#endif
{
fd_twoints u;
double a,hz,z,r,qx;
int ix;
ix = __HI(x)&0x7fffffff; /* ix = |x|'s high word*/
u.d = x;
ix = __HI(u)&0x7fffffff; /* ix = |x|'s high word*/
if(ix<0x3e400000) { /* if x < 2**27 */
if(((int)x)==0) return one; /* generate inexact */
}
@ -115,8 +117,10 @@ C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
if(ix > 0x3fe90000) { /* x > 0.78125 */
qx = 0.28125;
} else {
__HI(qx) = ix-0x00200000; /* x/4 */
__LO(qx) = 0;
u.d = qx;
__HI(u) = ix-0x00200000; /* x/4 */
__LO(u) = 0;
qx = u.d;
}
hz = 0.5*z-qx;
a = one-qx;

Просмотреть файл

@ -94,9 +94,11 @@ S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
double x,y; int iy; /* iy=0 if y is zero */
#endif
{
fd_twoints u;
double z,r,v;
int ix;
ix = __HI(x)&0x7fffffff; /* high word of x */
u.d = x;
ix = __HI(u)&0x7fffffff; /* high word of x */
if(ix<0x3e400000) /* |x| < 2**-27 */
{if((int)x==0) return x;} /* generate inexact */
z = x*x;

Просмотреть файл

@ -127,8 +127,11 @@ static double zero = 0.0; /* used as const */
#ifndef HUGE_VAL /* this is the only routine that uses HUGE_VAL */
#define HUGE_VAL inf
double inf = 0.0;
fd_twoints u;
__HI(inf) = 0x7ff00000; /* set inf to infinite */
u.d = inf;
__HI(u) = 0x7ff00000; /* set inf to infinite */
inf = u.d;
#endif
#ifdef _USE_WRITE

Просмотреть файл

@ -110,13 +110,16 @@ T[] = {
double x,y; int iy;
#endif
{
fd_twoints u;
double z,r,v,w,s;
int ix,hx;
hx = __HI(x); /* high word of x */
u.d = x;
hx = __HI(u); /* high word of x */
ix = hx&0x7fffffff; /* high word of |x| */
if(ix<0x3e300000) /* x < 2**-28 */
{if((int)x==0) { /* generate inexact */
if(((ix|__LO(x))|(iy+1))==0) return one/fd_fabs(x);
u.d =x;
if(((ix|__LO(u))|(iy+1))==0) return one/fd_fabs(x);
else return (iy==1)? x: -one/x;
}
}
@ -148,10 +151,14 @@ T[] = {
/* compute -1.0/(x+r) accurately */
double a,t;
z = w;
__LO(z) = 0;
u.d = z;
__LO(u) = 0;
z = u.d;
v = r-(z - x); /* z+v = r+x */
t = a = -1.0/w; /* a = -1.0/w */
__LO(t) = 0;
u.d = t;
__LO(u) = 0;
t = u.d;
s = 1.0+t*z;
return t+a*(s+t*v);
}

Просмотреть файл

@ -73,9 +73,11 @@ really_big= 1.00000000000000000000e+300;
double x;
#endif
{
fd_twoints u;
double t,w;
int hx,ix;
hx = __HI(x);
u.d = x;
hx = __HI(u);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */
if(ix< 0x3e300000) { /* |x|<2**-28 */

Просмотреть файл

@ -122,14 +122,17 @@ really_big = 1.0e300;
double x;
#endif
{
fd_twoints u;
double w,s1,s2,z;
int ix,hx,id;
hx = __HI(x);
u.d = x;
hx = __HI(u);
ix = hx&0x7fffffff;
if(ix>=0x44100000) { /* if |x| >= 2^66 */
u.d = x;
if(ix>0x7ff00000||
(ix==0x7ff00000&&(__LO(x)!=0)))
(ix==0x7ff00000&&(__LO(u)!=0)))
return x+x; /* NaN */
if(hx>0) return atanhi[3]+atanlo[3];
else return -atanhi[3]-atanlo[3];

Просмотреть файл

@ -76,26 +76,31 @@ G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */
double x;
#endif
{
fd_twoints u;
int hx;
double r,s,t=0.0,w;
unsigned sign;
hx = __HI(x); /* high word of x */
u.d = x;
hx = __HI(u); /* high word of x */
sign=hx&0x80000000; /* sign= sign(x) */
hx ^=sign;
if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
if((hx|__LO(x))==0)
if((hx|__LO(u))==0) {
x = u.d;
return(x); /* cbrt(0) is itself */
__HI(x) = hx; /* x <- |x| */
}
u.d = x;
__HI(u) = hx; /* x <- |x| */
x = u.d;
/* rough cbrt to 5 bits */
if(hx<0x00100000) /* subnormal number */
{__HI(t)=0x43500000; /* set t= 2**54 */
t*=x; __HI(t)=__HI(t)/3+B2;
{u.d = t; __HI(u)=0x43500000; t=u.d; /* set t= 2**54 */
t*=x; __HI(u)=__HI(u)/3+B2;
}
else
__HI(t)=hx/3+B1;
else {
u.d = t; __HI(u)=hx/3+B1; t = u.d;
}
/* new cbrt to 23 bits, may be implemented in single precision */
@ -104,8 +109,9 @@ G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */
t*=G+F/(s+E+D/s);
/* chopped to 20 bits and make it larger than cbrt(x) */
__LO(t)=0; __HI(t)+=0x00000001;
u.d = t;
__LO(u)=0; __HI(u)+=0x00000001;
t = u.d;
/* one step newton iteration to 53 bits with error less than 0.667 ulps */
s=t*t; /* t*t is exact */
@ -115,6 +121,8 @@ G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */
t=t+t*r;
/* retore the sign bit */
__HI(t) |= sign;
u.d = t;
__HI(u) |= sign;
t = u.d;
return(t);
}

Просмотреть файл

@ -68,10 +68,12 @@ static double really_big = 1.0e300;
double x;
#endif
{
fd_twoints u;
int i0,i1,j0;
unsigned i,j;
i0 = __HI(x);
i1 = __LO(x);
u.d = x;
i0 = __HI(u);
i1 = __LO(u);
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<20) {
if(j0<0) { /* raise inexact if x != 0 */
@ -105,7 +107,9 @@ static double really_big = 1.0e300;
i1 &= (~i);
}
}
__HI(x) = i0;
__LO(x) = i1;
u.d = x;
__HI(u) = i0;
__LO(u) = i1;
x = u.d;
return x;
}

Просмотреть файл

@ -59,6 +59,9 @@
double x,y;
#endif
{
__HI(x) = (__HI(x)&0x7fffffff)|(__HI(y)&0x80000000);
fd_twoints ux, uy;
ux.d = x; uy.d = y;
__HI(ux) = (__HI(ux)&0x7fffffff)|(__HI(uy)&0x80000000);
x = ux.d;
return x;
}

Просмотреть файл

@ -84,11 +84,13 @@
double x;
#endif
{
fd_twoints u;
double y[2],z=0.0;
int n, ix;
/* High word of x. */
ix = __HI(x);
u.d = x;
ix = __HI(u);
/* |x| ~< pi/4 */
ix &= 0x7fffffff;

Просмотреть файл

@ -227,9 +227,11 @@ sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
double x;
#endif
{
fd_twoints u;
int hx,ix,i;
double R,S,P,Q,s,y,z,r;
hx = __HI(x);
u.d = x;
hx = __HI(u);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) { /* erf(nan)=nan */
i = ((unsigned)hx>>31)<<1;
@ -271,7 +273,9 @@ sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
sb5+s*(sb6+s*sb7))))));
}
z = x;
__LO(z) = 0;
u.d = z;
__LO(u) = 0;
z = u.d;
r = __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
if(hx>=0) return one-r/x; else return r/x-one;
}
@ -283,9 +287,11 @@ sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
double x;
#endif
{
fd_twoints u;
int hx,ix;
double R,S,P,Q,s,y,z,r;
hx = __HI(x);
u.d = x;
hx = __HI(u);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) { /* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */
@ -333,7 +339,9 @@ sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
sb5+s*(sb6+s*sb7))))));
}
z = x;
__LO(z) = 0;
u.d = z;
__LO(u) = 0;
z = u.d;
r = __ieee754_exp(-z*z-0.5625)*
__ieee754_exp((z-x)*(z+x)+R/S);
if(hx>0) return r/x; else return two-r/x;

Просмотреть файл

@ -167,11 +167,13 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
double x;
#endif
{
fd_twoints u;
double y,hi,lo,c,t,e,hxs,hfx,r1;
int k,xsb;
unsigned hx;
hx = __HI(x); /* high word of x */
u.d = x;
hx = __HI(u); /* high word of x */
xsb = hx&0x80000000; /* sign bit of x */
if(xsb==0) y=x; else y= -x; /* y = |x| */
hx &= 0x7fffffff; /* high word of |x| */
@ -180,7 +182,8 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */
if(hx >= 0x40862E42) { /* if |x|>=709.78... */
if(hx>=0x7ff00000) {
if(((hx&0xfffff)|__LO(x))!=0)
u.d = x;
if(((hx&0xfffff)|__LO(u))!=0)
return x+x; /* NaN */
else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
}
@ -230,19 +233,29 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
else return one+2.0*(x-e);
if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
y = one-(e-x);
__HI(y) += (k<<20); /* add k to y's exponent */
u.d = y;
__HI(u) += (k<<20); /* add k to y's exponent */
y = u.d;
return y-one;
}
t = one;
if(k<20) {
__HI(t) = 0x3ff00000 - (0x200000>>k); /* t=1-2^-k */
u.d = t;
__HI(u) = 0x3ff00000 - (0x200000>>k); /* t=1-2^-k */
t = u.d;
y = t-(e-x);
__HI(y) += (k<<20); /* add k to y's exponent */
u.d = y;
__HI(u) += (k<<20); /* add k to y's exponent */
y = u.d;
} else {
__HI(t) = ((0x3ff-k)<<20); /* 2^-k */
u.d = t;
__HI(u) = ((0x3ff-k)<<20); /* 2^-k */
t = u.d;
y = x-(e+t);
y += one;
__HI(y) += (k<<20); /* add k to y's exponent */
u.d = y;
__HI(u) += (k<<20); /* add k to y's exponent */
y = u.d;
}
}
return y;

Просмотреть файл

@ -57,6 +57,9 @@
double x;
#endif
{
__HI(x) &= 0x7fffffff;
fd_twoints u;
u.d = x;
__HI(u) &= 0x7fffffff;
x = u.d;
return x;
}

Просмотреть файл

@ -58,7 +58,9 @@
double x;
#endif
{
fd_twoints u;
int hx;
hx = __HI(x);
u.d = x;
hx = __HI(u);
return (unsigned)((hx&0x7fffffff)-0x7ff00000)>>31;
}

Просмотреть файл

@ -68,10 +68,12 @@ static double really_big = 1.0e300;
double x;
#endif
{
fd_twoints u;
int i0,i1,j0;
unsigned i,j;
i0 = __HI(x);
i1 = __LO(x);
u.d = x;
i0 = __HI(u);
i1 = __LO(u);
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<20) {
if(j0<0) { /* raise inexact if x != 0 */
@ -106,7 +108,9 @@ static double really_big = 1.0e300;
i1 &= (~i);
}
}
__HI(x) = i0;
__LO(x) = i1;
u.d = x;
__HI(u) = i0;
__LO(u) = i1;
x = u.d;
return x;
}

Просмотреть файл

@ -71,19 +71,24 @@ two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
#endif
{
int hx, ix, lx;
hx = __HI(x);
fd_twoints u;
u.d = x;
hx = __HI(u);
ix = 0x7fffffff&hx;
lx = __LO(x);
lx = __LO(u);
*eptr = 0;
if(ix>=0x7ff00000||((ix|lx)==0)) return x; /* 0,inf,nan */
if (ix<0x00100000) { /* subnormal */
x *= two54;
hx = __HI(x);
u.d = x;
hx = __HI(u);
ix = hx&0x7fffffff;
*eptr = -54;
}
*eptr += (ix>>20)-1022;
hx = (hx&0x800fffff)|0x3fe00000;
__HI(x) = hx;
u.d = x;
__HI(u) = hx;
x = u.d;
return x;
}

Просмотреть файл

@ -60,10 +60,11 @@
#endif
{
int hx,lx,ix;
hx = (__HI(x))&0x7fffffff; /* high word of x */
fd_twoints u;
u.d = x;
hx = (__HI(u))&0x7fffffff; /* high word of x */
if(hx<0x00100000) {
lx = __LO(x);
lx = __LO(u);
if((hx|lx)==0)
return 0x80000001; /* ilogb(0) = 0x80000001 */
else /* subnormal x */

Просмотреть файл

@ -58,9 +58,11 @@
double x;
#endif
{
fd_twoints u;
int hx,lx;
hx = (__HI(x)&0x7fffffff);
lx = __LO(x);
u.d = x;
hx = (__HI(u)&0x7fffffff);
lx = __LO(u);
hx |= (unsigned)(lx|(-lx))>>31;
hx = 0x7ff00000 - hx;
return ((unsigned)(hx))>>31;

Просмотреть файл

@ -138,8 +138,10 @@ static double zero = 0.0;
{
double hfsq,f,c,s,z,R,u;
int k,hx,hu,ax;
fd_twoints un;
hx = __HI(x); /* high word of x */
un.d = x;
hx = __HI(un); /* high word of x */
ax = hx&0x7fffffff;
k = 1;
@ -162,22 +164,28 @@ static double zero = 0.0;
if(k!=0) {
if(hx<0x43400000) {
u = 1.0+x;
hu = __HI(u); /* high word of u */
un.d = u;
hu = __HI(un); /* high word of u */
k = (hu>>20)-1023;
c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */
c /= u;
} else {
u = x;
hu = __HI(u); /* high word of u */
un.d = u;
hu = __HI(un); /* high word of u */
k = (hu>>20)-1023;
c = 0;
}
hu &= 0x000fffff;
if(hu<0x6a09e) {
__HI(u) = hu|0x3ff00000; /* normalize u */
un.d = u;
__HI(un) = hu|0x3ff00000; /* normalize u */
u = un.d;
} else {
k += 1;
__HI(u) = hu|0x3fe00000; /* normalize u/2 */
un.d = u;
__HI(un) = hu|0x3fe00000; /* normalize u/2 */
u = un.d;
hu = (0x00100000-hu)>>2;
}
f = u-1.0;

Просмотреть файл

@ -60,8 +60,11 @@
#endif
{
int lx,ix;
ix = (__HI(x))&0x7fffffff; /* high |x| */
lx = __LO(x); /* low x */
fd_twoints u;
u.d = x;
ix = (__HI(u))&0x7fffffff; /* high |x| */
lx = __LO(u); /* low x */
if((ix|lx)==0) return -1.0/fd_fabs(x);
if(ix>=0x7ff00000) return x*x;
if((ix>>=20)==0) /* IEEE 754 logb */

Просмотреть файл

@ -71,42 +71,56 @@ static double one = 1.0;
{
int i0,i1,j0;
unsigned i;
i0 = __HI(x); /* high x */
i1 = __LO(x); /* low x */
fd_twoints u;
u.d = x;
i0 = __HI(u); /* high x */
i1 = __LO(u); /* low x */
j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent of x */
if(j0<20) { /* integer part in high x */
if(j0<0) { /* |x|<1 */
__HIp(iptr) = i0&0x80000000;
__LOp(iptr) = 0; /* *iptr = +-0 */
u.d = *iptr;
__HI(u) = i0&0x80000000;
__LO(u) = 0; /* *iptr = +-0 */
*iptr = u.d;
return x;
} else {
i = (0x000fffff)>>j0;
if(((i0&i)|i1)==0) { /* x is integral */
*iptr = x;
__HI(x) &= 0x80000000;
__LO(x) = 0; /* return +-0 */
u.d = x;
__HI(u) &= 0x80000000;
__LO(u) = 0; /* return +-0 */
x = u.d;
return x;
} else {
__HIp(iptr) = i0&(~i);
__LOp(iptr) = 0;
u.d = *iptr;
__HI(u) = i0&(~i);
__LO(u) = 0;
*iptr = u.d;
return x - *iptr;
}
}
} else if (j0>51) { /* no fraction part */
*iptr = x*one;
__HI(x) &= 0x80000000;
__LO(x) = 0; /* return +-0 */
u.d = x;
__HI(u) &= 0x80000000;
__LO(u) = 0; /* return +-0 */
x = u.d;
return x;
} else { /* fraction part in low x */
i = ((unsigned)(0xffffffff))>>(j0-20);
if((i1&i)==0) { /* x is integral */
*iptr = x;
__HI(x) &= 0x80000000;
__LO(x) = 0; /* return +-0 */
u.d = x;
__HI(u) &= 0x80000000;
__LO(u) = 0; /* return +-0 */
x = u.d;
return x;
} else {
__HIp(iptr) = i0;
__LOp(iptr) = i1&(~i);
u.d = *iptr;
__HI(u) = i0;
__LO(u) = i1&(~i);
*iptr = u.d;
return x - *iptr;
}
}

Просмотреть файл

@ -62,11 +62,13 @@
{
int hx,hy,ix,iy;
unsigned lx,ly;
fd_twoints ux, uy;
hx = __HI(x); /* high word of x */
lx = __LO(x); /* low word of x */
hy = __HI(y); /* high word of y */
ly = __LO(y); /* low word of y */
ux.d = x; uy.d = y;
hx = __HI(ux); /* high word of x */
lx = __LO(ux); /* low word of x */
hy = __HI(uy); /* high word of y */
ly = __LO(uy); /* low word of y */
ix = hx&0x7fffffff; /* |x| */
iy = hy&0x7fffffff; /* |y| */
@ -75,8 +77,10 @@
return x+y;
if(x==y) return x; /* x=y, return x */
if((ix|lx)==0) { /* x == 0 */
__HI(x) = hy&0x80000000; /* return +-minsubnormal */
__LO(x) = 1;
ux.d = x;
__HI(ux) = hy&0x80000000; /* return +-minsubnormal */
__LO(ux) = 1;
x = ux.d;
y = x*x;
if(y==x) return y; else return x; /* raise underflow flag */
}
@ -102,10 +106,14 @@
if(hy<0x00100000) { /* underflow */
y = x*x;
if(y!=x) { /* raise underflow flag */
__HI(y) = hx; __LO(y) = lx;
uy.d = y;
__HI(uy) = hx; __LO(uy) = lx;
y = uy.d;
return y;
}
}
__HI(x) = hx; __LO(x) = lx;
ux.d = x;
__HI(ux) = hx; __LO(ux) = lx;
x = ux.d;
return x;
}

Просмотреть файл

@ -76,9 +76,12 @@ TWO52[2]={
int i0,j0,sx;
unsigned i,i1;
double w,t;
i0 = __HI(x);
fd_twoints u;
u.d = x;
i0 = __HI(u);
sx = (i0>>31)&1;
i1 = __LO(x);
i1 = __LO(u);
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<20) {
if(j0<0) {
@ -86,11 +89,15 @@ TWO52[2]={
i1 |= (i0&0x0fffff);
i0 &= 0xfffe0000;
i0 |= ((i1|-(int)i1)>>12)&0x80000;
__HI(x)=i0;
u.d = x;
__HI(u)=i0;
x = u.d;
w = TWO52[sx]+x;
t = w-TWO52[sx];
i0 = __HI(t);
__HI(t) = (i0&0x7fffffff)|(sx<<31);
u.d = t;
i0 = __HI(u);
__HI(u) = (i0&0x7fffffff)|(sx<<31);
t = u.d;
return t;
} else {
i = (0x000fffff)>>j0;
@ -110,8 +117,10 @@ TWO52[2]={
i>>=1;
if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
}
__HI(x) = i0;
__LO(x) = i1;
u.d = x;
__HI(u) = i0;
__LO(u) = i1;
x = u.d;
w = TWO52[sx]+x;
return w-TWO52[sx];
}

Просмотреть файл

@ -70,14 +70,17 @@ tiny = 1.0e-300;
double x; int n;
#endif
{
fd_twoints u;
int k,hx,lx;
hx = __HI(x);
lx = __LO(x);
u.d = x;
hx = __HI(u);
lx = __LO(u);
k = (hx&0x7ff00000)>>20; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
x *= two54;
hx = __HI(x);
x *= two54;
u.d = x;
hx = __HI(u);
k = ((hx&0x7ff00000)>>20) - 54;
if (n< -50000) return tiny*x; /*underflow*/
}
@ -85,13 +88,15 @@ tiny = 1.0e-300;
k = k+n;
if (k > 0x7fe) return really_big*fd_copysign(really_big,x); /* overflow */
if (k > 0) /* normal result */
{__HI(x) = (hx&0x800fffff)|(k<<20); return x;}
{u.d = x; __HI(u) = (hx&0x800fffff)|(k<<20); x = u.d; return x;}
if (k <= -54) {
if (n > 50000) /* in case integer overflow in n+k */
return really_big*fd_copysign(really_big,x); /*overflow*/
else return tiny*fd_copysign(tiny,x); /*underflow*/
}
k += 54; /* subnormal result */
__HI(x) = (hx&0x800fffff)|(k<<20);
u.d = x;
__HI(u) = (hx&0x800fffff)|(k<<20);
x = u.d;
return x*twom54;
}

Просмотреть файл

@ -84,11 +84,13 @@
double x;
#endif
{
fd_twoints u;
double y[2],z=0.0;
int n, ix;
/* High word of x. */
ix = __HI(x);
u.d = x;
ix = __HI(u);
/* |x| ~< pi/4 */
ix &= 0x7fffffff;

Просмотреть файл

@ -83,11 +83,13 @@
double x;
#endif
{
fd_twoints u;
double y[2],z=0.0;
int n, ix;
/* High word of x. */
ix = __HI(x);
u.d = x;
ix = __HI(u);
/* |x| ~< pi/4 */
ix &= 0x7fffffff;

Просмотреть файл

@ -85,9 +85,11 @@ static double one=1.0, two=2.0, tiny = 1.0e-300;
{
double t,z;
int jx,ix;
fd_twoints u;
/* High word of |x|. */
jx = __HI(x);
u.d = x;
jx = __HI(u);
ix = jx&0x7fffffff;
/* x is INF or NaN */