提交 c6383b7b 编写于 作者: S Szabolcs Nagy

math: use 0x1p-120f and 0x1p120f for tiny and huge values

previously 0x1p-1000 and 0x1p1000 was used for raising inexact
exception like x+tiny (when x is big) or x+huge (when x is small)

the rational is that these float consts are large enough
(0x1p-120 + 1 raises inexact even on ld128 which has 113 mant bits)
and float consts maybe smaller or easier to load on some platforms
(on i386 this reduced the object file size by 4bytes in some cases)
上级 d8a7619e
......@@ -72,7 +72,7 @@ double acos(double x)
if ((ix-0x3ff00000 | lx) == 0) {
/* acos(1)=0, acos(-1)=pi */
if (hx >> 31)
return 2*pio2_hi + 0x1p-1000;
return 2*pio2_hi + 0x1p-120f;
return 0;
}
return 0/(x-x);
......@@ -80,7 +80,7 @@ double acos(double x)
/* |x| < 0.5 */
if (ix < 0x3fe00000) {
if (ix <= 0x3c600000) /* |x| < 2**-57 */
return pio2_hi + 0x1p-1000;
return pio2_hi + 0x1p-120f;
return pio2_hi - (x - (pio2_lo-x*R(x*x)));
}
/* x < -0.5 */
......
......@@ -38,14 +38,14 @@ long double acosl(long double x)
((u.bits.manh & ~LDBL_NBIT) | u.bits.manl) == 0) {
if (expsign > 0)
return 0; /* acos(1) = 0 */
return 2*pio2_hi + 0x1p-1000; /* acos(-1)= pi */
return 2*pio2_hi + 0x1p-120f; /* acos(-1)= pi */
}
return 0/(x-x); /* acos(|x|>1) is NaN */
}
/* |x| < 0.5 */
if (expt < 0x3fff - 1) {
if (expt < 0x3fff - 65)
return pio2_hi + 0x1p-1000; /* x < 0x1p-65: acosl(x)=pi/2 */
return pio2_hi + 0x1p-120f; /* x < 0x1p-65: acosl(x)=pi/2 */
return pio2_hi - (x - (pio2_lo - x * __invtrigl_R(x*x)));
}
/* x < -0.5 */
......
......@@ -77,14 +77,14 @@ double asin(double x)
GET_LOW_WORD(lx, x);
if ((ix-0x3ff00000 | lx) == 0)
/* asin(1) = +-pi/2 with inexact */
return x*pio2_hi + 0x1p-1000;
return x*pio2_hi + 0x1p-120f;
return 0/(x-x);
}
/* |x| < 0.5 */
if (ix < 0x3fe00000) {
if (ix < 0x3e500000) {
/* |x|<0x1p-26, return x with inexact if x!=0*/
FORCE_EVAL(x + 0x1p1000);
FORCE_EVAL(x + 0x1p120f);
return x;
}
return x + x*R(x*x);
......
......@@ -22,7 +22,7 @@ double asinh(double x)
x = log1p(x + x*x/(sqrt(x*x+1)+1));
} else {
/* |x| < 0x1p-26, raise inexact if x != 0 */
FORCE_EVAL(x + 0x1p1000);
FORCE_EVAL(x + 0x1p120f);
}
return s ? -x : x;
}
......@@ -31,7 +31,7 @@ long double asinhl(long double x)
x = log1pl(x + x*x/(sqrtl(x*x+1)+1));
} else {
/* |x| < 0x1p-32, raise inexact if x!=0 */
FORCE_EVAL(x + 0x1p1000);
FORCE_EVAL(x + 0x1p120f);
}
return s ? -x : x;
}
......
......@@ -39,13 +39,13 @@ long double asinl(long double x)
if (expt == 0x3fff &&
((u.bits.manh&~LDBL_NBIT)|u.bits.manl) == 0)
/* asin(+-1)=+-pi/2 with inexact */
return x*pio2_hi + 0x1p-1000;
return x*pio2_hi + 0x1p-120f;
return 0/(x-x);
}
if (expt < 0x3fff - 1) { /* |x| < 0.5 */
if (expt < 0x3fff - 32) { /* |x|<0x1p-32, asinl(x)=x */
/* return x with inexact if x!=0 */
FORCE_EVAL(x + 0x1p1000);
FORCE_EVAL(x + 0x1p120f);
return x;
}
return x + x*__invtrigl_R(x*x);
......
......@@ -72,13 +72,13 @@ double atan(double x)
if (ix >= 0x44100000) { /* if |x| >= 2^66 */
if (isnan(x))
return x;
z = atanhi[3] + 0x1p-1000;
z = atanhi[3] + 0x1p-120f;
return sign ? -z : z;
}
if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
if (ix < 0x3e400000) { /* |x| < 2^-27 */
/* raise inexact if x!=0 */
FORCE_EVAL(x + 0x1p1000);
FORCE_EVAL(x + 0x1p120f);
return x;
}
id = -1;
......
......@@ -52,39 +52,39 @@ long double atan2l(long double y, long double x)
switch(m) {
case 0:
case 1: return y; /* atan(+-0,+anything)=+-0 */
case 2: return 2*pio2_hi+0x1p-1000; /* atan(+0,-anything) = pi */
case 3: return -2*pio2_hi-0x1p-1000; /* atan(-0,-anything) =-pi */
case 2: return 2*pio2_hi+0x1p-120f; /* atan(+0,-anything) = pi */
case 3: return -2*pio2_hi-0x1p-120f; /* atan(-0,-anything) =-pi */
}
}
/* when x = 0 */
if (exptx==0 && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0)
return expsigny < 0 ? -pio2_hi-0x1p-1000 : pio2_hi+0x1p-1000;
return expsigny < 0 ? -pio2_hi-0x1p-120f : pio2_hi+0x1p-120f;
/* when x is INF */
if (exptx == 0x7fff) {
if (expty == 0x7fff) {
switch(m) {
case 0: return pio2_hi*0.5+0x1p-1000; /* atan(+INF,+INF) */
case 1: return -pio2_hi*0.5-0x1p-1000; /* atan(-INF,+INF) */
case 2: return 1.5*pio2_hi+0x1p-1000; /* atan(+INF,-INF) */
case 3: return -1.5*pio2_hi-0x1p-1000; /* atan(-INF,-INF) */
case 0: return pio2_hi*0.5+0x1p-120f; /* atan(+INF,+INF) */
case 1: return -pio2_hi*0.5-0x1p-120f; /* atan(-INF,+INF) */
case 2: return 1.5*pio2_hi+0x1p-120f; /* atan(+INF,-INF) */
case 3: return -1.5*pio2_hi-0x1p-120f; /* atan(-INF,-INF) */
}
} else {
switch(m) {
case 0: return 0.0; /* atan(+...,+INF) */
case 1: return -0.0; /* atan(-...,+INF) */
case 2: return 2*pio2_hi+0x1p-1000; /* atan(+...,-INF) */
case 3: return -2*pio2_hi-0x1p-1000; /* atan(-...,-INF) */
case 2: return 2*pio2_hi+0x1p-120f; /* atan(+...,-INF) */
case 3: return -2*pio2_hi-0x1p-120f; /* atan(-...,-INF) */
}
}
}
/* when y is INF */
if (expty == 0x7fff)
return expsigny < 0 ? -pio2_hi-0x1p-1000 : pio2_hi+0x1p-1000;
return expsigny < 0 ? -pio2_hi-0x1p-120f : pio2_hi+0x1p-120f;
/* compute y/x */
k = expty-exptx;
if(k > LDBL_MANT_DIG+2) { /* |y/x| huge */
z = pio2_hi+0x1p-1000;
z = pio2_hi+0x1p-120f;
m &= 1;
} else if (expsignx < 0 && k < -LDBL_MANT_DIG-2) /* |y/x| tiny, x<0 */
z = 0.0;
......
......@@ -80,7 +80,7 @@ long double atanl(long double x)
if (expt == 0x7fff &&
((u.bits.manh&~LDBL_NBIT)|u.bits.manl)!=0) /* NaN */
return x+x;
z = atanhi[3] + 0x1p-1000;
z = atanhi[3] + 0x1p-120f;
return expsign < 0 ? -z : z;
}
/* Extract the exponent and the first few bits of the mantissa. */
......@@ -89,7 +89,7 @@ long double atanl(long double x)
if (expman < ((0x3fff - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
if (expt < 0x3fff - 32) { /* if |x| is small, atanl(x)~=x */
/* raise inexact if x!=0 */
FORCE_EVAL(x + 0x1p1000);
FORCE_EVAL(x + 0x1p120f);
return x;
}
id = -1;
......
......@@ -97,11 +97,11 @@ float exp2f(float x)
return x;
}
if (x >= 128) {
STRICT_ASSIGN(float, x, x * 0x1p127);
STRICT_ASSIGN(float, x, x * 0x1p127f);
return x;
}
if (x <= -150) {
STRICT_ASSIGN(float, x, 0x1p-100*0x1p-100);
STRICT_ASSIGN(float, x, 0x1p-100f*0x1p-100f);
return x;
}
} else if (ix <= 0x33000000) { /* |x| <= 0x1p-25 */
......
Markdown is supported
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册