diff --git a/jerry-libm/acos.c b/jerry-libm/acos.c index d3b688df3..875d41ef5 100644 --- a/jerry-libm/acos.c +++ b/jerry-libm/acos.c @@ -70,7 +70,7 @@ double acos (double x) { - double z, p, q, r, w, s, c, df; + double z, p, q, r, w, s, c; int hx, ix; hx = __HI (x); @@ -114,16 +114,17 @@ acos (double x) } else /* x > 0.5 */ { + double_accessor df; z = (one - x) * 0.5; s = sqrt (z); - df = s; - __LO (df) = 0; - c = (z - df * df) / (s + df); + df.dbl = s; + df.as_int.lo = 0; + c = (z - df.dbl * df.dbl) / (s + df.dbl); p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5))))); q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4))); r = p / q; w = r * s + c; - return 2.0 * (df + w); + return 2.0 * (df.dbl + w); } } /* acos */ diff --git a/jerry-libm/asin.c b/jerry-libm/asin.c index 3db5e36cb..1af9c3c52 100644 --- a/jerry-libm/asin.c +++ b/jerry-libm/asin.c @@ -77,7 +77,8 @@ double asin (double x) { - double t, w, p, q, c, r, s; + double t, p, q, c, r, s; + double_accessor w; int hx, ix; hx = __HI (x); @@ -102,28 +103,28 @@ asin (double x) t = x * x; p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5))))); q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4))); - w = p / q; - return x + x * w; + w.dbl = p / q; + return x + x * w.dbl; } /* 1 > |x| >= 0.5 */ - w = one - fabs (x); - t = w * 0.5; + w.dbl = one - fabs (x); + t = w.dbl * 0.5; p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5))))); q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4))); s = sqrt (t); if (ix >= 0x3FEF3333) /* if |x| > 0.975 */ { - w = p / q; - t = pio2_hi - (2.0 * (s + s * w) - pio2_lo); + w.dbl = p / q; + t = pio2_hi - (2.0 * (s + s * w.dbl) - pio2_lo); } else { - w = s; - __LO (w) = 0; - c = (t - w * w) / (s + w); + w.dbl = s; + w.as_int.lo = 0; + c = (t - w.dbl * w.dbl) / (s + w.dbl); r = p / q; 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); } if (hx > 0) diff --git a/jerry-libm/atan2.c b/jerry-libm/atan2.c index 47764153b..025e128f0 100644 --- a/jerry-libm/atan2.c +++ b/jerry-libm/atan2.c @@ -64,7 +64,7 @@ double atan2 (double y, double x) { - double z; + double_accessor z; int k, m, hx, hy, ix, iy; unsigned lx, ly; @@ -168,35 +168,35 @@ atan2 (double y, double x) k = (iy - ix) >> 20; 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 */ { - z = 0.0; + z.dbl = 0.0; } else /* safe to do y / x */ { - z = atan (fabs (y / x)); + z.dbl = atan (fabs (y / x)); } switch (m) { case 0: /* atan(+,+) */ { - return z; + return z.dbl; } case 1: /* atan(-,+) */ { - __HI (z) ^= 0x80000000; - return z; + z.as_int.hi ^= 0x80000000; + return z.dbl; } case 2: /* atan(+,-) */ { - return pi - (z - pi_lo); + return pi - (z.dbl - pi_lo); } /* case 3: */ default: /* atan(-,-) */ { - return (z - pi_lo) - pi; + return (z.dbl - pi_lo) - pi; } } } /* atan2 */ diff --git a/jerry-libm/ceil.c b/jerry-libm/ceil.c index fda9b2be2..578cb62c3 100644 --- a/jerry-libm/ceil.c +++ b/jerry-libm/ceil.c @@ -123,9 +123,11 @@ ceil (double x) i1 &= (~i); } } - __HI (x) = i0; - __LO (x) = i1; - return x; + + double_accessor ret; + ret.as_int.hi = i0; + ret.as_int.lo = i1; + return ret.dbl; } /* ceil */ #undef huge diff --git a/jerry-libm/copysign.c b/jerry-libm/copysign.c index 48d35f74d..23e209a49 100644 --- a/jerry-libm/copysign.c +++ b/jerry-libm/copysign.c @@ -34,6 +34,8 @@ double copysign (double x, double y) { - __HI (x) = (__HI (x) & 0x7fffffff) | (__HI (y) & 0x80000000); - return x; + double_accessor ret; + ret.dbl = x; + ret.as_int.hi = (__HI (x) & 0x7fffffff) | (__HI (y) & 0x80000000); + return ret.dbl; } /* copysign */ diff --git a/jerry-libm/exp.c b/jerry-libm/exp.c index f20d7d5d1..8970528fc 100644 --- a/jerry-libm/exp.c +++ b/jerry-libm/exp.c @@ -120,7 +120,7 @@ static const double ln2LO[2] = double exp (double x) /* default IEEE double exp */ { - double y, hi, lo, c, t; + double hi, lo, c, t; int k = 0, xsb; unsigned hx; @@ -182,6 +182,8 @@ exp (double x) /* default IEEE double exp */ k = 0; } + double_accessor ret; + /* x is now in primary range */ t = x * x; c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); @@ -191,17 +193,17 @@ exp (double x) /* default IEEE double exp */ } else { - y = one - ((lo - (x * c) / (2.0 - c)) - hi); + ret.dbl = one - ((lo - (x * c) / (2.0 - c)) - hi); } if (k >= -1021) { - __HI (y) += (k << 20); /* add k to y's exponent */ - return y; + ret.as_int.hi += (k << 20); /* add k to y's exponent */ + return ret.dbl; } else { - __HI (y) += ((k + 1000) << 20); /* add k to y's exponent */ - return y * twom1000; + ret.as_int.hi += ((k + 1000) << 20); /* add k to y's exponent */ + return ret.dbl * twom1000; } } /* exp */ diff --git a/jerry-libm/fabs.c b/jerry-libm/fabs.c index ed208357f..6364d02a5 100644 --- a/jerry-libm/fabs.c +++ b/jerry-libm/fabs.c @@ -33,6 +33,8 @@ double fabs (double x) { - __HI (x) &= 0x7fffffff; - return x; + double_accessor ret; + ret.dbl = x; + ret.as_int.hi &= 0x7fffffff; + return ret.dbl; } /* fabs */ diff --git a/jerry-libm/floor.c b/jerry-libm/floor.c index b279cd0af..ade7921d8 100644 --- a/jerry-libm/floor.c +++ b/jerry-libm/floor.c @@ -122,9 +122,11 @@ floor (double x) i1 &= (~i); } } - __HI (x) = i0; - __LO (x) = i1; - return x; + + double_accessor ret; + ret.as_int.hi = i0; + ret.as_int.lo = i1; + return ret.dbl; } /* floor */ #undef huge diff --git a/jerry-libm/fmod.c b/jerry-libm/fmod.c index bb58bb163..7a74bf3f6 100644 --- a/jerry-libm/fmod.c +++ b/jerry-libm/fmod.c @@ -35,8 +35,6 @@ static const double Zero[] = { 0.0, -0.0, }; -#define one 1.0 - double fmod (double x, double y) { @@ -201,11 +199,13 @@ fmod (double x, double y) lx = lx + lx; iy -= 1; } + + double_accessor ret; if (iy >= -1022) /* normalize output */ { hx = ((hx - 0x00100000) | ((iy + 1023) << 20)); - __HI (x) = hx | sx; - __LO (x) = lx; + ret.as_int.hi = hx | sx; + ret.as_int.lo = lx; } else /* subnormal output */ { @@ -225,11 +225,8 @@ fmod (double x, double y) lx = hx >> (n - 32); hx = sx; } - __HI (x) = hx | sx; - __LO (x) = lx; - x *= one; /* create necessary signal */ + ret.as_int.hi = hx | sx; + ret.as_int.lo = lx; } - return x; /* exact output */ + return ret.dbl; /* exact output */ } /* fmod */ - -#undef one diff --git a/jerry-libm/jerry-libm-internal.h b/jerry-libm/jerry-libm-internal.h index 4c6e68c6f..fb29c7f44 100644 --- a/jerry-libm/jerry-libm-internal.h +++ b/jerry-libm/jerry-libm-internal.h @@ -48,13 +48,36 @@ #endif /* !__LITTLE_ENDIAN */ #ifdef __LITTLE_ENDIAN -#define __HI(x) *(1 + (int *) &x) -#define __LO(x) *(int *) &x +#define __HI(x) *(1 + (const int *) &x) +#define __LO(x) *(const int *) &x +typedef union +{ + double dbl; + struct + { + int lo; + int hi; + } as_int; +} double_accessor; #else /* !__LITTLE_ENDIAN */ -#define __HI(x) *(int *) &x -#define __LO(x) *(1 + (int *) &x) +#define __HI(x) *(const 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 */ +#ifndef NAN +#define NAN (0.0/0.0) +#endif + /* * ANSI/POSIX */ diff --git a/jerry-libm/log.c b/jerry-libm/log.c index 7c9e56944..ed7152ac7 100644 --- a/jerry-libm/log.c +++ b/jerry-libm/log.c @@ -122,9 +122,13 @@ log (double x) k += (hx >> 20) - 1023; hx &= 0x000fffff; 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); - f = x - 1.0; + f = temp.dbl - 1.0; + if ((0x000fffff & (2 + hx)) < 3) /* |f| < 2**-20 */ { if (f == zero) diff --git a/jerry-libm/nextafter.c b/jerry-libm/nextafter.c index df6f4476c..f9473d5a7 100644 --- a/jerry-libm/nextafter.c +++ b/jerry-libm/nextafter.c @@ -33,6 +33,7 @@ nextafter (double x, { int hx, hy, ix, iy; unsigned lx, ly; + double_accessor ret; hx = __HI (x); /* high word of x */ lx = __LO (x); /* low word of x */ @@ -54,16 +55,16 @@ nextafter (double x, if ((ix | lx) == 0) { /* x == 0 */ - __HI (x) = hy & 0x80000000; /* return +-minsubnormal */ - __LO (x) = 1; - y = x * x; - if (y == x) + ret.as_int.hi = hy & 0x80000000; /* return +-minsubnormal */ + ret.as_int.lo = 1; + y = ret.dbl * ret.dbl; + if (y == ret.dbl) { return y; } else { - return x; /* raise underflow flag */ + return ret.dbl; /* raise underflow flag */ } } @@ -122,13 +123,13 @@ nextafter (double x, y = x * x; if (y != x) { /* raise underflow flag */ - __HI (y) = hx; - __LO (y) = lx; - return y; + ret.as_int.hi = hx; + ret.as_int.lo = lx; + return ret.dbl; } } - __HI (x) = hx; - __LO (x) = lx; - return x; + ret.as_int.hi = hx; + ret.as_int.lo = lx; + return ret.dbl; } /* nextafter */ diff --git a/jerry-libm/pow.c b/jerry-libm/pow.c index 18f2bc17a..581876c1d 100644 --- a/jerry-libm/pow.c +++ b/jerry-libm/pow.c @@ -120,8 +120,9 @@ static const double dp_l[] = double pow (double x, double y) { - double z, ax, z_h, z_l, p_h, p_l; - double y1, t1, t2, r, s, t, u, v, w; + double_accessor t1, ax, p_h, y1, t, z; + 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; @@ -227,29 +228,29 @@ pow (double x, double y) } } - ax = fabs (x); + ax.dbl = fabs (x); /* special value of x */ if (lx == 0) { 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) { - z = one / z; /* z = (1 / |x|) */ + z.dbl = one / z.dbl; /* z = (1 / |x|) */ } if (hx < 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) { - 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 log(x) by x - x^2 / 2 + x^3 / 3 - x^4 / 4 */ - t = ax - one; /* t has 20 trailing zeros */ - w = (t * t) * (0.5 - t * (0.3333333333333333333333 - t * 0.25)); - u = ivln2_h * t; /* ivln2_h has 21 sig. bits */ - v = t * ivln2_l - w * ivln2; - t1 = u + v; - __LO (t1) = 0; - t2 = v - (t1 - u); + t.dbl = ax.dbl - one; /* t has 20 trailing zeros */ + w = (t.dbl * t.dbl) * (0.5 - t.dbl * (0.3333333333333333333333 - t.dbl * 0.25)); + u = ivln2_h * t.dbl; /* ivln2_h has 21 sig. bits */ + v = t.dbl * ivln2_l - w * ivln2; + t1.dbl = u + v; + t1.as_int.lo = 0; + t2 = v - (t1.dbl - u); } 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; /* take care subnormal number */ if (ix < 0x00100000) { - ax *= two53; + ax.dbl *= two53; n -= 53; - ix = __HI (ax); + ix = ax.as_int.hi; } n += ((ix) >> 20) - 0x3ff; j = ix & 0x000fffff; @@ -330,51 +332,51 @@ pow (double x, double y) n += 1; 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) */ - u = ax - bp[k]; /* bp[0] = 1.0, bp[1] = 1.5 */ - v = one / (ax + bp[k]); + u = ax.dbl - bp[k]; /* bp[0] = 1.0, bp[1] = 1.5 */ + v = one / (ax.dbl + bp[k]); ss = u * v; - s_h = ss; - __LO (s_h) = 0; + s_h.dbl = ss; + s_h.as_int.lo = 0; /* t_h = ax + bp[k] High */ - t_h = zero; - __HI (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); + t_h.dbl = zero; + t_h.as_int.hi = ((ix >> 1) | 0x20000000) + 0x00080000 + (k << 18); + t_l = ax.dbl - (t_h.dbl - bp[k]); + s_l = v * ((u - s_h.dbl * t_h.dbl) - s_h.dbl * t_l); /* compute log(ax) */ s2 = ss * ss; r = s2 * s2 * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6))))); - r += s_l * (s_h + ss); - s2 = s_h * s_h; - t_h = 3.0 + s2 + r; - __LO (t_h) = 0; - t_l = r - ((t_h - 3.0) - s2); + r += s_l * (s_h.dbl + ss); + s2 = s_h.dbl * s_h.dbl; + t_h.dbl = 3.0 + s2 + r; + t_h.as_int.lo = 0; + t_l = r - ((t_h.dbl - 3.0) - s2); /* u + v = ss * (1 + ...) */ - u = s_h * t_h; - v = s_l * t_h + t_l * ss; + u = s_h.dbl * t_h.dbl; + v = s_l * t_h.dbl + t_l * ss; /* 2 / (3 * log2) * (ss + ...) */ - p_h = u + v; - __LO (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]; + p_h.dbl = u + v; + p_h.as_int.lo = 0; + p_l = v - (p_h.dbl - u); + z_h = cp_h * p_h.dbl; /* cp_h + cp_l = 2 / (3 * log2) */ + 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 */ - t = (double) n; - t1 = (((z_h + z_l) + dp_h[k]) + t); - __LO (t1) = 0; - t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); + t.dbl = (double) n; + t1.dbl = (((z_h + z_l) + dp_h[k]) + t.dbl); + t1.as_int.lo = 0; + t2 = z_l - (((t1.dbl - t.dbl) - dp_h[k]) - z_h); } /* split up y into y1 + y2 and compute (y1 + y2) * (t1 + t2) */ - y1 = y; - __LO (y1) = 0; - p_l = (y - y1) * t1 + y * t2; - p_h = y1 * t1; - z = p_l + p_h; - j = __HI (z); - i = __LO (z); + y1.dbl = y; + y1.as_int.lo = 0; + p_l = (y - y1.dbl) * t1.dbl + y * t2; + p_h.dbl = y1.dbl * t1.dbl; + z.dbl = p_l + p_h.dbl; + j = z.as_int.hi; + i = z.as_int.lo; if (j >= 0x40900000) /* z >= 1024 */ { if (((j - 0x40900000) | i) != 0) /* if z > 1024 */ @@ -383,7 +385,7 @@ pow (double x, double y) } else { - if (p_l + ovt > z - p_h) + if (p_l + ovt > z.dbl - p_h.dbl) { return s * huge * huge; /* overflow */ } @@ -397,7 +399,7 @@ pow (double x, double y) } else { - if (p_l <= z - p_h) + if (p_l <= z.dbl - p_h.dbl) { return s * tiny * tiny; /* underflow */ } @@ -413,36 +415,36 @@ pow (double x, double y) { n = j + (0x00100000 >> (k + 1)); k = ((n & 0x7fffffff) >> 20) - 0x3ff; /* new k for n */ - t = zero; - __HI (t) = (n & ~(0x000fffff >> k)); + t.dbl = zero; + t.as_int.hi = (n & ~(0x000fffff >> k)); n = ((n & 0x000fffff) | 0x00100000) >> (20 - k); if (j < 0) { n = -n; } - p_h -= t; + p_h.dbl -= t.dbl; } - t = p_l + p_h; - __LO (t) = 0; - u = t * lg2_h; - v = (p_l - (t - p_h)) * lg2 + t * lg2_l; - z = u + v; - w = v - (z - u); - t = z * z; - 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); + t.dbl = p_l + p_h.dbl; + t.as_int.lo = 0; + u = t.dbl * lg2_h; + v = (p_l - (t.dbl - p_h.dbl)) * lg2 + t.dbl * lg2_l; + z.dbl = u + v; + w = v - (z.dbl - u); + t.dbl = z.dbl * z.dbl; + t1.dbl = z.dbl - t.dbl * (P1 + t.dbl * (P2 + t.dbl * (P3 + t.dbl * (P4 + t.dbl * P5)))); + r = (z.dbl * t1.dbl) / (t1.dbl - two) - (w + z.dbl * w); + z.dbl = one - (r - z.dbl); + j = z.as_int.hi; j += (n << 20); if ((j >> 20) <= 0) /* subnormal output */ { - z = scalbn (z, n); + z.dbl = scalbn (z.dbl, n); } else { - __HI (z) += (n << 20); + z.as_int.hi += (n << 20); } - return s * z; + return s * z.dbl; } /* pow */ #undef zero diff --git a/jerry-libm/scalbn.c b/jerry-libm/scalbn.c index 2ee4f355d..4c61b5af0 100644 --- a/jerry-libm/scalbn.c +++ b/jerry-libm/scalbn.c @@ -70,8 +70,10 @@ scalbn (double x, int n) } if (k > 0) /* normal result */ { - __HI (x) = (hx & 0x800fffff) | (k << 20); - return x; + double_accessor ret; + ret.dbl = x; + ret.as_int.hi = (hx & 0x800fffff) | (k << 20); + return ret.dbl; } if (k <= -54) { @@ -85,8 +87,10 @@ scalbn (double x, int n) } } k += 54; /* subnormal result */ - __HI (x) = (hx & 0x800fffff) | (k << 20); - return x * twom54; + double_accessor ret; + ret.dbl = x; + ret.as_int.hi = (hx & 0x800fffff) | (k << 20); + return ret.dbl * twom54; } /* scalbn */ #undef two54 diff --git a/jerry-libm/sqrt.c b/jerry-libm/sqrt.c index 43135351b..fd45a6f59 100644 --- a/jerry-libm/sqrt.c +++ b/jerry-libm/sqrt.c @@ -104,7 +104,6 @@ double sqrt (double x) { - double z; int sign = (int) 0x80000000; unsigned r, t1, s1, ix1, q1; int ix0, s0, q, m, t, i; @@ -201,19 +200,21 @@ sqrt (double x) r >>= 1; } + double_accessor ret; + /* use floating add to find out rounding direction */ if ((ix0 | ix1) != 0) { - z = one - tiny; /* trigger inexact flag */ - if (z >= one) + ret.dbl = one - tiny; /* trigger inexact flag */ + if (ret.dbl >= one) { - z = one + tiny; + ret.dbl = one + tiny; if (q1 == (unsigned) 0xffffffff) { q1 = 0; q += 1; } - else if (z > one) + else if (ret.dbl > one) { if (q1 == (unsigned) 0xfffffffe) { @@ -234,9 +235,9 @@ sqrt (double x) ix1 |= sign; } ix0 += (m << 20); - __HI (z) = ix0; - __LO (z) = ix1; - return z; + ret.as_int.hi = ix0; + ret.as_int.lo = ix1; + return ret.dbl; } /* sqrt */ #undef one diff --git a/jerry-libm/trig.c b/jerry-libm/trig.c index dd2440d22..4211571e0 100644 --- a/jerry-libm/trig.c +++ b/jerry-libm/trig.c @@ -475,7 +475,8 @@ static const int npio2_hw[] = static int __ieee754_rem_pio2 (double x, double *y) { - double z, w, t, r, fn; + double_accessor z; + double w, t, r, fn; double tx[3]; int e0, i, j, nx, n, ix, hx; @@ -491,33 +492,33 @@ __ieee754_rem_pio2 (double x, double *y) { if (hx > 0) { - z = x - pio2_1; + z.dbl = x - pio2_1; if (ix != 0x3ff921fb) /* 33 + 53 bit pi is good enough */ { - y[0] = z - pio2_1t; - y[1] = (z - y[0]) - pio2_1t; + y[0] = z.dbl - pio2_1t; + y[1] = (z.dbl - y[0]) - pio2_1t; } else /* near pi/2, use 33 + 33 + 53 bit pi */ { - z -= pio2_2; - y[0] = z - pio2_2t; - y[1] = (z - y[0]) - pio2_2t; + z.dbl -= pio2_2; + y[0] = z.dbl - pio2_2t; + y[1] = (z.dbl - y[0]) - pio2_2t; } return 1; } else /* negative x */ { - z = x + pio2_1; + z.dbl = x + pio2_1; if (ix != 0x3ff921fb) /* 33 + 53 bit pi is good enough */ { - y[0] = z + pio2_1t; - y[1] = (z - y[0]) + pio2_1t; + y[0] = z.dbl + pio2_1t; + y[1] = (z.dbl - y[0]) + pio2_1t; } else /* near pi/2, use 33 + 33 + 53 bit pi */ { - z += pio2_2; - y[0] = z + pio2_2t; - y[1] = (z - y[0]) + pio2_2t; + z.dbl += pio2_2; + y[0] = z.dbl + pio2_2t; + y[1] = (z.dbl - y[0]) + pio2_2t; } return -1; } @@ -577,15 +578,15 @@ __ieee754_rem_pio2 (double x, double *y) return 0; } /* 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; */ - __HI (z) = ix - (e0 << 20); + z.as_int.hi = ix - (e0 << 20); for (i = 0; i < 2; i++) { - tx[i] = (double) ((int) (z)); - z = (z - tx[i]) * two24; + tx[i] = (double) ((int) (z.dbl)); + z.dbl = (z.dbl - tx[i]) * two24; } - tx[2] = z; + tx[2] = z.dbl; nx = 3; while (tx[nx - 1] == zero) /* skip zero term */ { @@ -708,7 +709,7 @@ __kernel_sin (double x, double y, int iy) static double __kernel_cos (double x, double y) { - double a, hz, z, r, qx; + double a, hz, z, r; int ix; ix = __HI (x) & 0x7fffffff; /* ix = |x|'s high word */ @@ -727,17 +728,18 @@ __kernel_cos (double x, double y) } else { + double_accessor qx; if (ix > 0x3fe90000) /* x > 0.78125 */ { - qx = 0.28125; + qx.dbl = 0.28125; } else { - __HI (qx) = ix - 0x00200000; /* x / 4 */ - __LO (qx) = 0; + qx.as_int.hi = ix - 0x00200000; /* x / 4 */ + qx.as_int.lo = 0; } - hz = 0.5 * z - qx; - a = one - qx; + hz = 0.5 * z - qx.dbl; + a = one - qx.dbl; return a - (hz - (z * r - x * y)); } } /* __kernel_cos */ @@ -794,7 +796,8 @@ __kernel_cos (double x, double y) static double __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; 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 */ { - double a, t; + double a; + double_accessor t; - z = w = x + y; - __LO (z) = 0; - v = y - (z - x); - t = a = -one / w; - __LO (t) = 0; - s = one + t * z; - return t + a * (s + t * v); + z.dbl = w = x + y; + z.as_int.lo = 0; + v = y - (z.dbl - x); + t.dbl = a = -one / w; + t.as_int.lo = 0; + s = one + t.dbl * z.dbl; + return t.dbl + a * (s + t.dbl * v); } } } @@ -835,22 +839,22 @@ __kernel_tan (double x, double y, int iy) x = -x; y = -y; } - z = pio4 - x; + z.dbl = pio4 - x; w = pio4lo - y; - x = z + w; + x = z.dbl + w; y = 0.0; } - z = x * x; - w = z * z; + z.dbl = x * x; + w = z.dbl * z.dbl; /* * Break x^5 * (T[1] + x^2 * T[2] + ...) into * 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])) */ 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))))); - s = z * x; - r = y + z * (s * (r + v) + y); + v = z.dbl * (T2 + w * (T4 + w * (T6 + w * (T8 + w * (T10 + w * T12))))); + s = z.dbl * x; + r = y + z.dbl * (s * (r + v) + y); r += T0 * s; w = x + r; if (ix >= 0x3FE59428) @@ -869,15 +873,16 @@ __kernel_tan (double x, double y, int iy) * -1.0 / (x + r) here */ /* compute -1.0 / (x + r) accurately */ - double a, t; + double a; + double_accessor t; - z = w; - __LO (z) = 0; - v = r - (z - x); /* z + v = r + x */ - t = a = -1.0 / w; /* a = -1.0 / w */ - __LO (t) = 0; - s = 1.0 + t * z; - return t + a * (s + t * v); + z.dbl = w; + z.as_int.lo = 0; + v = r - (z.dbl - x); /* z + v = r + x */ + t.dbl = a = -1.0 / w; /* a = -1.0 / w */ + t.as_int.lo = 0; + s = 1.0 + t.dbl * z.dbl; + return t.dbl + a * (s + t.dbl * v); } } /* __kernel_tan */