Skip to content

Commit

Permalink
fixing a batch of fdlibm functions
Browse files Browse the repository at this point in the history
  • Loading branch information
@JeffBezanson
JeffBezanson committed Jun 16, 2011
1 parent 6c27dde commit b248040
Show file tree
Hide file tree
Showing 16 changed files with 233 additions and 93 deletions.
8 changes: 5 additions & 3 deletions external/fdlibm/e_asin.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,10 +75,12 @@ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
{
double t,w,p,q,c,r,s;
int hx,ix;
hx = __HI(x);
GET_HIGH_WORD(hx, x);
ix = hx&0x7fffffff;
if(ix>= 0x3ff00000) { /* |x|>= 1 */
if(((ix-0x3ff00000)|__LO(x))==0)
int _lx;
GET_LOW_WORD(_lx, x);
if(((ix-0x3ff00000)|_lx)==0)
/* asin(1)=+-pi/2 with inexact */
return x*pio2_hi+x*pio2_lo;
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
Expand All @@ -103,7 +105,7 @@ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
} else {
w = s;
__LO(w) = 0;
SET_LOW_WORD(w, 0);
c = (t-w*w)/(s+w);
r = p/q;
p = 2.0*s*r-(pio2_lo-2.0*c);
Expand Down
2 changes: 1 addition & 1 deletion external/fdlibm/e_cosh.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ static double one = 1.0, half=0.5, huge = 1.0e300;
unsigned lx;

/* High word of |x|. */
ix = __HI(x);
GET_HIGH_WORD(ix, x);
ix &= 0x7fffffff;

/* x is INF or NaN */
Expand Down
39 changes: 23 additions & 16 deletions external/fdlibm/e_hypot.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,61 +55,68 @@
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 */
GET_HIGH_WORD(ha, x); ha &= 0x7fffffff; /* high word of x */
GET_HIGH_WORD(hb, y); hb &= 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| */
SET_HIGH_WORD(a, ha); /* a <- |a| */
SET_HIGH_WORD(b, hb); /* b <- |b| */
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;
int _la, _lb;
GET_LOW_WORD(_la, a);
GET_LOW_WORD(_lb, b);
if(((ha&0xfffff)|_la)==0) w = a;
if(((hb^0x7ff00000)|_lb)==0) w = b;
return w;
}
/* scale a and b by 2**-600 */
ha -= 0x25800000; hb -= 0x25800000; k += 600;
__HI(a) = ha;
__HI(b) = hb;
SET_HIGH_WORD(a, ha);
SET_HIGH_WORD(b, hb);
}
if(hb < 0x20b00000) { /* b < 2**-500 */
if(hb <= 0x000fffff) { /* subnormal b or 0 */
if((hb|(__LO(b)))==0) return a;
int _lb;
GET_LOW_WORD(_lb, b);
if((hb|(_lb))==0) return a;
t1=0;
__HI(t1) = 0x7fd00000; /* t1=2^1022 */
SET_HIGH_WORD(t1, 0x7fd00000); /* t1=2^1022 */
b *= t1;
a *= t1;
k -= 1022;
} else { /* scale a and b by 2^600 */
ha += 0x25800000; /* a *= 2^600 */
hb += 0x25800000; /* b *= 2^600 */
k -= 600;
__HI(a) = ha;
__HI(b) = hb;
SET_HIGH_WORD(a, ha);
SET_HIGH_WORD(b, hb);
}
}
/* medium size a and b */
w = a-b;
if (w>b) {
t1 = 0;
__HI(t1) = ha;
SET_HIGH_WORD(t1, ha);
t2 = a-t1;
w = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
} else {
a = a+a;
y1 = 0;
__HI(y1) = hb;
SET_HIGH_WORD(y1, hb);
y2 = b - y1;
t1 = 0;
__HI(t1) = ha+0x00100000;
SET_HIGH_WORD(t1, ha+0x00100000);
t2 = a - t1;
w = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
}
if(k!=0) {
t1 = 1.0;
__HI(t1) += (k<<20);
int _ht1;
GET_HIGH_WORD(_ht1, t1);
SET_HIGH_WORD(t1, _ht1 + (k<<20));
return t1*w;
} else return w;
}
40 changes: 23 additions & 17 deletions external/fdlibm/e_pow.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,8 +111,10 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
unsigned lx,ly;

