Get rid of strict aliasing rule violations from libm (#3069)

JerryScript-DCO-1.0-Signed-off-by: Csaba Osztrogonác oszi@inf.u-szeged.hu
This commit is contained in:
Csaba Osztrogonác
2019-10-01 14:37:18 +02:00
committed by Robert Fancsik
parent f7391a94ae
commit 3763ac8371
16 changed files with 246 additions and 197 deletions
+6 -5
View File
@@ -70,7 +70,7 @@
double double
acos (double x) acos (double x)
{ {
double z, p, q, r, w, s, c, df; double z, p, q, r, w, s, c;
int hx, ix; int hx, ix;
hx = __HI (x); hx = __HI (x);
@@ -114,16 +114,17 @@ acos (double x)
} }
else /* x > 0.5 */ else /* x > 0.5 */
{ {
double_accessor df;
z = (one - x) * 0.5; z = (one - x) * 0.5;
s = sqrt (z); s = sqrt (z);
df = s; df.dbl = s;
__LO (df) = 0; df.as_int.lo = 0;
c = (z - df * df) / (s + df); c = (z - df.dbl * df.dbl) / (s + df.dbl);
p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5))))); p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4))); q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
r = p / q; r = p / q;
w = r * s + c; w = r * s + c;
return 2.0 * (df + w); return 2.0 * (df.dbl + w);
} }
} /* acos */ } /* acos */
+12 -11
View File
@@ -77,7 +77,8 @@
double double
asin (double x) asin (double x)
{ {
double t, w, p, q, c, r, s; double t, p, q, c, r, s;
double_accessor w;
int hx, ix; int hx, ix;
hx = __HI (x); hx = __HI (x);
@@ -102,28 +103,28 @@ asin (double x)
t = x * x; t = x * x;
p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5))))); p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4))); q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4)));
w = p / q; w.dbl = p / q;
return x + x * w; return x + x * w.dbl;
} }
/* 1 > |x| >= 0.5 */ /* 1 > |x| >= 0.5 */
w = one - fabs (x); w.dbl = one - fabs (x);
t = w * 0.5; t = w.dbl * 0.5;
p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5))))); p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4))); q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4)));
s = sqrt (t); s = sqrt (t);
if (ix >= 0x3FEF3333) /* if |x| > 0.975 */ if (ix >= 0x3FEF3333) /* if |x| > 0.975 */
{ {
w = p / q; w.dbl = p / q;
t = pio2_hi - (2.0 * (s + s * w) - pio2_lo); t = pio2_hi - (2.0 * (s + s * w.dbl) - pio2_lo);
} }
else else
{ {
w = s; w.dbl = s;
__LO (w) = 0; w.as_int.lo = 0;
c = (t - w * w) / (s + w); c = (t - w.dbl * w.dbl) / (s + w.dbl);
r = p / q; r = p / q;
p = 2.0 * s * r - (pio2_lo - 2.0 * c); p = 2.0 * s * r - (pio2_lo - 2.0 * c);
q = pio4_hi - 2.0 * w; q = pio4_hi - 2.0 * w.dbl;
t = pio4_hi - (p - q); t = pio4_hi - (p - q);
} }
if (hx > 0) if (hx > 0)
+9 -9
View File
@@ -64,7 +64,7 @@
double double
atan2 (double y, double x) atan2 (double y, double x)
{ {
double z; double_accessor z;
int k, m, hx, hy, ix, iy; int k, m, hx, hy, ix, iy;
unsigned lx, ly; unsigned lx, ly;
@@ -168,35 +168,35 @@ atan2 (double y, double x)
k = (iy - ix) >> 20; k = (iy - ix) >> 20;
if (k > 60) /* |y / x| > 2**60 */ if (k > 60) /* |y / x| > 2**60 */
{ {
z = pi_o_2 + 0.5 * pi_lo; z.dbl = pi_o_2 + 0.5 * pi_lo;
} }
else if (hx < 0 && k < -60) /* |y| / x < -2**60 */ else if (hx < 0 && k < -60) /* |y| / x < -2**60 */
{ {
z = 0.0; z.dbl = 0.0;
} }
else /* safe to do y / x */ else /* safe to do y / x */
{ {
z = atan (fabs (y / x)); z.dbl = atan (fabs (y / x));
} }
switch (m) switch (m)
{ {
case 0: /* atan(+,+) */ case 0: /* atan(+,+) */
{ {
return z; return z.dbl;
} }
case 1: /* atan(-,+) */ case 1: /* atan(-,+) */
{ {
__HI (z) ^= 0x80000000; z.as_int.hi ^= 0x80000000;
return z; return z.dbl;
} }
case 2: /* atan(+,-) */ case 2: /* atan(+,-) */
{ {
return pi - (z - pi_lo); return pi - (z.dbl - pi_lo);
} }
/* case 3: */ /* case 3: */
default: /* atan(-,-) */ default: /* atan(-,-) */
{ {
return (z - pi_lo) - pi; return (z.dbl - pi_lo) - pi;
} }
} }
} /* atan2 */ } /* atan2 */
+5 -3
View File
@@ -123,9 +123,11 @@ ceil (double x)
i1 &= (~i); i1 &= (~i);
} }
} }
__HI (x) = i0;
__LO (x) = i1; double_accessor ret;
return x; ret.as_int.hi = i0;
ret.as_int.lo = i1;
return ret.dbl;
} /* ceil */ } /* ceil */
#undef huge #undef huge
+4 -2
View File
@@ -34,6 +34,8 @@
double double
copysign (double x, double y) copysign (double x, double y)
{ {
__HI (x) = (__HI (x) & 0x7fffffff) | (__HI (y) & 0x80000000); double_accessor ret;
return x; ret.dbl = x;
ret.as_int.hi = (__HI (x) & 0x7fffffff) | (__HI (y) & 0x80000000);
return ret.dbl;
} /* copysign */ } /* copysign */
+8 -6
View File
@@ -120,7 +120,7 @@ static const double ln2LO[2] =
double double
exp (double x) /* default IEEE double exp */ exp (double x) /* default IEEE double exp */
{ {
double y, hi, lo, c, t; double hi, lo, c, t;
int k = 0, xsb; int k = 0, xsb;
unsigned hx; unsigned hx;
@@ -182,6 +182,8 @@ exp (double x) /* default IEEE double exp */
k = 0; k = 0;
} }
double_accessor ret;
/* x is now in primary range */ /* x is now in primary range */
t = x * x; t = x * x;
c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
@@ -191,17 +193,17 @@ exp (double x) /* default IEEE double exp */
} }
else else
{ {
y = one - ((lo - (x * c) / (2.0 - c)) - hi); ret.dbl = one - ((lo - (x * c) / (2.0 - c)) - hi);
} }
if (k >= -1021) if (k >= -1021)
{ {
__HI (y) += (k << 20); /* add k to y's exponent */ ret.as_int.hi += (k << 20); /* add k to y's exponent */
return y; return ret.dbl;
} }
else else
{ {
__HI (y) += ((k + 1000) << 20); /* add k to y's exponent */ ret.as_int.hi += ((k + 1000) << 20); /* add k to y's exponent */
return y * twom1000; return ret.dbl * twom1000;
} }
} /* exp */ } /* exp */
+4 -2
View File
@@ -33,6 +33,8 @@
double double
fabs (double x) fabs (double x)
{ {
__HI (x) &= 0x7fffffff; double_accessor ret;
return x; ret.dbl = x;
ret.as_int.hi &= 0x7fffffff;
return ret.dbl;
} /* fabs */ } /* fabs */
+5 -3
View File
@@ -122,9 +122,11 @@ floor (double x)
i1 &= (~i); i1 &= (~i);
} }
} }
__HI (x) = i0;
__LO (x) = i1; double_accessor ret;
return x; ret.as_int.hi = i0;
ret.as_int.lo = i1;
return ret.dbl;
} /* floor */ } /* floor */
#undef huge #undef huge
+7 -10
View File
@@ -35,8 +35,6 @@
static const double Zero[] = { 0.0, -0.0, }; static const double Zero[] = { 0.0, -0.0, };
#define one 1.0
double double
fmod (double x, double y) fmod (double x, double y)
{ {
@@ -201,11 +199,13 @@ fmod (double x, double y)
lx = lx + lx; lx = lx + lx;
iy -= 1; iy -= 1;
} }
double_accessor ret;
if (iy >= -1022) /* normalize output */ if (iy >= -1022) /* normalize output */
{ {
hx = ((hx - 0x00100000) | ((iy + 1023) << 20)); hx = ((hx - 0x00100000) | ((iy + 1023) << 20));
__HI (x) = hx | sx; ret.as_int.hi = hx | sx;
__LO (x) = lx; ret.as_int.lo = lx;
} }
else /* subnormal output */ else /* subnormal output */
{ {
@@ -225,11 +225,8 @@ fmod (double x, double y)
lx = hx >> (n - 32); lx = hx >> (n - 32);
hx = sx; hx = sx;
} }
__HI (x) = hx | sx; ret.as_int.hi = hx | sx;
__LO (x) = lx; ret.as_int.lo = lx;
x *= one; /* create necessary signal */
} }
return x; /* exact output */ return ret.dbl; /* exact output */
} /* fmod */ } /* fmod */
#undef one
+27 -4
View File
@@ -48,13 +48,36 @@
#endif /* !__LITTLE_ENDIAN */ #endif /* !__LITTLE_ENDIAN */
#ifdef __LITTLE_ENDIAN #ifdef __LITTLE_ENDIAN
#define __HI(x) *(1 + (int *) &x) #define __HI(x) *(1 + (const int *) &x)
#define __LO(x) *(int *) &x #define __LO(x) *(const int *) &x
typedef union
{
double dbl;
struct
{
int lo;
int hi;
} as_int;
} double_accessor;
#else /* !__LITTLE_ENDIAN */ #else /* !__LITTLE_ENDIAN */
#define __HI(x) *(int *) &x #define __HI(x) *(const int *) &x
#define __LO(x) *(1 + (int *) &x) #define __LO(x) *(1 + (const int *) &x)
typedef union
{
double dbl;
struct
{
int hi;
int lo;
} as_int;
} double_accessor;
#endif /* __LITTLE_ENDIAN */ #endif /* __LITTLE_ENDIAN */
#ifndef NAN
#define NAN (0.0/0.0)
#endif
/* /*
* ANSI/POSIX * ANSI/POSIX
*/ */
+6 -2
View File
@@ -122,9 +122,13 @@ log (double 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 */
double_accessor temp;
temp.dbl = x;
temp.as_int.hi = hx | (i ^ 0x3ff00000); /* normalize x or x / 2 */
k += (i >> 20); k += (i >> 20);
f = x - 1.0; f = temp.dbl - 1.0;
if ((0x000fffff & (2 + hx)) < 3) /* |f| < 2**-20 */ if ((0x000fffff & (2 + hx)) < 3) /* |f| < 2**-20 */
{ {
if (f == zero) if (f == zero)
+12 -11
View File
@@ -33,6 +33,7 @@ nextafter (double x,
{ {
int hx, hy, ix, iy; int hx, hy, ix, iy;
unsigned lx, ly; unsigned lx, ly;
double_accessor ret;
hx = __HI (x); /* high word of x */ hx = __HI (x); /* high word of x */
lx = __LO (x); /* low word of x */ lx = __LO (x); /* low word of x */
@@ -54,16 +55,16 @@ nextafter (double x,
if ((ix | lx) == 0) if ((ix | lx) == 0)
{ /* x == 0 */ { /* x == 0 */
__HI (x) = hy & 0x80000000; /* return +-minsubnormal */ ret.as_int.hi = hy & 0x80000000; /* return +-minsubnormal */
__LO (x) = 1; ret.as_int.lo = 1;
y = x * x; y = ret.dbl * ret.dbl;
if (y == x) if (y == ret.dbl)
{ {
return y; return y;
} }
else else
{ {
return x; /* raise underflow flag */ return ret.dbl; /* raise underflow flag */
} }
} }
@@ -122,13 +123,13 @@ nextafter (double x,
y = x * x; y = x * x;
if (y != x) if (y != x)
{ /* raise underflow flag */ { /* raise underflow flag */
__HI (y) = hx; ret.as_int.hi = hx;
__LO (y) = lx; ret.as_int.lo = lx;
return y; return ret.dbl;
} }
} }
__HI (x) = hx; ret.as_int.hi = hx;
__LO (x) = lx; ret.as_int.lo = lx;
return x; return ret.dbl;
} /* nextafter */ } /* nextafter */
+71 -69
View File
@@ -120,8 +120,9 @@ static const double dp_l[] =
double double
pow (double x, double y) pow (double x, double y)
{ {
double z, ax, z_h, z_l, p_h, p_l; double_accessor t1, ax, p_h, y1, t, z;
double y1, t1, t2, r, s, t, u, v, w; double z_h, z_l, p_l;
double t2, r, s, u, v, w;
int i, j, k, yisint, n; int i, j, k, yisint, n;
int hx, hy, ix, iy; int hx, hy, ix, iy;
unsigned lx, ly; unsigned lx, ly;
@@ -227,29 +228,29 @@ pow (double x, double y)
} }
} }
ax = fabs (x); ax.dbl = fabs (x);
/* special value of x */ /* special value of x */
if (lx == 0) if (lx == 0)
{ {
if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000)
{ {
z = ax; /* x is +-0,+-inf,+-1 */ z.dbl = ax.dbl; /* x is +-0,+-inf,+-1 */
if (hy < 0) if (hy < 0)
{ {
z = one / z; /* z = (1 / |x|) */ z.dbl = one / z.dbl; /* z = (1 / |x|) */
} }
if (hx < 0) if (hx < 0)
{ {
if (((ix - 0x3ff00000) | yisint) == 0) if (((ix - 0x3ff00000) | yisint) == 0)
{ {
z = (z - z) / (z - z); /* (-1)**non-int is NaN */ z.dbl = NAN; /* (-1)**non-int is NaN */
} }
else if (yisint == 1) else if (yisint == 1)
{ {
z = -z; /* (x<0)**odd = -(|x|**odd) */ z.dbl = -z.dbl; /* (x<0)**odd = -(|x|**odd) */
} }
} }
return z; return z.dbl;
} }
} }
@@ -292,25 +293,26 @@ pow (double x, double y)
} }
/* now |1 - x| is tiny <= 2**-20, suffice to compute /* now |1 - x| is tiny <= 2**-20, suffice to compute
log(x) by x - x^2 / 2 + x^3 / 3 - x^4 / 4 */ log(x) by x - x^2 / 2 + x^3 / 3 - x^4 / 4 */
t = ax - one; /* t has 20 trailing zeros */ t.dbl = ax.dbl - one; /* t has 20 trailing zeros */
w = (t * t) * (0.5 - t * (0.3333333333333333333333 - t * 0.25)); w = (t.dbl * t.dbl) * (0.5 - t.dbl * (0.3333333333333333333333 - t.dbl * 0.25));
u = ivln2_h * t; /* ivln2_h has 21 sig. bits */ u = ivln2_h * t.dbl; /* ivln2_h has 21 sig. bits */
v = t * ivln2_l - w * ivln2; v = t.dbl * ivln2_l - w * ivln2;
t1 = u + v; t1.dbl = u + v;
__LO (t1) = 0; t1.as_int.lo = 0;
t2 = v - (t1 - u); t2 = v - (t1.dbl - u);
} }
else else
{ {
double ss, s2, s_h, s_l, t_h, t_l; double_accessor s_h, t_h;
double ss, s2, s_l, t_l;
n = 0; n = 0;
/* take care subnormal number */ /* take care subnormal number */
if (ix < 0x00100000) if (ix < 0x00100000)
{ {
ax *= two53; ax.dbl *= two53;
n -= 53; n -= 53;
ix = __HI (ax); ix = ax.as_int.hi;
} }
n += ((ix) >> 20) - 0x3ff; n += ((ix) >> 20) - 0x3ff;
j = ix & 0x000fffff; j = ix & 0x000fffff;
@@ -330,51 +332,51 @@ pow (double x, double y)
n += 1; n += 1;
ix -= 0x00100000; ix -= 0x00100000;
} }
__HI (ax) = ix; ax.as_int.hi = ix;
/* compute ss = s_h + s_l = (x - 1) / (x + 1) or (x - 1.5) / (x + 1.5) */ /* 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 */ u = ax.dbl - bp[k]; /* bp[0] = 1.0, bp[1] = 1.5 */
v = one / (ax + bp[k]); v = one / (ax.dbl + bp[k]);
ss = u * v; ss = u * v;
s_h = ss; s_h.dbl = ss;
__LO (s_h) = 0; s_h.as_int.lo = 0;
/* t_h = ax + bp[k] High */ /* t_h = ax + bp[k] High */
t_h = zero; t_h.dbl = zero;
__HI (t_h) = ((ix >> 1) | 0x20000000) + 0x00080000 + (k << 18); t_h.as_int.hi = ((ix >> 1) | 0x20000000) + 0x00080000 + (k << 18);
t_l = ax - (t_h - bp[k]); t_l = ax.dbl - (t_h.dbl - bp[k]);
s_l = v * ((u - s_h * t_h) - s_h * t_l); s_l = v * ((u - s_h.dbl * t_h.dbl) - s_h.dbl * t_l);
/* compute log(ax) */ /* compute log(ax) */
s2 = ss * ss; s2 = ss * ss;
r = s2 * s2 * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6))))); r = s2 * s2 * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6)))));
r += s_l * (s_h + ss); r += s_l * (s_h.dbl + ss);
s2 = s_h * s_h; s2 = s_h.dbl * s_h.dbl;
t_h = 3.0 + s2 + r; t_h.dbl = 3.0 + s2 + r;
__LO (t_h) = 0; t_h.as_int.lo = 0;
t_l = r - ((t_h - 3.0) - s2); t_l = r - ((t_h.dbl - 3.0) - s2);
/* u + v = ss * (1 + ...) */ /* u + v = ss * (1 + ...) */
u = s_h * t_h; u = s_h.dbl * t_h.dbl;
v = s_l * t_h + t_l * ss; v = s_l * t_h.dbl + t_l * ss;
/* 2 / (3 * log2) * (ss + ...) */ /* 2 / (3 * log2) * (ss + ...) */
p_h = u + v; p_h.dbl = u + v;
__LO (p_h) = 0; p_h.as_int.lo = 0;
p_l = v - (p_h - u); p_l = v - (p_h.dbl - u);
z_h = cp_h * p_h; /* cp_h + cp_l = 2 / (3 * log2) */ z_h = cp_h * p_h.dbl; /* cp_h + cp_l = 2 / (3 * log2) */
z_l = cp_l * p_h + p_l * cp + dp_l[k]; z_l = cp_l * p_h.dbl + p_l * cp + dp_l[k];
/* log2(ax) = (ss + ...) * 2 / (3 * log2) = n + dp_h + z_h + z_l */ /* log2(ax) = (ss + ...) * 2 / (3 * log2) = n + dp_h + z_h + z_l */
t = (double) n; t.dbl = (double) n;
t1 = (((z_h + z_l) + dp_h[k]) + t); t1.dbl = (((z_h + z_l) + dp_h[k]) + t.dbl);
__LO (t1) = 0; t1.as_int.lo = 0;
t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); t2 = z_l - (((t1.dbl - t.dbl) - dp_h[k]) - z_h);
} }
/* split up y into y1 + y2 and compute (y1 + y2) * (t1 + t2) */ /* split up y into y1 + y2 and compute (y1 + y2) * (t1 + t2) */
y1 = y; y1.dbl = y;
__LO (y1) = 0; y1.as_int.lo = 0;
p_l = (y - y1) * t1 + y * t2; p_l = (y - y1.dbl) * t1.dbl + y * t2;
p_h = y1 * t1; p_h.dbl = y1.dbl * t1.dbl;
z = p_l + p_h; z.dbl = p_l + p_h.dbl;
j = __HI (z); j = z.as_int.hi;
i = __LO (z); i = z.as_int.lo;
if (j >= 0x40900000) /* z >= 1024 */ if (j >= 0x40900000) /* z >= 1024 */
{ {
if (((j - 0x40900000) | i) != 0) /* if z > 1024 */ if (((j - 0x40900000) | i) != 0) /* if z > 1024 */
@@ -383,7 +385,7 @@ pow (double x, double y)
} }
else else
{ {
if (p_l + ovt > z - p_h) if (p_l + ovt > z.dbl - p_h.dbl)
{ {
return s * huge * huge; /* overflow */ return s * huge * huge; /* overflow */
} }
@@ -397,7 +399,7 @@ pow (double x, double y)
} }
else else
{ {
if (p_l <= z - p_h) if (p_l <= z.dbl - p_h.dbl)
{ {
return s * tiny * tiny; /* underflow */ return s * tiny * tiny; /* underflow */
} }
@@ -413,36 +415,36 @@ pow (double x, double y)
{ {
n = j + (0x00100000 >> (k + 1)); n = j + (0x00100000 >> (k + 1));
k = ((n & 0x7fffffff) >> 20) - 0x3ff; /* new k for n */ k = ((n & 0x7fffffff) >> 20) - 0x3ff; /* new k for n */
t = zero; t.dbl = zero;
__HI (t) = (n & ~(0x000fffff >> k)); t.as_int.hi = (n & ~(0x000fffff >> k));
n = ((n & 0x000fffff) | 0x00100000) >> (20 - k); n = ((n & 0x000fffff) | 0x00100000) >> (20 - k);
if (j < 0) if (j < 0)
{ {
n = -n; n = -n;
} }
p_h -= t; p_h.dbl -= t.dbl;
} }
t = p_l + p_h; t.dbl = p_l + p_h.dbl;
__LO (t) = 0; t.as_int.lo = 0;
u = t * lg2_h; u = t.dbl * lg2_h;
v = (p_l - (t - p_h)) * lg2 + t * lg2_l; v = (p_l - (t.dbl - p_h.dbl)) * lg2 + t.dbl * lg2_l;
z = u + v; z.dbl = u + v;
w = v - (z - u); w = v - (z.dbl - u);
t = z * z; t.dbl = z.dbl * z.dbl;
t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); t1.dbl = z.dbl - t.dbl * (P1 + t.dbl * (P2 + t.dbl * (P3 + t.dbl * (P4 + t.dbl * P5))));
r = (z * t1) / (t1 - two) - (w + z * w); r = (z.dbl * t1.dbl) / (t1.dbl - two) - (w + z.dbl * w);
z = one - (r - z); z.dbl = one - (r - z.dbl);
j = __HI (z); j = z.as_int.hi;
j += (n << 20); j += (n << 20);
if ((j >> 20) <= 0) /* subnormal output */ if ((j >> 20) <= 0) /* subnormal output */
{ {
z = scalbn (z, n); z.dbl = scalbn (z.dbl, n);
} }
else else
{ {
__HI (z) += (n << 20); z.as_int.hi += (n << 20);
} }
return s * z; return s * z.dbl;
} /* pow */ } /* pow */
#undef zero #undef zero
+8 -4
View File
@@ -70,8 +70,10 @@ scalbn (double x, int n)
} }
if (k > 0) /* normal result */ if (k > 0) /* normal result */
{ {
__HI (x) = (hx & 0x800fffff) | (k << 20); double_accessor ret;
return x; ret.dbl = x;
ret.as_int.hi = (hx & 0x800fffff) | (k << 20);
return ret.dbl;
} }
if (k <= -54) if (k <= -54)
{ {
@@ -85,8 +87,10 @@ scalbn (double x, int n)
} }
} }
k += 54; /* subnormal result */ k += 54; /* subnormal result */
__HI (x) = (hx & 0x800fffff) | (k << 20); double_accessor ret;
return x * twom54; ret.dbl = x;
ret.as_int.hi = (hx & 0x800fffff) | (k << 20);
return ret.dbl * twom54;
} /* scalbn */ } /* scalbn */
#undef two54 #undef two54
+9 -8
View File
@@ -104,7 +104,6 @@
double double
sqrt (double x) sqrt (double x)
{ {
double z;
int sign = (int) 0x80000000; int sign = (int) 0x80000000;
unsigned r, t1, s1, ix1, q1; unsigned r, t1, s1, ix1, q1;
int ix0, s0, q, m, t, i; int ix0, s0, q, m, t, i;
@@ -201,19 +200,21 @@ sqrt (double x)
r >>= 1; r >>= 1;
} }
double_accessor ret;
/* use floating add to find out rounding direction */ /* use floating add to find out rounding direction */
if ((ix0 | ix1) != 0) if ((ix0 | ix1) != 0)
{ {
z = one - tiny; /* trigger inexact flag */ ret.dbl = one - tiny; /* trigger inexact flag */
if (z >= one) if (ret.dbl >= one)
{ {
z = one + tiny; ret.dbl = one + tiny;
if (q1 == (unsigned) 0xffffffff) if (q1 == (unsigned) 0xffffffff)
{ {
q1 = 0; q1 = 0;
q += 1; q += 1;
} }
else if (z > one) else if (ret.dbl > one)
{ {
if (q1 == (unsigned) 0xfffffffe) if (q1 == (unsigned) 0xfffffffe)
{ {
@@ -234,9 +235,9 @@ sqrt (double x)
ix1 |= sign; ix1 |= sign;
} }
ix0 += (m << 20); ix0 += (m << 20);
__HI (z) = ix0; ret.as_int.hi = ix0;
__LO (z) = ix1; ret.as_int.lo = ix1;
return z; return ret.dbl;
} /* sqrt */ } /* sqrt */
#undef one #undef one
+53 -48
View File
@@ -475,7 +475,8 @@ static const int npio2_hw[] =
static int static int
__ieee754_rem_pio2 (double x, double *y) __ieee754_rem_pio2 (double x, double *y)
{ {
double z, w, t, r, fn; double_accessor z;
double w, t, r, fn;
double tx[3]; double tx[3];
int e0, i, j, nx, n, ix, hx; int e0, i, j, nx, n, ix, hx;
@@ -491,33 +492,33 @@ __ieee754_rem_pio2 (double x, double *y)
{ {
if (hx > 0) if (hx > 0)
{ {
z = x - pio2_1; z.dbl = x - pio2_1;
if (ix != 0x3ff921fb) /* 33 + 53 bit pi is good enough */ if (ix != 0x3ff921fb) /* 33 + 53 bit pi is good enough */
{ {
y[0] = z - pio2_1t; y[0] = z.dbl - pio2_1t;
y[1] = (z - y[0]) - pio2_1t; y[1] = (z.dbl - y[0]) - pio2_1t;
} }
else /* near pi/2, use 33 + 33 + 53 bit pi */ else /* near pi/2, use 33 + 33 + 53 bit pi */
{ {
z -= pio2_2; z.dbl -= pio2_2;
y[0] = z - pio2_2t; y[0] = z.dbl - pio2_2t;
y[1] = (z - y[0]) - pio2_2t; y[1] = (z.dbl - y[0]) - pio2_2t;
} }
return 1; return 1;
} }
else /* negative x */ else /* negative x */
{ {
z = x + pio2_1; z.dbl = x + pio2_1;
if (ix != 0x3ff921fb) /* 33 + 53 bit pi is good enough */ if (ix != 0x3ff921fb) /* 33 + 53 bit pi is good enough */
{ {
y[0] = z + pio2_1t; y[0] = z.dbl + pio2_1t;
y[1] = (z - y[0]) + pio2_1t; y[1] = (z.dbl - y[0]) + pio2_1t;
} }
else /* near pi/2, use 33 + 33 + 53 bit pi */ else /* near pi/2, use 33 + 33 + 53 bit pi */
{ {
z += pio2_2; z.dbl += pio2_2;
y[0] = z + pio2_2t; y[0] = z.dbl + pio2_2t;
y[1] = (z - y[0]) + pio2_2t; y[1] = (z.dbl - y[0]) + pio2_2t;
} }
return -1; return -1;
} }
@@ -577,15 +578,15 @@ __ieee754_rem_pio2 (double x, double *y)
return 0; return 0;
} }
/* set z = scalbn(|x|, ilogb(x) - 23) */ /* set z = scalbn(|x|, ilogb(x) - 23) */
__LO (z) = __LO (x); z.as_int.lo = __LO (x);
e0 = (ix >> 20) - 1046; /* e0 = ilogb(z) - 23; */ e0 = (ix >> 20) - 1046; /* e0 = ilogb(z) - 23; */
__HI (z) = ix - (e0 << 20); z.as_int.hi = ix - (e0 << 20);
for (i = 0; i < 2; i++) for (i = 0; i < 2; i++)
{ {
tx[i] = (double) ((int) (z)); tx[i] = (double) ((int) (z.dbl));
z = (z - tx[i]) * two24; z.dbl = (z.dbl - tx[i]) * two24;
} }
tx[2] = z; tx[2] = z.dbl;
nx = 3; nx = 3;
while (tx[nx - 1] == zero) /* skip zero term */ while (tx[nx - 1] == zero) /* skip zero term */
{ {
@@ -708,7 +709,7 @@ __kernel_sin (double x, double y, int iy)
static double static double
__kernel_cos (double x, double y) __kernel_cos (double x, double y)
{ {
double a, hz, z, r, qx; double a, hz, z, r;
int ix; int ix;
ix = __HI (x) & 0x7fffffff; /* ix = |x|'s high word */ ix = __HI (x) & 0x7fffffff; /* ix = |x|'s high word */
@@ -727,17 +728,18 @@ __kernel_cos (double x, double y)
} }
else else
{ {
double_accessor qx;
if (ix > 0x3fe90000) /* x > 0.78125 */ if (ix > 0x3fe90000) /* x > 0.78125 */
{ {
qx = 0.28125; qx.dbl = 0.28125;
} }
else else
{ {
__HI (qx) = ix - 0x00200000; /* x / 4 */ qx.as_int.hi = ix - 0x00200000; /* x / 4 */
__LO (qx) = 0; qx.as_int.lo = 0;
} }
hz = 0.5 * z - qx; hz = 0.5 * z - qx.dbl;
a = one - qx; a = one - qx.dbl;
return a - (hz - (z * r - x * y)); return a - (hz - (z * r - x * y));
} }
} /* __kernel_cos */ } /* __kernel_cos */
@@ -794,7 +796,8 @@ __kernel_cos (double x, double y)
static double static double
__kernel_tan (double x, double y, int iy) __kernel_tan (double x, double y, int iy)
{ {
double z, r, v, w, s; double_accessor z;
double r, v, w, s;
int ix, hx; int ix, hx;
hx = __HI (x); /* high word of x */ hx = __HI (x); /* high word of x */
@@ -815,15 +818,16 @@ __kernel_tan (double x, double y, int iy)
} }
else /* compute -1 / (x + y) carefully */ else /* compute -1 / (x + y) carefully */
{ {
double a, t; double a;
double_accessor t;
z = w = x + y; z.dbl = w = x + y;
__LO (z) = 0; z.as_int.lo = 0;
v = y - (z - x); v = y - (z.dbl - x);
t = a = -one / w; t.dbl = a = -one / w;
__LO (t) = 0; t.as_int.lo = 0;
s = one + t * z; s = one + t.dbl * z.dbl;
return t + a * (s + t * v); return t.dbl + a * (s + t.dbl * v);
} }
} }
} }
@@ -835,22 +839,22 @@ __kernel_tan (double x, double y, int iy)
x = -x; x = -x;
y = -y; y = -y;
} }
z = pio4 - x; z.dbl = pio4 - x;
w = pio4lo - y; w = pio4lo - y;
x = z + w; x = z.dbl + w;
y = 0.0; y = 0.0;
} }
z = x * x; z.dbl = x * x;
w = z * z; w = z.dbl * z.dbl;
/* /*
* Break x^5 * (T[1] + x^2 * T[2] + ...) into * Break x^5 * (T[1] + x^2 * T[2] + ...) into
* x^5 (T[1] + x^4 * T[3] + ... + x^20 * T[11]) + * x^5 (T[1] + x^4 * T[3] + ... + x^20 * T[11]) +
* x^5 (x^2 * (T[2] + x^4 * T[4] + ... + x^22 * [T12])) * x^5 (x^2 * (T[2] + x^4 * T[4] + ... + x^22 * [T12]))
*/ */
r = T1 + w * (T3 + w * (T5 + w * (T7 + w * (T9 + w * T11)))); r = T1 + w * (T3 + w * (T5 + w * (T7 + w * (T9 + w * T11))));
v = z * (T2 + w * (T4 + w * (T6 + w * (T8 + w * (T10 + w * T12))))); v = z.dbl * (T2 + w * (T4 + w * (T6 + w * (T8 + w * (T10 + w * T12)))));
s = z * x; s = z.dbl * x;
r = y + z * (s * (r + v) + y); r = y + z.dbl * (s * (r + v) + y);
r += T0 * s; r += T0 * s;
w = x + r; w = x + r;
if (ix >= 0x3FE59428) if (ix >= 0x3FE59428)
@@ -869,15 +873,16 @@ __kernel_tan (double x, double y, int iy)
* -1.0 / (x + r) here * -1.0 / (x + r) here
*/ */
/* compute -1.0 / (x + r) accurately */ /* compute -1.0 / (x + r) accurately */
double a, t; double a;
double_accessor t;
z = w; z.dbl = w;
__LO (z) = 0; z.as_int.lo = 0;
v = r - (z - x); /* z + v = r + x */ v = r - (z.dbl - x); /* z + v = r + x */
t = a = -1.0 / w; /* a = -1.0 / w */ t.dbl = a = -1.0 / w; /* a = -1.0 / w */
__LO (t) = 0; t.as_int.lo = 0;
s = 1.0 + t * z; s = 1.0 + t.dbl * z.dbl;
return t + a * (s + t * v); return t.dbl + a * (s + t.dbl * v);
} }
} /* __kernel_tan */ } /* __kernel_tan */