i0 = ((*(int*)&one)>>29)^1; i1=1-i0;
hx = __HI(x); lx = __LO(x);
hy = __HI(y); ly = __LO(y);
GET_HIGH_WORD(hx, x); /* high word of x */
GET_LOW_WORD(lx, x); /* low word of x */
GET_HIGH_WORD(hy, y); /* high word of y */
GET_LOW_WORD(ly, y); /* low word of y */
ix = hx&0x7fffffff; iy = hy&0x7fffffff;

/* y==zero: x**0 = 1 */
Expand Down Expand Up @@ -203,32 +205,32 @@ 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;
SET_LOW_WORD(t1, 0);
t2 = v-(t1-u);
} else {
double ss,s2,s_h,s_l,t_h,t_l;
n = 0;
/* take care subnormal number */
if(ix<0x00100000)
{ax *= two53; n -= 53; ix = __HI(ax); }
{ax *= two53; n -= 53; GET_HIGH_WORD(ix, ax); }
n += ((ix)>>20)-0x3ff;
j = ix&0x000fffff;
/* determine interval */
ix = j|0x3ff00000; /* normalize ix */
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;
SET_HIGH_WORD(ax, ix);

/* compute ss = 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]);
ss = u*v;
s_h = ss;
__LO(s_h) = 0;
SET_LOW_WORD(s_h, 0);
/* t_h=ax+bp[k] High */
t_h = zero;
__HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18);
SET_HIGH_WORD(t_h, ((ix>>1)|0x20000000)+0x00080000+(k<<18));
t_l = ax - (t_h-bp[k]);
s_l = v*((u-s_h*t_h)-s_h*t_l);
/* compute log(ax) */
Expand All @@ -237,32 +239,32 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
r += s_l*(s_h+ss);
s2 = s_h*s_h;
t_h = 3.0+s2+r;
__LO(t_h) = 0;
SET_LOW_WORD(t_h, 0);
t_l = r-((t_h-3.0)-s2);
/* u+v = ss*(1+...) */
u = s_h*t_h;
v = s_l*t_h+t_l*ss;
/* 2/(3log2)*(ss+...) */
p_h = u+v;
__LO(p_h) = 0;
SET_LOW_WORD(p_h, 0);
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) = (ss+..)*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;
SET_LOW_WORD(t1, 0);
t2 = z_l-(((t1-t)-dp_h[k])-z_h);
}

/* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
y1 = y;
__LO(y1) = 0;
SET_LOW_WORD(y1, 0);
p_l = (y-y1)*t1+y*t2;
p_h = y1*t1;
z = p_l+p_h;
j = __HI(z);
i = __LO(z);
GET_HIGH_WORD(j, z);
GET_LOW_WORD(i, z);
if (j>=0x40900000) { /* z >= 1024 */
if(((j-0x40900000)|i)!=0) /* if z > 1024 */
return s*huge*huge; /* overflow */
Expand All @@ -286,13 +288,13 @@ 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));
SET_HIGH_WORD(t, (n&~(0x000fffff>>k)));
n = ((n&0x000fffff)|0x00100000)>>(20-k);
if(j<0) n = -n;
p_h -= t;
}
t = p_l+p_h;
__LO(t) = 0;
SET_LOW_WORD(t, 0);
u = t*lg2_h;
v = (p_l-(t-p_h))*lg2+t*lg2_l;
z = u+v;
Expand All @@ -301,9 +303,13 @@ 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);
GET_HIGH_WORD(j, z);
j += (n<<20);
if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */
else __HI(z) += (n<<20);
else {
int _hz;
GET_HIGH_WORD(_hz, z);
SET_HIGH_WORD(z, _hz + (n<<20));
}
return s*z;
}
12 changes: 7 additions & 5 deletions external/fdlibm/e_remainder.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,10 +40,10 @@ static double zero = 0.0;
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 */
GET_HIGH_WORD(hx, x); /* high word of x */
GET_LOW_WORD(lx, x); /* low word of x */
GET_HIGH_WORD(hp, p); /* high word of p */
GET_LOW_WORD(lp, p); /* low word of p */
sx = hx&0x80000000;
hp &= 0x7fffffff;
hx &= 0x7fffffff;
Expand Down Expand Up @@ -72,6 +72,8 @@ static double zero = 0.0;
if(x>=p_half) x -= p;
}
}
__HI(x) ^= sx;
int _hx;
GET_HIGH_WORD(_hx, x);
SET_HIGH_WORD(x, _hx ^ sx);
return x;
}
Loading

0 comments on commit b248040

Please sign in to comment.