diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm')
63 files changed, 907 insertions, 833 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c b/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c index abc78a35bd..8a4a5bb7b9 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c @@ -36,8 +36,12 @@ __ieee754_acoshl(long double x) { long double t; int64_t hx; - u_int64_t lx; - GET_LDOUBLE_WORDS64(hx,lx,x); + uint64_t lx; + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); if(hx<0x3ff0000000000000LL) { /* x < 1 */ return (x-x)/(x-x); } else if(hx >=0x41b0000000000000LL) { /* x > 2**28 */ diff --git a/sysdeps/ieee754/ldbl-128ibm/e_acosl.c b/sysdeps/ieee754/ldbl-128ibm/e_acosl.c index 5d2af30346..2cb288238b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_acosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_acosl.c @@ -151,26 +151,27 @@ static const long double long double __ieee754_acosl (long double x) { - long double z, r, w, p, q, s, t, f2; - ieee854_long_double_shape_type u; + long double a, z, r, w, p, q, s, t, f2; - u.value = __builtin_fabsl (x); - if (u.value == 1.0L) + if (__glibc_unlikely (__isnanl (x))) + return x + x; + a = __builtin_fabsl (x); + if (a == 1.0L) { if (x > 0.0L) return 0.0; /* acos(1) = 0 */ else return (2.0 * pio2_hi) + (2.0 * pio2_lo); /* acos(-1)= pi */ } - else if (u.value > 1.0L) + else if (a > 1.0L) { return (x - x) / (x - x); /* acos(|x| > 1) is NaN */ } - if (u.value < 0.5L) + if (a < 0.5L) { - if (u.value < 6.938893903907228e-18L) /* |x| < 2**-57 */ + if (a < 6.938893903907228e-18L) /* |x| < 2**-57 */ return pio2_hi + pio2_lo; - if (u.value < 0.4375L) + if (a < 0.4375L) { /* Arcsine of x. */ z = x * x; @@ -199,7 +200,7 @@ __ieee754_acosl (long double x) return z; } /* .4375 <= |x| < .5 */ - t = u.value - 0.4375L; + t = a - 0.4375L; p = ((((((((((P10 * t + P9) * t + P8) * t @@ -230,9 +231,9 @@ __ieee754_acosl (long double x) r = acosr4375 + r; return r; } - else if (u.value < 0.625L) + else if (a < 0.625L) { - t = u.value - 0.5625L; + t = a - 0.5625L; p = ((((((((((rS10 * t + rS9) * t + rS8) * t @@ -264,7 +265,9 @@ __ieee754_acosl (long double x) } else { /* |x| >= .625 */ - z = (one - u.value) * 0.5; + double shi, slo; + + z = (one - a) * 0.5; s = __ieee754_sqrtl (z); /* Compute an extended precision square root from the Newton iteration s -> 0.5 * (s + z / s). @@ -273,12 +276,11 @@ __ieee754_acosl (long double x) Express s = f1 + f2 where f1 * f1 is exactly representable. w = (z - s^2)/2s = (z - f1^2 - 2 f1 f2 - f2^2)/2s . s + w has extended precision. */ - u.value = s; - u.parts32.w2 = 0; - u.parts32.w3 = 0; - f2 = s - u.value; - w = z - u.value * u.value; - w = w - 2.0 * u.value * f2; + ldbl_unpack (s, &shi, &slo); + a = shi; + f2 = slo; + w = z - a * a; + w = w - 2.0 * a * f2; w = w - f2 * f2; w = w / (2.0 * s); /* Arcsine of s. */ diff --git a/sysdeps/ieee754/ldbl-128ibm/e_asinl.c b/sysdeps/ieee754/ldbl-128ibm/e_asinl.c index b395439495..dece11875b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_asinl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_asinl.c @@ -131,19 +131,20 @@ static const long double long double __ieee754_asinl (long double x) { - long double t, w, p, q, c, r, s; + long double a, t, w, p, q, c, r, s; int flag; - ieee854_long_double_shape_type u; + if (__glibc_unlikely (__isnanl (x))) + return x + x; flag = 0; - u.value = __builtin_fabsl (x); - if (u.value == 1.0L) /* |x|>= 1 */ + a = __builtin_fabsl (x); + if (a == 1.0L) /* |x|>= 1 */ return x * pio2_hi + x * pio2_lo; /* asin(1)=+-pi/2 with inexact */ - else if (u.value >= 1.0L) + else if (a >= 1.0L) return (x - x) / (x - x); /* asin(|x|>1) is NaN */ - else if (u.value < 0.5L) + else if (a < 0.5L) { - if (u.value < 6.938893903907228e-18L) /* |x| < 2**-57 */ + if (a < 6.938893903907228e-18L) /* |x| < 2**-57 */ { if (huge + x > one) return x; /* return x with inexact if x!=0 */ @@ -155,9 +156,9 @@ __ieee754_asinl (long double x) flag = 1; } } - else if (u.value < 0.625L) + else if (a < 0.625L) { - t = u.value - 0.5625; + t = a - 0.5625; p = ((((((((((rS10 * t + rS9) * t + rS8) * t @@ -190,7 +191,7 @@ __ieee754_asinl (long double x) else { /* 1 > |x| >= 0.625 */ - w = one - u.value; + w = one - a; t = w * 0.5; } @@ -223,17 +224,14 @@ __ieee754_asinl (long double x) } s = __ieee754_sqrtl (t); - if (u.value > 0.975L) + if (a > 0.975L) { w = p / q; t = pio2_hi - (2.0 * (s + s * w) - pio2_lo); } else { - u.value = s; - u.parts32.w3 = 0; - u.parts32.w2 = 0; - w = u.value; + w = ldbl_high (s); c = (t - w * w) / (s + w); r = p / q; p = 2.0 * s * r - (pio2_lo - 2.0 * c); diff --git a/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c b/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c index 3e0535561c..b625323df3 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c @@ -56,11 +56,15 @@ __ieee754_atan2l(long double y, long double x) { long double z; int64_t k,m,hx,hy,ix,iy; - u_int64_t lx,ly; + uint64_t lx; + double xhi, xlo, yhi; - GET_LDOUBLE_WORDS64(hx,lx,x); + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); ix = hx&0x7fffffffffffffffLL; - GET_LDOUBLE_WORDS64(hy,ly,y); + yhi = ldbl_high (y); + EXTRACT_WORDS64 (hy, yhi); iy = hy&0x7fffffffffffffffLL; if(((ix)>0x7ff0000000000000LL)|| ((iy)>0x7ff0000000000000LL)) /* x or y is NaN */ @@ -70,7 +74,7 @@ __ieee754_atan2l(long double y, long double x) m = ((hy>>63)&1)|((hx>>62)&2); /* 2*sign(x)+sign(y) */ /* when y = 0 */ - if((iy|(ly&0x7fffffffffffffffLL))==0) { + if(iy==0) { switch(m) { case 0: case 1: return y; /* atan(+-0,+anything)=+-0 */ @@ -79,7 +83,7 @@ __ieee754_atan2l(long double y, long double x) } } /* when x = 0 */ - if((ix|(lx&0x7fffffffffffffff))==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; + if(ix==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; /* when x is INF */ if(ix==0x7ff0000000000000LL) { diff --git a/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c b/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c index f35182f03e..29f2e92072 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c @@ -40,8 +40,10 @@ __ieee754_atanhl(long double x) { long double t; int64_t hx,ix; - u_int64_t lx __attribute__ ((unused)); - GET_LDOUBLE_WORDS64(hx,lx,x); + double xhi; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); ix = hx&0x7fffffffffffffffLL; if (ix >= 0x3ff0000000000000LL) { /* |x|>=1 */ if (ix > 0x3ff0000000000000LL) diff --git a/sysdeps/ieee754/ldbl-128ibm/e_coshl.c b/sysdeps/ieee754/ldbl-128ibm/e_coshl.c index 3e8e1875c6..05683bc02f 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_coshl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_coshl.c @@ -41,9 +41,11 @@ __ieee754_coshl (long double x) { long double t,w; int64_t ix; + double xhi; /* High word of |x|. */ - GET_LDOUBLE_MSW64(ix,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (ix, xhi); ix &= 0x7fffffffffffffffLL; /* x is INF or NaN */ diff --git a/sysdeps/ieee754/ldbl-128ibm/e_exp10l.c b/sysdeps/ieee754/ldbl-128ibm/e_exp10l.c index 1eaf2fe398..49121ca315 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_exp10l.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_exp10l.c @@ -36,9 +36,9 @@ __ieee754_exp10l (long double arg) else if (arg > LDBL_MAX_10_EXP + 1) return LDBL_MAX * LDBL_MAX; - u.d = arg; - arg_high = u.dd[0]; - arg_low = u.dd[1]; + u.ld = arg; + arg_high = u.d[0].d; + arg_low = u.d[1].d; exp_high = arg_high * log10_high; exp_low = arg_high * log10_low + arg_low * M_LN10l; return __ieee754_expl (exp_high) * __ieee754_expl (exp_low); diff --git a/sysdeps/ieee754/ldbl-128ibm/e_expl.c b/sysdeps/ieee754/ldbl-128ibm/e_expl.c index 1b994cd7a9..65ef18532d 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_expl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_expl.c @@ -134,18 +134,17 @@ static const long double C[] = { long double __ieee754_expl (long double x) { + long double result, x22; + union ibm_extended_long_double ex2_u, scale_u; + int unsafe; + /* Check for usual case. */ if (isless (x, himark) && isgreater (x, lomark)) { - int tval1, tval2, unsafe, n_i, exponent2; - long double x22, n, result, xl; - union ibm_extended_long_double ex2_u, scale_u; - fenv_t oldenv; - - feholdexcept (&oldenv); -#ifdef FE_TONEAREST - fesetround (FE_TONEAREST); -#endif + int tval1, tval2, n_i, exponent2; + long double n, xl; + + SET_RESTORE_ROUND (FE_TONEAREST); n = __roundl (x*M_1_LN2); x = x-n*M_LN2_0; @@ -162,50 +161,45 @@ __ieee754_expl (long double x) x = x + xl; /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]). */ - ex2_u.d = __expl_table[T_EXPL_RES1 + tval1] - * __expl_table[T_EXPL_RES2 + tval2]; + ex2_u.ld = (__expl_table[T_EXPL_RES1 + tval1] + * __expl_table[T_EXPL_RES2 + tval2]); n_i = (int)n; /* 'unsafe' is 1 iff n_1 != 0. */ unsafe = fabsl(n_i) >= -LDBL_MIN_EXP - 1; - ex2_u.ieee.exponent += n_i >> unsafe; + ex2_u.d[0].ieee.exponent += n_i >> unsafe; /* Fortunately, there are no subnormal lowpart doubles in __expl_table, only normal values and zeros. But after scaling it can be subnormal. */ - exponent2 = ex2_u.ieee.exponent2 + (n_i >> unsafe); - if (ex2_u.ieee.exponent2 == 0) - /* assert ((ex2_u.ieee.mantissa2|ex2_u.ieee.mantissa3) == 0) */; + exponent2 = ex2_u.d[1].ieee.exponent + (n_i >> unsafe); + if (ex2_u.d[1].ieee.exponent == 0) + /* assert ((ex2_u.d[1].ieee.mantissa0|ex2_u.d[1].ieee.mantissa1) == 0) */; else if (exponent2 > 0) - ex2_u.ieee.exponent2 = exponent2; + ex2_u.d[1].ieee.exponent = exponent2; else if (exponent2 <= -54) { - ex2_u.ieee.exponent2 = 0; - ex2_u.ieee.mantissa2 = 0; - ex2_u.ieee.mantissa3 = 0; + ex2_u.d[1].ieee.exponent = 0; + ex2_u.d[1].ieee.mantissa0 = 0; + ex2_u.d[1].ieee.mantissa1 = 0; } else { static const double two54 = 1.80143985094819840000e+16, /* 4350000000000000 */ twom54 = 5.55111512312578270212e-17; /* 3C90000000000000 */ - ex2_u.dd[1] *= two54; - ex2_u.ieee.exponent2 += n_i >> unsafe; - ex2_u.dd[1] *= twom54; + ex2_u.d[1].d *= two54; + ex2_u.d[1].ieee.exponent += n_i >> unsafe; + ex2_u.d[1].d *= twom54; } /* Compute scale = 2^n_1. */ - scale_u.d = 1.0L; - scale_u.ieee.exponent += n_i - (n_i >> unsafe); + scale_u.ld = 1.0L; + scale_u.d[0].ieee.exponent += n_i - (n_i >> unsafe); /* Approximate e^x2 - 1, using a seventh-degree polynomial, with maximum error in [-2^-16-2^-53,2^-16+2^-53] less than 4.8e-39. */ x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6))))); - /* Return result. */ - fesetenv (&oldenv); - - result = x22 * ex2_u.d + ex2_u.d; - /* Now we can test whether the result is ultimate or if we are unsure. In the later case we should probably call a mpn based routine to give the ultimate result. @@ -235,10 +229,6 @@ __ieee754_expl (long double x) return __ieee754_expl_proc2 (origx); } */ - if (!unsafe) - return result; - else - return result * scale_u.d; } /* Exceptional cases: */ else if (isless (x, himark)) @@ -253,5 +243,10 @@ __ieee754_expl (long double x) else /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ return TWO1023*x; + + result = x22 * ex2_u.ld + ex2_u.ld; + if (!unsafe) + return result; + return result * scale_u.ld; } strong_alias (__ieee754_expl, __expl_finite) diff --git a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c b/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c index a60963c84d..a140fb322d 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c @@ -27,76 +27,83 @@ static const long double one = 1.0, Zero[] = {0.0, -0.0,}; long double __ieee754_fmodl (long double x, long double y) { - int64_t n,hx,hy,hz,ix,iy,sx, i; - u_int64_t lx,ly,lz; - int temp; + int64_t hx, hy, hz, sx, sy; + uint64_t lx, ly, lz; + int n, ix, iy; + double xhi, xlo, yhi, ylo; - GET_LDOUBLE_WORDS64(hx,lx,x); - GET_LDOUBLE_WORDS64(hy,ly,y); + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); + ldbl_unpack (y, &yhi, &ylo); + EXTRACT_WORDS64 (hy, yhi); + EXTRACT_WORDS64 (ly, ylo); sx = hx&0x8000000000000000ULL; /* sign of x */ - hx ^=sx; /* |x| */ - hy &= 0x7fffffffffffffffLL; /* |y| */ + hx ^= sx; /* |x| */ + sy = hy&0x8000000000000000ULL; /* sign of y */ + hy ^= sy; /* |y| */ /* purge off exception values */ - if(__builtin_expect((hy|(ly&0x7fffffffffffffff))==0 || + if(__builtin_expect(hy==0 || (hx>=0x7ff0000000000000LL)|| /* y=0,or x not finite */ (hy>0x7ff0000000000000LL),0)) /* or y is NaN */ return (x*y)/(x*y); - if(__builtin_expect(hx<=hy,0)) { - if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */ - if(lx==ly) - return Zero[(u_int64_t)sx>>63]; /* |x|=|y| return x*0*/ + if (__builtin_expect (hx <= hy, 0)) + { + /* If |x| < |y| return x. */ + if (hx < hy) + return x; + /* At this point the absolute value of the high doubles of + x and y must be equal. */ + /* If the low double of y is the same sign as the high + double of y (ie. the low double increases |y|)... */ + if (((ly ^ sy) & 0x8000000000000000LL) == 0 + /* ... then a different sign low double to high double + for x or same sign but lower magnitude... */ + && (int64_t) (lx ^ sx) < (int64_t) (ly ^ sy)) + /* ... means |x| < |y|. */ + return x; + /* If the low double of x differs in sign to the high + double of x (ie. the low double decreases |x|)... */ + if (((lx ^ sx) & 0x8000000000000000LL) != 0 + /* ... then a different sign low double to high double + for y with lower magnitude (we've already caught + the same sign for y case above)... */ + && (int64_t) (lx ^ sx) > (int64_t) (ly ^ sy)) + /* ... means |x| < |y|. */ + return x; + /* If |x| == |y| return x*0. */ + if ((lx ^ sx) == (ly ^ sy)) + return Zero[(uint64_t) sx >> 63]; } - /* determine ix = ilogb(x) */ - if(__builtin_expect(hx<0x0010000000000000LL,0)) { /* subnormal x */ - if(hx==0) { - for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; - } else { - for (ix = -1022, i=(hx<<11); i>0; i<<=1) ix -=1; - } - } else ix = (hx>>52)-0x3ff; - - /* determine iy = ilogb(y) */ - if(__builtin_expect(hy<0x0010000000000000LL,0)) { /* subnormal y */ - if(hy==0) { - for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; - } else { - for (iy = -1022, i=(hy<<11); i>0; i<<=1) iy -=1; - } - } else iy = (hy>>52)-0x3ff; - /* Make the IBM extended format 105 bit mantissa look like the ieee854 112 bit mantissa so the following operations will give the correct result. */ - ldbl_extract_mantissa(&hx, &lx, &temp, x); - ldbl_extract_mantissa(&hy, &ly, &temp, y); + ldbl_extract_mantissa(&hx, &lx, &ix, x); + ldbl_extract_mantissa(&hy, &ly, &iy, y); - /* set up {hx,lx}, {hy,ly} and align y to x */ - if(__builtin_expect(ix >= -1022, 1)) - hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx); - else { /* subnormal x, shift x to normal */ - n = -1022-ix; - if(n<=63) { - hx = (hx<<n)|(lx>>(64-n)); - lx <<= n; - } else { - hx = lx<<(n-64); - lx = 0; - } - } - if(__builtin_expect(iy >= -1022, 1)) - hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy); - else { /* subnormal y, shift y to normal */ - n = -1022-iy; - if(n<=63) { - hy = (hy<<n)|(ly>>(64-n)); - ly <<= n; - } else { - hy = ly<<(n-64); - ly = 0; - } - } + if (__builtin_expect (ix == -IEEE754_DOUBLE_BIAS, 0)) + { + /* subnormal x, shift x to normal. */ + while ((hx & (1LL << 48)) == 0) + { + hx = (hx << 1) | (lx >> 63); + lx = lx << 1; + ix -= 1; + } + } + + if (__builtin_expect (iy == -IEEE754_DOUBLE_BIAS, 0)) + { + /* subnormal y, shift y to normal. */ + while ((hy & (1LL << 48)) == 0) + { + hy = (hy << 1) | (ly >> 63); + ly = ly << 1; + iy -= 1; + } + } /* fix point fmod */ n = ix - iy; @@ -104,7 +111,7 @@ __ieee754_fmodl (long double x, long double y) hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;} else { - if((hz|(lz&0x7fffffffffffffff))==0) /* return sign(x)*0 */ + if((hz|lz)==0) /* return sign(x)*0 */ return Zero[(u_int64_t)sx>>63]; hx = hz+hz+(lz>>63); lx = lz+lz; } @@ -113,7 +120,7 @@ __ieee754_fmodl (long double x, long double y) if(hz>=0) {hx=hz;lx=lz;} /* convert back to floating value and restore the sign */ - if((hx|(lx&0x7fffffffffffffff))==0) /* return sign(x)*0 */ + if((hx|lx)==0) /* return sign(x)*0 */ return Zero[(u_int64_t)sx>>63]; while(hx<0x0001000000000000LL) { /* normalize x */ hx = hx+hx+(lx>>63); lx = lx+lx; diff --git a/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c b/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c index 90d8e3f0d2..84c13de9b8 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c @@ -122,11 +122,12 @@ long double __ieee754_gammal_r (long double x, int *signgamp) { int64_t hx; - u_int64_t lx; + double xhi; - GET_LDOUBLE_WORDS64 (hx, lx, x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); - if (((hx | lx) & 0x7fffffffffffffffLL) == 0) + if ((hx & 0x7fffffffffffffffLL) == 0) { /* Return value for x == 0 is Inf with divide by zero exception. */ *signgamp = 0; diff --git a/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c b/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c index 768bd3b06c..3b07a47b40 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c @@ -45,76 +45,84 @@ #include <math.h> #include <math_private.h> -static const long double two600 = 0x1.0p+600L; -static const long double two1022 = 0x1.0p+1022L; - long double __ieee754_hypotl(long double x, long double y) { - long double a,b,t1,t2,y1,y2,w,kld; + long double a,b,a1,a2,b1,b2,w,kld; int64_t j,k,ha,hb; + double xhi, yhi, hi, lo; - GET_LDOUBLE_MSW64(ha,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (ha, xhi); + yhi = ldbl_high (y); + EXTRACT_WORDS64 (hb, yhi); ha &= 0x7fffffffffffffffLL; - GET_LDOUBLE_MSW64(hb,y); hb &= 0x7fffffffffffffffLL; if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} a = fabsl(a); /* a <- |a| */ b = fabsl(b); /* b <- |b| */ - if((ha-hb)>0x780000000000000LL) {return a+b;} /* x/y > 2**120 */ + if((ha-hb)>0x0780000000000000LL) {return a+b;} /* x/y > 2**120 */ k=0; kld = 1.0L; if(ha > 0x5f30000000000000LL) { /* a>2**500 */ if(ha >= 0x7ff0000000000000LL) { /* Inf or NaN */ - u_int64_t low; w = a+b; /* for sNaN */ - GET_LDOUBLE_LSW64(low,a); - if(((ha&0xfffffffffffffLL)|(low&0x7fffffffffffffffLL))==0) + if(ha == 0x7ff0000000000000LL) w = a; - GET_LDOUBLE_LSW64(low,b); - if(((hb^0x7ff0000000000000LL)|(low&0x7fffffffffffffffLL))==0) + if(hb == 0x7ff0000000000000LL) w = b; return w; } /* scale a and b by 2**-600 */ - ha -= 0x2580000000000000LL; hb -= 0x2580000000000000LL; k += 600; - a /= two600; - b /= two600; - k += 600; - kld = two600; + a *= 0x1p-600L; + b *= 0x1p-600L; + k = 600; + kld = 0x1p+600L; } - if(hb < 0x23d0000000000000LL) { /* b < 2**-450 */ + else if(hb < 0x23d0000000000000LL) { /* b < 2**-450 */ if(hb <= 0x000fffffffffffffLL) { /* subnormal b or 0 */ - u_int64_t low; - GET_LDOUBLE_LSW64(low,b); - if((hb|(low&0x7fffffffffffffffLL))==0) return a; - t1=two1022; /* t1=2^1022 */ - b *= t1; - a *= t1; - k -= 1022; - kld = kld / two1022; + if(hb==0) return a; + a *= 0x1p+1022L; + b *= 0x1p+1022L; + k = -1022; + kld = 0x1p-1022L; } else { /* scale a and b by 2^600 */ - ha += 0x2580000000000000LL; /* a *= 2^600 */ - hb += 0x2580000000000000LL; /* b *= 2^600 */ - k -= 600; - a *= two600; - b *= two600; - kld = kld / two600; + a *= 0x1p+600L; + b *= 0x1p+600L; + k = -600; + kld = 0x1p-600L; } } /* medium size a and b */ w = a-b; if (w>b) { - SET_LDOUBLE_WORDS64(t1,ha,0); - t2 = a-t1; - w = __ieee754_sqrtl(t1*t1-(b*(-b)-t2*(a+t1))); + ldbl_unpack (a, &hi, &lo); + a1 = hi; + a2 = lo; + /* a*a + b*b + = (a1+a2)*a + b*b + = a1*a + a2*a + b*b + = a1*(a1+a2) + a2*a + b*b + = a1*a1 + a1*a2 + a2*a + b*b + = a1*a1 + a2*(a+a1) + b*b */ + w = __ieee754_sqrtl(a1*a1-(b*(-b)-a2*(a+a1))); } else { a = a+a; - SET_LDOUBLE_WORDS64(y1,hb,0); - y2 = b - y1; - SET_LDOUBLE_WORDS64(t1,ha+0x0010000000000000LL,0); - t2 = a - t1; - w = __ieee754_sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b))); + ldbl_unpack (b, &hi, &lo); + b1 = hi; + b2 = lo; + ldbl_unpack (a, &hi, &lo); + a1 = hi; + a2 = lo; + /* a*a + b*b + = a*a + (a-b)*(a-b) - (a-b)*(a-b) + b*b + = a*a + w*w - (a*a - 2*a*b + b*b) + b*b + = w*w + 2*a*b + = w*w + (a1+a2)*b + = w*w + a1*b + a2*b + = w*w + a1*(b1+b2) + a2*b + = w*w + a1*b1 + a1*b2 + a2*b */ + w = __ieee754_sqrtl(a1*b1-(w*(-w)-(a1*b2+a2*b))); } if(k!=0) return w*kld; diff --git a/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c b/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c index 55f87ed422..aeace7c977 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c @@ -31,26 +31,24 @@ static char rcsid[] = "$NetBSD: $"; int __ieee754_ilogbl(long double x) { - int64_t hx,lx; + int64_t hx; int ix; + double xhi; - GET_LDOUBLE_WORDS64(hx,lx,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); hx &= 0x7fffffffffffffffLL; if(hx <= 0x0010000000000000LL) { - if((hx|(lx&0x7fffffffffffffffLL))==0) + if(hx==0) return FP_ILOGB0; /* ilogbl(0) = FP_ILOGB0 */ else /* subnormal x */ - if(hx==0) { - for (ix = -1043; lx>0; lx<<=1) ix -=1; - } else { - for (ix = -1022, hx<<=11; hx>0; hx<<=1) ix -=1; - } + for (ix = -1022, hx<<=11; hx>0; hx<<=1) ix -=1; return ix; } else if (hx<0x7ff0000000000000LL) return (hx>>52)-0x3ff; else if (FP_ILOGBNAN != INT_MAX) { /* ISO C99 requires ilogbl(+-Inf) == INT_MAX. */ - if (((hx^0x7ff0000000000000LL)|lx) == 0) + if (hx==0x7ff0000000000000LL) return INT_MAX; } return FP_ILOGBNAN; diff --git a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c index 40012e41e1..6761a0d26f 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c @@ -70,26 +70,25 @@ static const long double long double __ieee754_jnl (int n, long double x) { - u_int32_t se; + uint32_t se, lx; int32_t i, ix, sgn; long double a, b, temp, di; long double z, w; - ieee854_long_double_shape_type u; + double xhi; /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) * Thus, J(-n,x) = J(n,-x) */ - u.value = x; - se = u.parts32.w0; + xhi = ldbl_high (x); + EXTRACT_WORDS (se, lx, xhi); ix = se & 0x7fffffff; /* if J(n,NaN) is NaN */ if (ix >= 0x7ff00000) { - if ((u.parts32.w0 & 0xfffff) | u.parts32.w1 - | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) + if (((ix - 0x7ff00000) | lx) != 0) return x + x; } @@ -298,27 +297,26 @@ strong_alias (__ieee754_jnl, __jnl_finite) long double __ieee754_ynl (int n, long double x) { - u_int32_t se; + uint32_t se, lx; int32_t i, ix; int32_t sign; long double a, b, temp; - ieee854_long_double_shape_type u; + double xhi; - u.value = x; - se = u.parts32.w0; + xhi = ldbl_high (x); + EXTRACT_WORDS (se, lx, xhi); ix = se & 0x7fffffff; /* if Y(n,NaN) is NaN */ if (ix >= 0x7ff00000) { - if ((u.parts32.w0 & 0xfffff) | u.parts32.w1 - | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) + if (((ix - 0x7ff00000) | lx) != 0) return x + x; } if (x <= 0.0L) { if (x == 0.0L) - return -HUGE_VALL + x; + return ((n < 0 && (n & 1) != 0) ? 1.0L : -1.0L) / 0.0L; if (se & 0x80000000) return zero / (zero * x); } @@ -377,14 +375,16 @@ __ieee754_ynl (int n, long double x) a = __ieee754_y0l (x); b = __ieee754_y1l (x); /* quit if b is -inf */ - u.value = b; - se = u.parts32.w0 & 0xfff00000; + xhi = ldbl_high (b); + GET_HIGH_WORD (se, xhi); + se &= 0xfff00000; for (i = 1; i < n && se != 0xfff00000; i++) { temp = b; b = ((long double) (i + i) / x) * b - a; - u.value = b; - se = u.parts32.w0 & 0xfff00000; + xhi = ldbl_high (b); + GET_HIGH_WORD (se, xhi); + se &= 0xfff00000; a = temp; } } diff --git a/sysdeps/ieee754/ldbl-128ibm/e_log10l.c b/sysdeps/ieee754/ldbl-128ibm/e_log10l.c index fae774cea8..1a6a4a0fa3 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_log10l.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_log10l.c @@ -182,11 +182,13 @@ __ieee754_log10l (long double x) long double z; long double y; int e; - int64_t hx, lx; + int64_t hx; + double xhi; /* Test for domain */ - GET_LDOUBLE_WORDS64 (hx, lx, x); - if (((hx & 0x7fffffffffffffffLL) | (lx & 0x7fffffffffffffffLL)) == 0) + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); + if ((hx & 0x7fffffffffffffffLL) == 0) return (-1.0L / (x - x)); if (hx < 0) return (x - x) / (x - x); diff --git a/sysdeps/ieee754/ldbl-128ibm/e_log2l.c b/sysdeps/ieee754/ldbl-128ibm/e_log2l.c index f0098f6c73..323ded0c0f 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_log2l.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_log2l.c @@ -177,11 +177,13 @@ __ieee754_log2l (x) long double z; long double y; int e; - int64_t hx, lx; + int64_t hx; + double xhi; /* Test for domain */ - GET_LDOUBLE_WORDS64 (hx, lx, x); - if (((hx & 0x7fffffffffffffffLL) | (lx & 0x7fffffffffffffffLL)) == 0) + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); + if ((hx & 0x7fffffffffffffffLL) == 0) return (-1.0L / (x - x)); if (hx < 0) return (x - x) / (x - x); diff --git a/sysdeps/ieee754/ldbl-128ibm/e_logl.c b/sysdeps/ieee754/ldbl-128ibm/e_logl.c index 15b5edfab3..b7db2b9784 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_logl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_logl.c @@ -188,18 +188,20 @@ static const long double long double __ieee754_logl(long double x) { - long double z, y, w; - ieee854_long_double_shape_type u, t; + long double z, y, w, t; unsigned int m; int k, e; + double xhi; + uint32_t hx, lx; - u.value = x; - m = u.parts32.w0; + xhi = ldbl_high (x); + EXTRACT_WORDS (hx, lx, xhi); + m = hx; /* Check for IEEE special cases. */ k = m & 0x7fffffff; /* log(0) = -infinity. */ - if ((k | u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) == 0) + if ((k | lx) == 0) { return -0.5L / ZERO; } @@ -219,7 +221,7 @@ __ieee754_logl(long double x) { z = x - 1.0L; k = 64; - t.value = 1.0L; + t = 1.0L; e = 0; } else @@ -236,10 +238,8 @@ __ieee754_logl(long double x) k = (m - 0xff000) >> 13; /* t is the argument 0.5 + (k+26)/128 of the nearest item to u in the lookup table. */ - t.parts32.w0 = 0x3ff00000 + (k << 13); - t.parts32.w1 = 0; - t.parts32.w2 = 0; - t.parts32.w3 = 0; + INSERT_WORDS (xhi, 0x3ff00000 + (k << 13), 0); + t = xhi; w0 += 0x100000; e -= 1; k += 64; @@ -247,17 +247,15 @@ __ieee754_logl(long double x) else { k = (m - 0xfe000) >> 14; - t.parts32.w0 = 0x3fe00000 + (k << 14); - t.parts32.w1 = 0; - t.parts32.w2 = 0; - t.parts32.w3 = 0; + INSERT_WORDS (xhi, 0x3fe00000 + (k << 14), 0); + t = xhi; } - u.value = __scalbnl (u.value, ((int) ((w0 - u.parts32.w0) * 2)) >> 21); + x = __scalbnl (x, ((int) ((w0 - hx) * 2)) >> 21); /* log(u) = log( t u/t ) = log(t) + log(u/t) log(t) is tabulated in the lookup table. Express log(u/t) = log(1+z), where z = u/t - 1 = (u-t)/t. cf. Cody & Waite. */ - z = (u.value - t.value) / t.value; + z = (x - t) / t; } /* Series expansion of log(1+z). */ w = z * z; @@ -284,7 +282,7 @@ __ieee754_logl(long double x) y += e * ln2b; /* Base 2 exponent offset times ln(2). */ y += z; y += logtbl[k-26]; /* log(t) - (t-1) */ - y += (t.value - 1.0L); + y += (t - 1.0L); y += e * ln2a; return y; } diff --git a/sysdeps/ieee754/ldbl-128ibm/e_powl.c b/sysdeps/ieee754/ldbl-128ibm/e_powl.c index 8bd35d0c88..c942f2f249 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_powl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_powl.c @@ -151,37 +151,32 @@ __ieee754_powl (long double x, long double y) long double y1, t1, t2, r, s, t, u, v, w; long double s2, s_h, s_l, t_h, t_l, ay; int32_t i, j, k, yisint, n; - u_int32_t ix, iy; - int32_t hx, hy; - ieee854_long_double_shape_type o, p, q; + uint32_t ix, iy; + int32_t hx, hy, hax; + double ohi, xhi, xlo, yhi, ylo; + uint32_t lx, ly, lj; - p.value = x; - hx = p.parts32.w0; + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS (hx, lx, xhi); ix = hx & 0x7fffffff; - q.value = y; - hy = q.parts32.w0; + ldbl_unpack (y, &yhi, &ylo); + EXTRACT_WORDS (hy, ly, yhi); iy = hy & 0x7fffffff; - /* y==zero: x**0 = 1 */ - if ((iy | q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) == 0) + if ((iy | ly) == 0) return one; /* 1.0**y = 1; -1.0**+-Inf = 1 */ if (x == one) return one; - if (x == -1.0L && iy == 0x7ff00000 - && (q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) == 0) + if (x == -1.0L && ((iy - 0x7ff00000) | ly) == 0) return one; /* +-NaN return x+y */ - if ((ix > 0x7ff00000) - || ((ix == 0x7ff00000) - && ((p.parts32.w1 | (p.parts32.w2 & 0x7fffffff) | p.parts32.w3) != 0)) - || (iy > 0x7ff00000) - || ((iy == 0x7ff00000) - && ((q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) != 0))) + if ((ix >= 0x7ff00000 && ((ix - 0x7ff00000) | lx) != 0) + || (iy >= 0x7ff00000 && ((iy - 0x7ff00000) | ly) != 0)) return x + y; /* determine if y is an odd int when x < 0 @@ -192,7 +187,10 @@ __ieee754_powl (long double x, long double y) yisint = 0; if (hx < 0) { - if ((q.parts32.w2 & 0x7fffffff) >= 0x43400000) /* Low part >= 2^53 */ + uint32_t low_ye; + + GET_HIGH_WORD (low_ye, ylo); + if ((low_ye & 0x7fffffff) >= 0x43400000) /* Low part >= 2^53 */ yisint = 2; /* even integer y */ else if (iy >= 0x3ff00000) /* 1.0 */ { @@ -207,42 +205,43 @@ __ieee754_powl (long double x, long double y) } } + ax = fabsl (x); + /* special value of y */ - if ((q.parts32.w1 | (q.parts32.w2 & 0x7fffffff) | q.parts32.w3) == 0) + if (ly == 0) { - if (iy == 0x7ff00000 && q.parts32.w1 == 0) /* y is +-inf */ + if (iy == 0x7ff00000) /* y is +-inf */ { - if (((ix - 0x3ff00000) | p.parts32.w1 - | (p.parts32.w2 & 0x7fffffff) | p.parts32.w3) == 0) - return y - y; /* inf**+-1 is NaN */ - else if (ix > 0x3ff00000 || fabsl (x) > 1.0L) + if (ax > one) /* (|x|>1)**+-inf = inf,0 */ return (hy >= 0) ? y : zero; else /* (|x|<1)**-,+inf = inf,0 */ return (hy < 0) ? -y : zero; } - if (iy == 0x3ff00000) - { /* y is +-1 */ - if (hy < 0) - return one / x; - else - return x; - } - if (hy == 0x40000000) - return x * x; /* y is 2 */ - if (hy == 0x3fe00000) - { /* y is 0.5 */ - if (hx >= 0) /* x >= +0 */ - return __ieee754_sqrtl (x); + if (ylo == 0.0) + { + if (iy == 0x3ff00000) + { /* y is +-1 */ + if (hy < 0) + return one / x; + else + return x; + } + if (hy == 0x40000000) + return x * x; /* y is 2 */ + if (hy == 0x3fe00000) + { /* y is 0.5 */ + if (hx >= 0) /* x >= +0 */ + return __ieee754_sqrtl (x); + } } } - ax = fabsl (x); /* special value of x */ - if ((p.parts32.w1 | (p.parts32.w2 & 0x7fffffff) | p.parts32.w3) == 0) + if (lx == 0) { - if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) + if (ix == 0x7ff00000 || ix == 0 || (ix == 0x3ff00000 && xlo == 0.0)) { z = ax; /*x is +-0,+-inf,+-1 */ if (hy < 0) @@ -294,8 +293,8 @@ __ieee754_powl (long double x, long double y) { ax *= two113; n -= 113; - o.value = ax; - ix = o.parts32.w0; + ohi = ldbl_high (ax); + GET_HIGH_WORD (ix, ohi); } n += ((ix) >> 20) - 0x3ff; j = ix & 0x000fffff; @@ -312,26 +311,19 @@ __ieee754_powl (long double x, long double y) ix -= 0x00100000; } - o.value = ax; - o.value = __scalbnl (o.value, ((int) ((ix - o.parts32.w0) * 2)) >> 21); - ax = o.value; + ohi = ldbl_high (ax); + GET_HIGH_WORD (hax, ohi); + ax = __scalbnl (ax, ((int) ((ix - hax) * 2)) >> 21); /* 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; + s_h = ldbl_high (s); - o.value = s_h; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - s_h = o.value; /* t_h=ax+bp[k] High */ t_h = ax + bp[k]; - o.value = t_h; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - t_h = o.value; + t_h = ldbl_high (t_h); t_l = ax - (t_h - bp[k]); s_l = v * ((u - s_h * t_h) - s_h * t_l); /* compute log(ax) */ @@ -342,30 +334,21 @@ __ieee754_powl (long double x, long double y) r += s_l * (s_h + s); s2 = s_h * s_h; t_h = 3.0 + s2 + r; - o.value = t_h; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - t_h = o.value; + t_h = ldbl_high (t_h); 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; - o.value = p_h; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - p_h = o.value; + p_h = ldbl_high (p_h); 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 = (long double) n; t1 = (((z_h + z_l) + dp_h[k]) + t); - o.value = t1; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - t1 = o.value; + t1 = ldbl_high (t1); t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); /* s (sign of result -ve**odd) = -1 else = 1 */ @@ -374,21 +357,16 @@ __ieee754_powl (long double x, long double y) s = -one; /* (-ve)**(odd int) */ /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ - y1 = y; - o.value = y1; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - y1 = o.value; + y1 = ldbl_high (y); p_l = (y - y1) * t1 + y * t2; p_h = y1 * t1; z = p_l + p_h; - o.value = z; - j = o.parts32.w0; + ohi = ldbl_high (z); + EXTRACT_WORDS (j, lj, ohi); if (j >= 0x40d00000) /* z >= 16384 */ { /* if z > 16384 */ - if (((j - 0x40d00000) | o.parts32.w1 - | (o.parts32.w2 & 0x7fffffff) | o.parts32.w3) != 0) + if (((j - 0x40d00000) | lj) != 0) return s * huge * huge; /* overflow */ else { @@ -399,8 +377,7 @@ __ieee754_powl (long double x, long double y) else if ((j & 0x7fffffff) >= 0x40d01b90) /* z <= -16495 */ { /* z < -16495 */ - if (((j - 0xc0d01bc0) | o.parts32.w1 - | (o.parts32.w2 & 0x7fffffff) | o.parts32.w3) != 0) + if (((j - 0xc0d01bc0) | lj) != 0) return s * tiny * tiny; /* underflow */ else { @@ -419,10 +396,7 @@ __ieee754_powl (long double x, long double y) p_h -= t; } t = p_l + p_h; - o.value = t; - o.parts32.w3 = 0; - o.parts32.w2 = 0; - t = o.value; + t = ldbl_high (t); u = t * lg2_h; v = (p_l - (t - p_h)) * lg2 + t * lg2_l; z = u + v; diff --git a/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c b/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c index 6a72d6a853..36bc03226b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c @@ -200,10 +200,11 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y) double tx[8]; int exp; int64_t n, ix, hx, ixd; - u_int64_t lx __attribute__ ((unused)); u_int64_t lxd; + double xhi; - GET_LDOUBLE_WORDS64 (hx, lx, x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); ix = hx & 0x7fffffffffffffffLL; if (ix <= 0x3fe921fb54442d10LL) /* x in <-pi/4, pi/4> */ { @@ -243,7 +244,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y) We split the 113 bits of the mantissa into 5 24bit integers stored in a double array. */ /* Make the IBM extended format 105 bit mantissa look like the ieee854 112 - bit mantissa so the next operatation will give the correct result. */ + bit mantissa so the next operation will give the correct result. */ ldbl_extract_mantissa (&ixd, &lxd, &exp, x); exp = exp - 23; /* This is faster than doing this in floating point, because we diff --git a/sysdeps/ieee754/ldbl-128ibm/e_remainderl.c b/sysdeps/ieee754/ldbl-128ibm/e_remainderl.c index 67d7db7fbe..800416f29a 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_remainderl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_remainderl.c @@ -33,18 +33,22 @@ __ieee754_remainderl(long double x, long double p) int64_t hx,hp; u_int64_t sx,lx,lp; long double p_half; + double xhi, xlo, phi, plo; - GET_LDOUBLE_WORDS64(hx,lx,x); - GET_LDOUBLE_WORDS64(hp,lp,p); + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); + ldbl_unpack (p, &phi, &plo); + EXTRACT_WORDS64 (hp, phi); + EXTRACT_WORDS64 (lp, plo); sx = hx&0x8000000000000000ULL; hp &= 0x7fffffffffffffffLL; hx &= 0x7fffffffffffffffLL; /* purge off exception values */ - if((hp|(lp&0x7fffffffffffffff))==0) return (x*p)/(x*p); /* p = 0 */ + if(hp==0) return (x*p)/(x*p); /* p = 0 */ if((hx>=0x7ff0000000000000LL)|| /* x not finite */ - ((hp>=0x7ff0000000000000LL)&& /* p is NaN */ - (((hp-0x7ff0000000000000LL)|lp)!=0))) + (hp>0x7ff0000000000000LL)) /* p is NaN */ return (x*p)/(x*p); @@ -64,8 +68,8 @@ __ieee754_remainderl(long double x, long double p) if(x>=p_half) x -= p; } } - GET_LDOUBLE_MSW64(hx,x); - SET_LDOUBLE_MSW64(x,hx^sx); + if (sx) + x = -x; return x; } strong_alias (__ieee754_remainderl, __remainderl_finite) diff --git a/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c b/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c index 4e8481c419..1790bef87e 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c @@ -38,9 +38,11 @@ __ieee754_sinhl(long double x) { long double t,w,h; int64_t ix,jx; + double xhi; /* High word of |x|. */ - GET_LDOUBLE_MSW64(jx,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (jx, xhi); ix = jx&0x7fffffffffffffffLL; /* x is INF or NaN */ diff --git a/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c b/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c index 2b0f7c62e4..61feb367f7 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c @@ -34,15 +34,13 @@ #include <math_private.h> -typedef unsigned int int4; -typedef union {int4 i[4]; long double x; double d[2]; } mynumber; +typedef union {int64_t i[2]; long double x; double d[2]; } mynumber; -static const mynumber - t512 = {{0x5ff00000, 0x00000000, 0x00000000, 0x00000000 }}, /* 2^512 */ - tm256 = {{0x2ff00000, 0x00000000, 0x00000000, 0x00000000 }}; /* 2^-256 */ static const double -two54 = 1.80143985094819840000e+16, /* 0x4350000000000000 */ -twom54 = 5.55111512312578270212e-17; /* 0x3C90000000000000 */ + t512 = 0x1p512, + tm256 = 0x1p-256, + two54 = 0x1p54, /* 0x4350000000000000 */ + twom54 = 0x1p-54; /* 0x3C90000000000000 */ /*********************************************************************/ /* An ultimate sqrt routine. Given an IEEE double machine number x */ @@ -54,56 +52,53 @@ long double __ieee754_sqrtl(long double x) static const long double big = 134217728.0, big1 = 134217729.0; long double t,s,i; mynumber a,c; - int4 k, l, m; - int n; + uint64_t k, l; + int64_t m, n; double d; a.x=x; - k=a.i[0] & 0x7fffffff; + k=a.i[0] & INT64_C(0x7fffffffffffffff); /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/ - if (k>0x000fffff && k<0x7ff00000) { + if (k>INT64_C(0x000fffff00000000) && k<INT64_C(0x7ff0000000000000)) { if (x < 0) return (big1-big1)/(big-big); - l = (k&0x001fffff)|0x3fe00000; - if (((a.i[2] & 0x7fffffff) | a.i[3]) != 0) { - n = (int) ((l - k) * 2) >> 21; - m = (a.i[2] >> 20) & 0x7ff; + l = (k&INT64_C(0x001fffffffffffff))|INT64_C(0x3fe0000000000000); + if ((a.i[1] & INT64_C(0x7fffffffffffffff)) != 0) { + n = (int64_t) ((l - k) * 2) >> 53; + m = (a.i[1] >> 52) & 0x7ff; if (m == 0) { a.d[1] *= two54; - m = ((a.i[2] >> 20) & 0x7ff) - 54; + m = ((a.i[1] >> 52) & 0x7ff) - 54; } m += n; - if ((int) m > 0) - a.i[2] = (a.i[2] & 0x800fffff) | (m << 20); - else if ((int) m <= -54) { - a.i[2] &= 0x80000000; - a.i[3] = 0; + if (m > 0) + a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52); + else if (m <= -54) { + a.i[1] &= INT64_C(0x8000000000000000); } else { m += 54; - a.i[2] = (a.i[2] & 0x800fffff) | (m << 20); + a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52); a.d[1] *= twom54; } } a.i[0] = l; s = a.x; d = __ieee754_sqrt (a.d[0]); - c.i[0] = 0x20000000+((k&0x7fe00000)>>1); + c.i[0] = INT64_C(0x2000000000000000)+((k&INT64_C(0x7fe0000000000000))>>1); c.i[1] = 0; - c.i[2] = 0; - c.i[3] = 0; i = d; t = 0.5L * (i + s / i); i = 0.5L * (t + s / t); return c.x * i; } else { - if (k>=0x7ff00000) { - if (a.i[0] == 0xfff00000 && a.i[1] == 0) + if (k>=INT64_C(0x7ff0000000000000)) { + if (a.i[0] == INT64_C(0xfff0000000000000)) return (big1-big1)/(big-big); /* sqrt (-Inf) = NaN. */ return x; /* sqrt (NaN) = NaN, sqrt (+Inf) = +Inf. */ } if (x == 0) return x; if (x < 0) return (big1-big1)/(big-big); - return tm256.x*__ieee754_sqrtl(x*t512.x); + return tm256*__ieee754_sqrtl(x*t512); } } strong_alias (__ieee754_sqrtl, __sqrtl_finite) diff --git a/sysdeps/ieee754/ldbl-128ibm/ieee754.h b/sysdeps/ieee754/ldbl-128ibm/ieee754.h index 9e94f53b04..0c97a99207 100644 --- a/sysdeps/ieee754/ldbl-128ibm/ieee754.h +++ b/sysdeps/ieee754/ldbl-128ibm/ieee754.h @@ -111,61 +111,6 @@ union ieee754_double #define IEEE754_DOUBLE_BIAS 0x3ff /* Added to exponent. */ -union ieee854_long_double - { - long double d; - - /* This is the IEEE 854 quad-precision format. */ - struct - { -#if __BYTE_ORDER == __BIG_ENDIAN - unsigned int negative:1; - unsigned int exponent:15; - /* Together these comprise the mantissa. */ - unsigned int mantissa0:16; - unsigned int mantissa1:32; - unsigned int mantissa2:32; - unsigned int mantissa3:32; -#endif /* Big endian. */ -#if __BYTE_ORDER == __LITTLE_ENDIAN - /* Together these comprise the mantissa. */ - unsigned int mantissa3:32; - unsigned int mantissa2:32; - unsigned int mantissa1:32; - unsigned int mantissa0:16; - unsigned int exponent:15; - unsigned int negative:1; -#endif /* Little endian. */ - } ieee; - - /* This format makes it easier to see if a NaN is a signalling NaN. */ - struct - { -#if __BYTE_ORDER == __BIG_ENDIAN - unsigned int negative:1; - unsigned int exponent:15; - unsigned int quiet_nan:1; - /* Together these comprise the mantissa. */ - unsigned int mantissa0:15; - unsigned int mantissa1:32; - unsigned int mantissa2:32; - unsigned int mantissa3:32; -#else - /* Together these comprise the mantissa. */ - unsigned int mantissa3:32; - unsigned int mantissa2:32; - unsigned int mantissa1:32; - unsigned int mantissa0:15; - unsigned int quiet_nan:1; - unsigned int exponent:15; - unsigned int negative:1; -#endif - } ieee_nan; - }; - -#define IEEE854_LONG_DOUBLE_BIAS 0x3fff /* Added to exponent. */ - - /* IBM extended format for long double. Each long double is made up of two IEEE doubles. The value of the @@ -179,49 +124,10 @@ union ieee854_long_double union ibm_extended_long_double { - long double d; - double dd[2]; - - /* This is the IBM extended format long double. */ - struct - { /* Big endian. There is no other. */ - - unsigned int negative:1; - unsigned int exponent:11; - /* Together Mantissa0-3 comprise the mantissa. */ - unsigned int mantissa0:20; - unsigned int mantissa1:32; - - unsigned int negative2:1; - unsigned int exponent2:11; - /* There is an implied 1 here? */ - /* Together these comprise the mantissa. */ - unsigned int mantissa2:20; - unsigned int mantissa3:32; - } ieee; - - /* This format makes it easier to see if a NaN is a signalling NaN. */ - struct - { /* Big endian. There is no other. */ - - unsigned int negative:1; - unsigned int exponent:11; - unsigned int quiet_nan:1; - /* Together Mantissa0-3 comprise the mantissa. */ - unsigned int mantissa0:19; - unsigned int mantissa1:32; - - unsigned int negative2:1; - unsigned int exponent2:11; - /* There is an implied 1 here? */ - /* Together these comprise the mantissa. */ - unsigned int mantissa2:20; - unsigned int mantissa3:32; - } ieee_nan; + long double ld; + union ieee754_double d[2]; }; -#define IBM_EXTENDED_LONG_DOUBLE_BIAS 0x3ff /* Added to exponent. */ - __END_DECLS #endif /* ieee754.h */ diff --git a/sysdeps/ieee754/ldbl-128ibm/k_cosl.c b/sysdeps/ieee754/ldbl-128ibm/k_cosl.c index 0b81782fdb..046f3b573c 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_cosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/k_cosl.c @@ -81,8 +81,11 @@ __kernel_cosl(long double x, long double y) { long double h, l, z, sin_l, cos_l_m1; int64_t ix; - u_int32_t tix, hix, index; - GET_LDOUBLE_MSW64 (ix, x); + uint32_t tix, hix, index; + double xhi, hhi; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (ix, xhi); tix = ((u_int64_t)ix) >> 32; tix &= ~0x80000000; /* tix = |x|'s high 32 bits */ if (tix < 0x3fc30000) /* |x| < 0.1484375 */ @@ -136,7 +139,8 @@ __kernel_cosl(long double x, long double y) case 2: index = (hix - 0x3fc30000) >> 14; break; } */ - SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0); + INSERT_WORDS64 (hhi, ((uint64_t)hix) << 32); + h = hhi; l = y - (h - x); z = l * l; sin_l = l*(ONE+z*(SSIN1+z*(SSIN2+z*(SSIN3+z*(SSIN4+z*SSIN5))))); diff --git a/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c b/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c index fc1ead6597..3ba9d7e907 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c @@ -100,9 +100,12 @@ __kernel_sincosl(long double x, long double y, long double *sinx, long double *c { long double h, l, z, sin_l, cos_l_m1; int64_t ix; - u_int32_t tix, hix, index; - GET_LDOUBLE_MSW64 (ix, x); - tix = ((u_int64_t)ix) >> 32; + uint32_t tix, hix, index; + double xhi, hhi; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (ix, xhi); + tix = ((uint64_t)ix) >> 32; tix &= ~0x80000000; /* tix = |x|'s high 32 bits */ if (tix < 0x3fc30000) /* |x| < 0.1484375 */ { @@ -164,7 +167,8 @@ __kernel_sincosl(long double x, long double y, long double *sinx, long double *c case 2: index = (hix - 0x3fc30000) >> 14; break; } */ - SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0); + INSERT_WORDS64 (hhi, ((uint64_t)hix) << 32); + h = hhi; if (iy) l = y - (h - x); else diff --git a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c index f17c0ae5db..b12ea134d5 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c +++ b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c @@ -82,7 +82,10 @@ __kernel_sinl(long double x, long double y, int iy) long double h, l, z, sin_l, cos_l_m1; int64_t ix; u_int32_t tix, hix, index; - GET_LDOUBLE_MSW64 (ix, x); + double xhi, hhi; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (ix, xhi); tix = ((u_int64_t)ix) >> 32; tix &= ~0x80000000; /* tix = |x|'s high 32 bits */ if (tix < 0x3fc30000) /* |x| < 0.1484375 */ @@ -132,7 +135,8 @@ __kernel_sinl(long double x, long double y, int iy) case 2: index = (hix - 0x3fc30000) >> 14; break; } */ - SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0); + INSERT_WORDS64 (hhi, ((uint64_t)hix) << 32); + h = hhi; if (iy) l = (ix < 0 ? -y : y) - (h - x); else diff --git a/sysdeps/ieee754/ldbl-128ibm/k_tanl.c b/sysdeps/ieee754/ldbl-128ibm/k_tanl.c index 1f6bad241b..bcf8b5e7d6 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_tanl.c +++ b/sysdeps/ieee754/ldbl-128ibm/k_tanl.c @@ -85,17 +85,17 @@ long double __kernel_tanl (long double x, long double y, int iy) { long double z, r, v, w, s; - int32_t ix, sign; - ieee854_long_double_shape_type u, u1; + int32_t ix, sign, hx, lx; + double xhi; - u.value = x; - ix = u.parts32.w0 & 0x7fffffff; + xhi = ldbl_high (x); + EXTRACT_WORDS (hx, lx, xhi); + ix = hx & 0x7fffffff; if (ix < 0x3c600000) /* x < 2**-57 */ { - if ((int) x == 0) - { /* generate inexact */ - if ((ix | u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3 - | (iy + 1)) == 0) + if ((int) x == 0) /* generate inexact */ + { + if ((ix | lx | (iy + 1)) == 0) return one / fabs (x); else return (iy == 1) ? x : -one / x; @@ -103,7 +103,7 @@ __kernel_tanl (long double x, long double y, int iy) } if (ix >= 0x3fe59420) /* |x| >= 0.6743316650390625 */ { - if ((u.parts32.w0 & 0x80000000) != 0) + if ((hx & 0x80000000) != 0) { x = -x; y = -y; @@ -139,15 +139,13 @@ __kernel_tanl (long double x, long double y, int iy) { /* if allow error up to 2 ulp, simply return -1.0/(x+r) here */ /* compute -1.0/(x+r) accurately */ - u1.value = w; - u1.parts32.w2 = 0; - u1.parts32.w3 = 0; - v = r - (u1.value - x); /* u1+v = r+x */ + long double u1, z1; + + u1 = ldbl_high (w); + v = r - (u1 - x); /* u1+v = r+x */ z = -1.0 / w; - u.value = z; - u.parts32.w2 = 0; - u.parts32.w3 = 0; - s = 1.0 + u.value * u1.value; - return u.value + z * (s + u.value * v); + z1 = ldbl_high (z); + s = 1.0 + z1 * u1; + return z1 + z * (s + z1 * v); } } diff --git a/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c b/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c index 00e44b8b92..e46fde74fc 100644 --- a/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c +++ b/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c @@ -36,34 +36,44 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, union ibm_extended_long_double u; unsigned long long hi, lo; int ediff; - u.d = value; - *is_neg = u.ieee.negative; - *expt = (int) u.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS; + u.ld = value; - lo = ((long long) u.ieee.mantissa2 << 32) | u.ieee.mantissa3; - hi = ((long long) u.ieee.mantissa0 << 32) | u.ieee.mantissa1; - /* If the lower double is not a denomal or zero then set the hidden + *is_neg = u.d[0].ieee.negative; + *expt = (int) u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS; + + lo = ((long long) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1; + hi = ((long long) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1; + + /* If the lower double is not a denormal or zero then set the hidden 53rd bit. */ - if (u.ieee.exponent2 > 0) - { - lo |= 1LL << 52; + if (u.d[1].ieee.exponent != 0) + lo |= 1ULL << 52; + else + lo = lo << 1; - /* The lower double is normalized separately from the upper. We may - need to adjust the lower manitissa to reflect this. */ - ediff = u.ieee.exponent - u.ieee.exponent2; - if (ediff > 53) - lo = lo >> (ediff-53); + /* The lower double is normalized separately from the upper. We may + need to adjust the lower manitissa to reflect this. */ + ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53; + if (ediff > 0) + { + if (ediff < 64) + lo = lo >> ediff; + else + lo = 0; } + else if (ediff < 0) + lo = lo << -ediff; + /* The high double may be rounded and the low double reflects the difference between the long double and the rounded high double value. This is indicated by a differnce between the signs of the high and low doubles. */ - if ((u.ieee.negative != u.ieee.negative2) - && ((u.ieee.exponent2 != 0) && (lo != 0L))) + if (u.d[0].ieee.negative != u.d[1].ieee.negative + && lo != 0) { lo = (1ULL << 53) - lo; - if (hi == 0LL) + if (hi == 0) { /* we have a borrow from the hidden bit, so shift left 1. */ hi = 0x0ffffffffffffeLL | (lo >> 51); @@ -92,7 +102,7 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, #define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \ - (LDBL_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB))) - if (u.ieee.exponent == 0) + if (u.d[0].ieee.exponent == 0) { /* A biased exponent of zero is a special case. Either it is a zero or it is a denormal number. */ diff --git a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h b/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h index 046293e598..1b6e27a9ff 100644 --- a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h +++ b/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h @@ -2,10 +2,13 @@ #error "Never use <math_ldbl.h> directly; include <math_private.h> instead." #endif -#include <sysdeps/ieee754/ldbl-128/math_ldbl.h> #include <ieee754.h> #include <stdint.h> +/* To suit our callers we return *hi64 and *lo64 as if they came from + an ieee854 112 bit mantissa, that is, 48 bits in *hi64 (plus one + implicit bit) and 64 bits in *lo64. */ + static inline void ldbl_extract_mantissa (int64_t *hi64, uint64_t *lo64, int *exp, long double x) { @@ -14,77 +17,119 @@ ldbl_extract_mantissa (int64_t *hi64, uint64_t *lo64, int *exp, long double x) the number before the decimal point and the second implicit bit as bit 53 of the mantissa. */ uint64_t hi, lo; - int ediff; - union ibm_extended_long_double eldbl; - eldbl.d = x; - *exp = eldbl.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS; - - lo = ((int64_t)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3; - hi = ((int64_t)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1; - /* If the lower double is not a denomal or zero then set the hidden - 53rd bit. */ - if (eldbl.ieee.exponent2 > 0x001) - { - lo |= (1ULL << 52); - lo = lo << 7; /* pre-shift lo to match ieee854. */ - /* The lower double is normalized separately from the upper. We - may need to adjust the lower manitissa to reflect this. */ - ediff = eldbl.ieee.exponent - eldbl.ieee.exponent2; - if (ediff > 53) - lo = lo >> (ediff-53); - hi |= (1ULL << 52); - } + union ibm_extended_long_double u; - if ((eldbl.ieee.negative != eldbl.ieee.negative2) - && ((eldbl.ieee.exponent2 != 0) && (lo != 0LL))) + u.ld = x; + *exp = u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS; + + lo = ((uint64_t) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1; + hi = ((uint64_t) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1; + + if (u.d[0].ieee.exponent != 0) { - hi--; - lo = (1ULL << 60) - lo; - if (hi < (1ULL << 52)) + int ediff; + + /* If not a denormal or zero then we have an implicit 53rd bit. */ + hi |= (uint64_t) 1 << 52; + + if (u.d[1].ieee.exponent != 0) + lo |= (uint64_t) 1 << 52; + else + /* A denormal is to be interpreted as having a biased exponent + of 1. */ + lo = lo << 1; + + /* We are going to shift 4 bits out of hi later, because we only + want 48 bits in *hi64. That means we want 60 bits in lo, but + we currently only have 53. Shift the value up. */ + lo = lo << 7; + + /* The lower double is normalized separately from the upper. + We may need to adjust the lower mantissa to reflect this. + The difference between the exponents can be larger than 53 + when the low double is much less than 1ULP of the upper + (in which case there are significant bits, all 0's or all + 1's, between the two significands). The difference between + the exponents can be less than 53 when the upper double + exponent is nearing its minimum value (in which case the low + double is denormal ie. has an exponent of zero). */ + ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53; + if (ediff > 0) + { + if (ediff < 64) + lo = lo >> ediff; + else + lo = 0; + } + else if (ediff < 0) + lo = lo << -ediff; + + if (u.d[0].ieee.negative != u.d[1].ieee.negative + && lo != 0) { - /* we have a borrow from the hidden bit, so shift left 1. */ - hi = (hi << 1) | (lo >> 59); - lo = 0xfffffffffffffffLL & (lo << 1); - *exp = *exp - 1; + hi--; + lo = ((uint64_t) 1 << 60) - lo; + if (hi < (uint64_t) 1 << 52) + { + /* We have a borrow from the hidden bit, so shift left 1. */ + hi = (hi << 1) | (lo >> 59); + lo = (((uint64_t) 1 << 60) - 1) & (lo << 1); + *exp = *exp - 1; + } } } + else + /* If the larger magnitude double is denormal then the smaller + one must be zero. */ + hi = hi << 1; + *lo64 = (hi << 60) | lo; *hi64 = hi >> 4; } static inline long double -ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64) +ldbl_insert_mantissa (int sign, int exp, int64_t hi64, uint64_t lo64) { union ibm_extended_long_double u; - unsigned long hidden2, lzcount; - unsigned long long hi, lo; + int expnt2; + uint64_t hi, lo; + + u.d[0].ieee.negative = sign; + u.d[1].ieee.negative = sign; + u.d[0].ieee.exponent = exp + IEEE754_DOUBLE_BIAS; + u.d[1].ieee.exponent = 0; + expnt2 = exp - 53 + IEEE754_DOUBLE_BIAS; - u.ieee.negative = sign; - u.ieee.negative2 = sign; - u.ieee.exponent = exp + IBM_EXTENDED_LONG_DOUBLE_BIAS; - u.ieee.exponent2 = exp-53 + IBM_EXTENDED_LONG_DOUBLE_BIAS; /* Expect 113 bits (112 bits + hidden) right justified in two longs. The low order 53 bits (52 + hidden) go into the lower double */ - lo = (lo64 >> 7)& ((1ULL << 53) - 1); - hidden2 = (lo64 >> 59) & 1ULL; + lo = (lo64 >> 7) & (((uint64_t) 1 << 53) - 1); /* The high order 53 bits (52 + hidden) go into the upper double */ - hi = (lo64 >> 60) & ((1ULL << 11) - 1); - hi |= (hi64 << 4); + hi = lo64 >> 60; + hi |= hi64 << 4; - if (lo != 0LL) + if (lo != 0) { - /* hidden2 bit of low double controls rounding of the high double. - If hidden2 is '1' then round up hi and adjust lo (2nd mantissa) + int lzcount; + + /* hidden bit of low double controls rounding of the high double. + If hidden is '1' and either the explicit mantissa is non-zero + or hi is odd, then round up hi and adjust lo (2nd mantissa) plus change the sign of the low double to compensate. */ - if (hidden2) + if ((lo & ((uint64_t) 1 << 52)) != 0 + && ((hi & 1) != 0 || (lo & (((uint64_t) 1 << 52) - 1)) != 0)) { hi++; - u.ieee.negative2 = !sign; - lo = (1ULL << 53) - lo; + if ((hi & ((uint64_t) 1 << 53)) != 0) + { + hi = hi >> 1; + u.d[0].ieee.exponent++; + } + u.d[1].ieee.negative = !sign; + lo = ((uint64_t) 1 << 53) - lo; } - /* The hidden bit of the lo mantissa is zero so we need to - normalize the it for the low double. Shift it left until the - hidden bit is '1' then adjust the 2nd exponent accordingly. */ + + /* Normalize the low double. Shift the mantissa left until + the hidden bit is '1' and adjust the exponent accordingly. */ if (sizeof (lo) == sizeof (long)) lzcount = __builtin_clzl (lo); @@ -92,35 +137,31 @@ ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64) lzcount = __builtin_clzl ((long) (lo >> 32)); else lzcount = __builtin_clzl ((long) lo) + 32; - lzcount = lzcount - 11; - if (lzcount > 0) + lzcount = lzcount - (64 - 53); + lo <<= lzcount; + expnt2 -= lzcount; + + if (expnt2 >= 1) + /* Not denormal. */ + u.d[1].ieee.exponent = expnt2; + else { - int expnt2 = u.ieee.exponent2 - lzcount; - if (expnt2 >= 1) - { - /* Not denormal. Normalize and set low exponent. */ - lo = lo << lzcount; - u.ieee.exponent2 = expnt2; - } + /* Is denormal. Note that biased exponent of 0 is treated + as if it was 1, hence the extra shift. */ + if (expnt2 > -53) + lo >>= 1 - expnt2; else - { - /* Is denormal. */ - lo = lo << (lzcount + expnt2); - u.ieee.exponent2 = 0; - } + lo = 0; } } else - { - u.ieee.negative2 = 0; - u.ieee.exponent2 = 0; - } + u.d[1].ieee.negative = 0; - u.ieee.mantissa3 = lo & ((1ULL << 32) - 1); - u.ieee.mantissa2 = (lo >> 32) & ((1ULL << 20) - 1); - u.ieee.mantissa1 = hi & ((1ULL << 32) - 1); - u.ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 1); - return u.d; + u.d[1].ieee.mantissa1 = lo; + u.d[1].ieee.mantissa0 = lo >> 32; + u.d[0].ieee.mantissa1 = hi; + u.d[0].ieee.mantissa0 = hi >> 32; + return u.ld; } /* Handy utility functions to pack/unpack/cononicalize and find the nearbyint @@ -129,18 +170,18 @@ static inline long double default_ldbl_pack (double a, double aa) { union ibm_extended_long_double u; - u.dd[0] = a; - u.dd[1] = aa; - return u.d; + u.d[0].d = a; + u.d[1].d = aa; + return u.ld; } static inline void default_ldbl_unpack (long double l, double *a, double *aa) { union ibm_extended_long_double u; - u.d = l; - *a = u.dd[0]; - *aa = u.dd[1]; + u.ld = l; + *a = u.d[0].d; + *aa = u.d[1].d; } #ifndef ldbl_pack @@ -150,6 +191,9 @@ default_ldbl_unpack (long double l, double *a, double *aa) # define ldbl_unpack default_ldbl_unpack #endif +/* Extract high double. */ +#define ldbl_high(x) ((double) x) + /* Convert a finite long double to canonical form. Does not handle +/-Inf properly. */ static inline void @@ -163,13 +207,13 @@ ldbl_canonicalize (double *a, double *aa) *aa = xl; } -/* Simple inline nearbyint (double) function . +/* Simple inline nearbyint (double) function. Only works in the default rounding mode but is useful in long double rounding functions. */ static inline double ldbl_nearbyint (double a) { - double two52 = 0x10000000000000LL; + double two52 = 0x1p52; if (__builtin_expect ((__builtin_fabs (a) < two52), 1)) { diff --git a/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c b/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c index 3df42c5666..c96852dfd9 100644 --- a/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c +++ b/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c @@ -33,11 +33,11 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign) unsigned long long hi, lo; int exponent2; - u.ieee.negative = sign; - u.ieee.negative2 = sign; - u.ieee.exponent = expt + IBM_EXTENDED_LONG_DOUBLE_BIAS; - u.ieee.exponent2 = 0; - exponent2 = expt - 53 + IBM_EXTENDED_LONG_DOUBLE_BIAS; + u.d[0].ieee.negative = sign; + u.d[1].ieee.negative = sign; + u.d[0].ieee.exponent = expt + IEEE754_DOUBLE_BIAS; + u.d[1].ieee.exponent = 0; + exponent2 = expt - 53 + IEEE754_DOUBLE_BIAS; #if BITS_PER_MP_LIMB == 32 /* The low order 53 bits (52 + hidden) go into the lower double */ @@ -69,19 +69,19 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign) else lzcount = __builtin_clzl ((long) val) + 32; if (hi) - lzcount = lzcount - 11; + lzcount = lzcount - (64 - 53); else - lzcount = lzcount + 42; + lzcount = lzcount + 53 - (64 - 53); - if (lzcount > u.ieee.exponent) + if (lzcount > u.d[0].ieee.exponent) { - lzcount = u.ieee.exponent; - u.ieee.exponent = 0; + lzcount = u.d[0].ieee.exponent; + u.d[0].ieee.exponent = 0; exponent2 -= lzcount; } else { - u.ieee.exponent -= (lzcount - 1); + u.d[0].ieee.exponent -= (lzcount - 1); exponent2 -= (lzcount - 1); } @@ -97,29 +97,27 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign) } } - if (lo != 0L) + if (lo != 0) { - /* hidden2 bit of low double controls rounding of the high double. - If hidden2 is '1' and either the explicit mantissa is non-zero + /* hidden bit of low double controls rounding of the high double. + If hidden is '1' and either the explicit mantissa is non-zero or hi is odd, then round up hi and adjust lo (2nd mantissa) plus change the sign of the low double to compensate. */ if ((lo & (1LL << 52)) != 0 - && ((hi & 1) != 0 || (lo & ((1LL << 52) - 1)))) + && ((hi & 1) != 0 || (lo & ((1LL << 52) - 1)) != 0)) { hi++; - if ((hi & ((1LL << 52) - 1)) == 0) + if ((hi & (1LL << 53)) != 0) { - if ((hi & (1LL << 53)) != 0) - hi -= 1LL << 52; - u.ieee.exponent++; + hi >>= 1; + u.d[0].ieee.exponent++; } - u.ieee.negative2 = !sign; + u.d[1].ieee.negative = !sign; lo = (1LL << 53) - lo; } - /* The hidden bit of the lo mantissa is zero so we need to normalize - it for the low double. Shift it left until the hidden bit is '1' - then adjust the 2nd exponent accordingly. */ + /* Normalize the low double. Shift the mantissa left until + the hidden bit is '1' and adjust the exponent accordingly. */ if (sizeof (lo) == sizeof (long)) lzcount = __builtin_clzl (lo); @@ -127,24 +125,24 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign) lzcount = __builtin_clzl ((long) (lo >> 32)); else lzcount = __builtin_clzl ((long) lo) + 32; - lzcount = lzcount - 11; - if (lzcount > 0) - { - lo = lo << lzcount; - exponent2 = exponent2 - lzcount; - } + lzcount = lzcount - (64 - 53); + lo <<= lzcount; + exponent2 -= lzcount; + if (exponent2 > 0) - u.ieee.exponent2 = exponent2; - else + u.d[1].ieee.exponent = exponent2; + else if (exponent2 > -53) lo >>= 1 - exponent2; + else + lo = 0; } else - u.ieee.negative2 = 0; + u.d[1].ieee.negative = 0; - u.ieee.mantissa3 = lo & 0xffffffffLL; - u.ieee.mantissa2 = (lo >> 32) & 0xfffff; - u.ieee.mantissa1 = hi & 0xffffffffLL; - u.ieee.mantissa0 = (hi >> 32) & ((1LL << (LDBL_MANT_DIG - 86)) - 1); + u.d[1].ieee.mantissa1 = lo; + u.d[1].ieee.mantissa0 = lo >> 32; + u.d[0].ieee.mantissa1 = hi; + u.d[0].ieee.mantissa0 = hi >> 32; - return u.d; + return u.ld; } diff --git a/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c b/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c index 247dc20e5a..e0ec422b01 100644 --- a/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c +++ b/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c @@ -26,31 +26,31 @@ do { \ unsigned long long int num0, num1; \ unsigned long long hi, lo; \ int ediff; \ - union ibm_extended_long_double eldbl; \ - eldbl.d = fpnum.ldbl.d; \ + union ibm_extended_long_double u; \ + u.ld = fpnum.ldbl; \ \ assert (sizeof (long double) == 16); \ \ - lo = ((long long)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3; \ - hi = ((long long)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1; \ + lo = ((long long)u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1; \ + hi = ((long long)u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1; \ lo <<= 7; /* pre-shift lo to match ieee854. */ \ - /* If the lower double is not a denomal or zero then set the hidden \ + /* If the lower double is not a denormal or zero then set the hidden \ 53rd bit. */ \ - if (eldbl.ieee.exponent2 != 0) \ + if (u.d[1].ieee.exponent != 0) \ lo |= (1ULL << (52 + 7)); \ else \ lo <<= 1; \ /* The lower double is normalized separately from the upper. We \ may need to adjust the lower manitissa to reflect this. */ \ - ediff = eldbl.ieee.exponent - eldbl.ieee.exponent2; \ - if (ediff > 53 + 63) \ + ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53; \ + if (ediff > 63) \ lo = 0; \ - else if (ediff > 53) \ - lo = lo >> (ediff - 53); \ - else if (eldbl.ieee.exponent2 == 0 && ediff < 53) \ - lo = lo << (53 - ediff); \ - if (eldbl.ieee.negative != eldbl.ieee.negative2 \ - && (eldbl.ieee.exponent2 != 0 || lo != 0L)) \ + else if (ediff > 0) \ + lo = lo >> ediff; \ + else if (ediff < 0) \ + lo = lo << -ediff; \ + if (u.d[0].ieee.negative != u.d[1].ieee.negative \ + && lo != 0) \ { \ lo = (1ULL << 60) - lo; \ if (hi == 0L) \ @@ -58,7 +58,7 @@ do { \ /* we have a borrow from the hidden bit, so shift left 1. */ \ hi = 0xffffffffffffeLL | (lo >> 59); \ lo = 0xfffffffffffffffLL & (lo << 1); \ - eldbl.ieee.exponent--; \ + u.d[0].ieee.exponent--; \ } \ else \ hi--; \ @@ -109,9 +109,9 @@ do { \ *--wnumstr = L'0'; \ } \ \ - leading = eldbl.ieee.exponent == 0 ? '0' : '1'; \ + leading = u.d[0].ieee.exponent == 0 ? '0' : '1'; \ \ - exponent = eldbl.ieee.exponent; \ + exponent = u.d[0].ieee.exponent; \ \ if (exponent == 0) \ { \ @@ -121,18 +121,18 @@ do { \ { \ /* This is a denormalized number. */ \ expnegative = 1; \ - exponent = IBM_EXTENDED_LONG_DOUBLE_BIAS - 1; \ + exponent = IEEE754_DOUBLE_BIAS - 1; \ } \ } \ - else if (exponent >= IBM_EXTENDED_LONG_DOUBLE_BIAS) \ + else if (exponent >= IEEE754_DOUBLE_BIAS) \ { \ expnegative = 0; \ - exponent -= IBM_EXTENDED_LONG_DOUBLE_BIAS; \ + exponent -= IEEE754_DOUBLE_BIAS; \ } \ else \ { \ expnegative = 1; \ - exponent = -(exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS); \ + exponent = -(exponent - IEEE754_DOUBLE_BIAS); \ } \ } while (0) diff --git a/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c b/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c index a833457eab..63c6edbb1a 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_asinhl.c @@ -38,7 +38,10 @@ long double __asinhl(long double x) { long double t,w; int64_t hx,ix; - GET_LDOUBLE_MSW64(hx,x); + double xhi; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); ix = hx&0x7fffffffffffffffLL; if(ix>=0x7ff0000000000000LL) return x+x; /* x is inf or NaN */ if(ix< 0x3e20000000000000LL) { /* |x|<2**-29 */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_atanl.c b/sysdeps/ieee754/ldbl-128ibm/s_atanl.c index 2a36d16bb4..41dde23998 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_atanl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_atanl.c @@ -173,23 +173,20 @@ static const long double long double __atanl (long double x) { - int k, sign; + int32_t k, sign, lx; long double t, u, p, q; - ieee854_long_double_shape_type s; + double xhi; - s.value = x; - k = s.parts32.w0; - if (k & 0x80000000) - sign = 1; - else - sign = 0; + xhi = ldbl_high (x); + EXTRACT_WORDS (k, lx, xhi); + sign = k & 0x80000000; /* Check for IEEE special cases. */ k &= 0x7fffffff; if (k >= 0x7ff00000) { /* NaN. */ - if ((k & 0xfffff) | s.parts32.w1 ) + if (((k - 0x7ff00000) | lx) != 0) return (x + x); /* Infinity. */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_cosl.c b/sysdeps/ieee754/ldbl-128ibm/s_cosl.c index 23148392f1..54c6cc77d2 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_cosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_cosl.c @@ -53,9 +53,11 @@ long double __cosl(long double x) { long double y[2],z=0.0L; int64_t n, ix; + double xhi; /* High word of x. */ - GET_LDOUBLE_MSW64(ix,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (ix, xhi); /* |x| ~< pi/4 */ ix &= 0x7fffffffffffffffLL; diff --git a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c index 6a4475ed6b..95dc415845 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c @@ -101,6 +101,7 @@ * erfc/erf(NaN) is NaN */ +#include <errno.h> #include <math.h> #include <math_private.h> #include <math_ldbl_opt.h> @@ -760,16 +761,16 @@ long double __erfl (long double x) { long double a, y, z; - int32_t i, ix, sign; - ieee854_long_double_shape_type u; + int32_t i, ix, hx; + double xhi; - u.value = x; - sign = u.parts32.w0; - ix = sign & 0x7fffffff; + xhi = ldbl_high (x); + GET_HIGH_WORD (hx, xhi); + ix = hx & 0x7fffffff; if (ix >= 0x7ff00000) { /* erf(nan)=nan */ - i = ((sign & 0xfff00000) >> 31) << 1; + i = ((uint32_t) hx >> 31) << 1; return (long double) (1 - i) + one / x; /* erf(+-inf)=+-1 */ } @@ -778,7 +779,7 @@ __erfl (long double x) if (ix >= 0x4039A0DE) { /* __erfcl (x) underflows if x > 25.6283 */ - if (sign) + if ((hx & 0x80000000) == 0) return one-tiny; else return tiny-one; @@ -789,8 +790,9 @@ __erfl (long double x) return (one - y); } } - u.parts32.w0 = ix; - a = u.value; + a = x; + if ((hx & 0x80000000) != 0) + a = -a; z = x * x; if (ix < 0x3fec0000) /* a < 0.875 */ { @@ -814,7 +816,7 @@ __erfl (long double x) y = erf_const + neval (a, TN2, NTN2) / deval (a, TD2, NTD2); } - if (sign & 0x80000000) /* x < 0 */ + if (hx & 0x80000000) /* x < 0 */ y = -y; return( y ); } @@ -824,18 +826,18 @@ long double __erfcl (long double x) { long double y, z, p, r; - int32_t i, ix, sign; - ieee854_long_double_shape_type u; + int32_t i, ix; + uint32_t hx; + double xhi; - u.value = x; - sign = u.parts32.w0; - ix = sign & 0x7fffffff; - u.parts32.w0 = ix; + xhi = ldbl_high (x); + GET_HIGH_WORD (hx, xhi); + ix = hx & 0x7fffffff; if (ix >= 0x7ff00000) { /* erfc(nan)=nan */ /* erfc(+-inf)=0,2 */ - return (long double) (((u_int32_t) sign >> 31) << 1) + one / x; + return (long double) ((hx >> 31) << 1) + one / x; } if (ix < 0x3fd00000) /* |x| <1/4 */ @@ -846,7 +848,8 @@ __erfcl (long double x) } if (ix < 0x3ff40000) /* 1.25 */ { - x = u.value; + if ((hx & 0x80000000) != 0) + x = -x; i = 8.0 * x; switch (i) { @@ -891,7 +894,7 @@ __erfcl (long double x) y += C20a; break; } - if (sign & 0x80000000) + if (hx & 0x80000000) y = 2.0L - y; return y; } @@ -899,10 +902,11 @@ __erfcl (long double x) if (ix < 0x405ac000) { /* x < -9 */ - if ((ix >= 0x40220000) && (sign & 0x80000000)) + if (hx >= 0xc0220000) return two - tiny; - x = fabsl (x); + if ((hx & 0x80000000) != 0) + x = -x; z = one / (x * x); i = 8.0 / x; switch (i) @@ -933,22 +937,26 @@ __erfcl (long double x) p = neval (z, RNr8, NRNr8) / deval (z, RDr8, NRDr8); break; } - u.value = x; - u.parts32.w3 = 0; - u.parts32.w2 = 0; - u.parts32.w1 &= 0xf8000000; - z = u.value; + z = (float) x; r = __ieee754_expl (-z * z - 0.5625) * __ieee754_expl ((z - x) * (z + x) + p); - if ((sign & 0x80000000) == 0) - return r / x; + if ((hx & 0x80000000) == 0) + { + long double ret = r / x; + if (ret == 0) + __set_errno (ERANGE); + return ret; + } else return two - r / x; } else { - if ((sign & 0x80000000) == 0) - return tiny * tiny; + if ((hx & 0x80000000) == 0) + { + __set_errno (ERANGE); + return tiny * tiny; + } else return two - tiny; } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c b/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c index 8808dcd896..007e785346 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_expm1l.c @@ -92,19 +92,19 @@ long double __expm1l (long double x) { long double px, qx, xx; - int32_t ix, sign; - ieee854_long_double_shape_type u; + int32_t ix, lx, sign; int k; + double xhi; /* Detect infinity and NaN. */ - u.value = x; - ix = u.parts32.w0; + xhi = ldbl_high (x); + EXTRACT_WORDS (ix, lx, xhi); sign = ix & 0x80000000; ix &= 0x7fffffff; if (ix >= 0x7ff00000) { /* Infinity. */ - if (((ix & 0xfffff) | u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0) + if (((ix - 0x7ff00000) | lx) == 0) { if (sign) return -1.0L; @@ -116,7 +116,7 @@ __expm1l (long double x) } /* expm1(+- 0) = +- 0. */ - if ((ix == 0) && (u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0) + if ((ix | lx) == 0) return x; /* Overflow. */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_fabsl.c b/sysdeps/ieee754/ldbl-128ibm/s_fabsl.c index 99146d8021..c801c97065 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_fabsl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_fabsl.c @@ -29,10 +29,16 @@ static char rcsid[] = "$NetBSD: $"; long double __fabsl(long double x) { u_int64_t hx, lx; - GET_LDOUBLE_WORDS64(hx,lx,x); + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); lx = lx ^ ( hx & 0x8000000000000000LL ); hx = hx & 0x7fffffffffffffffLL; - SET_LDOUBLE_WORDS64(x,hx,lx); + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x; } long_double_symbol (libm, __fabsl, fabsl); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_finitel.c b/sysdeps/ieee754/ldbl-128ibm/s_finitel.c index 8edb34154d..7b4655fadb 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_finitel.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_finitel.c @@ -29,10 +29,14 @@ static char rcsid[] = "$NetBSD: $"; int ___finitel (long double x) { - int64_t hx; - GET_LDOUBLE_MSW64(hx,x); - return (int)((u_int64_t)((hx&0x7fffffffffffffffLL) - -0x7ff0000000000000LL)>>63); + uint64_t hx; + double xhi; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); + hx &= 0x7fffffffffffffffLL; + hx -= 0x7ff0000000000000LL; + return hx >> 63; } hidden_ver (___finitel, __finitel) weak_alias (___finitel, ____finitel) diff --git a/sysdeps/ieee754/ldbl-128ibm/s_fpclassifyl.c b/sysdeps/ieee754/ldbl-128ibm/s_fpclassifyl.c index f4a90b08c7..90586e822e 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_fpclassifyl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_fpclassifyl.c @@ -46,8 +46,10 @@ ___fpclassifyl (long double x) { u_int64_t hx, lx; int retval = FP_NORMAL; + double xhi, xlo; - GET_LDOUBLE_WORDS64 (hx, lx, x); + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); if ((hx & 0x7ff0000000000000ULL) == 0x7ff0000000000000ULL) { /* +/-NaN or +/-Inf */ if (hx & 0x000fffffffffffffULL) { @@ -65,6 +67,7 @@ ___fpclassifyl (long double x) retval = FP_NORMAL; } else { if ((hx & 0x7ff0000000000000ULL) == 0x0360000000000000ULL) { + EXTRACT_WORDS64 (lx, xlo); if ((lx & 0x7fffffffffffffff) /* lower is non-zero */ && ((lx^hx) & 0x8000000000000000ULL)) { /* and sign differs */ /* +/- denormal */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c b/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c index 3ac5374116..7e40663fd4 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c @@ -36,16 +36,21 @@ two107 = 162259276829213363391578010288128.0; /* 0x4670000000000000, 0 */ long double __frexpl(long double x, int *eptr) { - u_int64_t hx, lx, ix, ixl; + uint64_t hx, lx, ix, ixl; int64_t explo; - GET_LDOUBLE_WORDS64(hx,lx,x); + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); ixl = 0x7fffffffffffffffULL&lx; ix = 0x7fffffffffffffffULL&hx; *eptr = 0; - if(ix>=0x7ff0000000000000ULL||((ix|ixl)==0)) return x; /* 0,inf,nan */ + if(ix>=0x7ff0000000000000ULL||ix==0) return x; /* 0,inf,nan */ if (ix<0x0010000000000000ULL) { /* subnormal */ x *= two107; - GET_LDOUBLE_MSW64(hx,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); ix = hx&0x7fffffffffffffffULL; *eptr = -107; } @@ -54,7 +59,7 @@ long double __frexpl(long double x, int *eptr) if (ixl != 0ULL) { explo = (ixl>>52) - (ix>>52) + 0x3fe; if ((ixl&0x7ff0000000000000ULL) == 0LL) { - /* the lower double is a denomal so we need to correct its + /* the lower double is a denormal so we need to correct its mantissa and perhaps its exponent. */ int cnt; @@ -73,7 +78,9 @@ long double __frexpl(long double x, int *eptr) lx = 0ULL; hx = (hx&0x800fffffffffffffULL) | 0x3fe0000000000000ULL; - SET_LDOUBLE_WORDS64(x,hx,lx); + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x; } #ifdef IS_IN_libm diff --git a/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c b/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c index c8dd9ff98a..54e72c9166 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c @@ -1,6 +1,7 @@ /* * __isinf_nsl(x) returns != 0 if x is ±inf, else 0; * no branching! + * slightly dodgy in relying on signed shift right copying sign bit */ #include <math.h> @@ -9,8 +10,14 @@ int __isinf_nsl (long double x) { - int64_t hx,lx; - GET_LDOUBLE_WORDS64(hx,lx,x); - return !((lx & 0x7fffffffffffffffLL) - | ((hx & 0x7fffffffffffffffLL) ^ 0x7ff0000000000000LL)); + double xhi; + int64_t hx, mask; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); + + mask = (hx & 0x7fffffffffffffffLL) ^ 0x7ff0000000000000LL; + mask |= -mask; + mask >>= 63; + return ~mask; } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_isinfl.c b/sysdeps/ieee754/ldbl-128ibm/s_isinfl.c index 5f5b0144b2..6a728221fc 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_isinfl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_isinfl.c @@ -11,6 +11,7 @@ static char rcsid[] = "$NetBSD: $"; /* * isinfl(x) returns 1 if x is inf, -1 if x is -inf, else 0; * no branching! + * slightly dodgy in relying on signed shift right copying sign bit */ #include <math.h> @@ -20,12 +21,16 @@ static char rcsid[] = "$NetBSD: $"; int ___isinfl (long double x) { - int64_t hx,lx; - GET_LDOUBLE_WORDS64(hx,lx,x); - lx = (lx & 0x7fffffffffffffffLL); - lx |= (hx & 0x7fffffffffffffffLL) ^ 0x7ff0000000000000LL; - lx |= -lx; - return ~(lx >> 63) & (hx >> 62); + double xhi; + int64_t hx, mask; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); + + mask = (hx & 0x7fffffffffffffffLL) ^ 0x7ff0000000000000LL; + mask |= -mask; + mask >>= 63; + return ~mask & (hx >> 62); } hidden_ver (___isinfl, __isinfl) #ifndef IS_IN_libm diff --git a/sysdeps/ieee754/ldbl-128ibm/s_isnanl.c b/sysdeps/ieee754/ldbl-128ibm/s_isnanl.c index 264dec745e..d12f1d3bf5 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_isnanl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_isnanl.c @@ -29,12 +29,14 @@ static char rcsid[] = "$NetBSD: $"; int ___isnanl (long double x) { - int64_t hx; - int64_t lx __attribute__ ((unused)); - GET_LDOUBLE_WORDS64(hx,lx,x); - hx &= 0x7fffffffffffffffLL; - hx = 0x7ff0000000000000LL - hx; - return (int)((u_int64_t)hx>>63); + uint64_t hx; + double xhi; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); + hx &= 0x7fffffffffffffffLL; + hx = 0x7ff0000000000000LL - hx; + return (int) (hx >> 63); } hidden_ver (___isnanl, __isnanl) #ifndef IS_IN_libm diff --git a/sysdeps/ieee754/ldbl-128ibm/s_issignalingl.c b/sysdeps/ieee754/ldbl-128ibm/s_issignalingl.c index 96fab1aff1..bdd58f8f25 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_issignalingl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_issignalingl.c @@ -22,10 +22,13 @@ int __issignalingl (long double x) { - u_int64_t xi; + uint64_t xi; /* For inspecting NaN status, we only have to look at the first of the pair of IEEE 754 64-bit precision numbers. */ - GET_LDOUBLE_MSW64 (xi, x); + double xhi; + + xhi = ldbl_high (x); + EXTRACT_WORDS64 (xi, xhi); #ifdef HIGH_ORDER_BIT_IS_SET_FOR_SNAN # error untested /* We only have to care about the high-order bit of x's significand, because diff --git a/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c b/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c index 8560349631..35039737bf 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_llrintl.c @@ -43,7 +43,7 @@ __llrintl (long double x) #endif ) { - save_round = fegetround (); + save_round = __fegetround (); if (__builtin_expect ((xh == -(double) (-__LONG_LONG_MAX__ - 1)), 0)) { diff --git a/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c b/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c index 77c4fdea84..a346383052 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c @@ -126,19 +126,18 @@ long double __log1pl (long double xm1) { long double x, y, z, r, s; - ieee854_long_double_shape_type u; - int32_t hx; + double xhi; + int32_t hx, lx; int e; /* Test for NaN or infinity input. */ - u.value = xm1; - hx = u.parts32.w0; + xhi = ldbl_high (xm1); + EXTRACT_WORDS (hx, lx, xhi); if (hx >= 0x7ff00000) return xm1; /* log1p(+- 0) = +- 0. */ - if (((hx & 0x7fffffff) == 0) - && (u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) == 0) + if (((hx & 0x7fffffff) | lx) == 0) return xm1; x = xm1 + 1.0L; diff --git a/sysdeps/ieee754/ldbl-128ibm/s_logbl.c b/sysdeps/ieee754/ldbl-128ibm/s_logbl.c index 6cbfcfa1cc..da8d71bdec 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_logbl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_logbl.c @@ -27,9 +27,10 @@ long double __logbl (long double x) { int64_t hx, rhx; - int64_t lx __attribute__ ((unused)); + double xhi; - GET_LDOUBLE_WORDS64 (hx, lx, x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (hx, xhi); hx &= 0x7fffffffffffffffLL; /* high |x| */ if (hx == 0) return -1.0 / fabs (x); @@ -43,5 +44,6 @@ __logbl (long double x) } return (long double) (rhx - 1023); } - +#ifndef __logbl long_double_symbol (libm, __logbl, logbl); +#endif diff --git a/sysdeps/ieee754/ldbl-128ibm/s_lrintl.c b/sysdeps/ieee754/ldbl-128ibm/s_lrintl.c index 588098d090..49dbd42f5b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_lrintl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_lrintl.c @@ -49,7 +49,7 @@ __lrintl (long double x) #endif ) { - save_round = fegetround (); + save_round = __fegetround (); #if __LONG_MAX__ == 2147483647 long long llhi = (long long) xh; diff --git a/sysdeps/ieee754/ldbl-128ibm/s_modfl.c b/sysdeps/ieee754/ldbl-128ibm/s_modfl.c index 39de9d4bfb..ed03ce236c 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_modfl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_modfl.c @@ -37,43 +37,54 @@ long double __modfl(long double x, long double *iptr) { int64_t i0,i1,j0; u_int64_t i; - GET_LDOUBLE_WORDS64(i0,i1,x); + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (i0, xhi); + EXTRACT_WORDS64 (i1, xlo); i1 &= 0x000fffffffffffffLL; j0 = ((i0>>52)&0x7ff)-0x3ff; /* exponent of x */ if(j0<52) { /* integer part in high x */ if(j0<0) { /* |x|<1 */ /* *iptr = +-0 */ - SET_LDOUBLE_WORDS64(*iptr,i0&0x8000000000000000ULL,0); + INSERT_WORDS64 (xhi, i0&0x8000000000000000ULL); + *iptr = xhi; return x; } else { i = (0x000fffffffffffffLL)>>j0; if(((i0&i)|(i1&0x7fffffffffffffffLL))==0) { /* x is integral */ *iptr = x; /* return +-0 */ - SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0); + INSERT_WORDS64 (xhi, i0&0x8000000000000000ULL); + x = xhi; return x; } else { - SET_LDOUBLE_WORDS64(*iptr,i0&(~i),0); + INSERT_WORDS64 (xhi, i0&(~i)); + *iptr = xhi; return x - *iptr; } } } else if (j0>103) { /* no fraction part */ *iptr = x*one; /* We must handle NaNs separately. */ - if (j0 == 0x400 && ((i0 & 0x000fffffffffffffLL) | i1)) + if ((i0 & 0x7fffffffffffffffLL) > 0x7ff0000000000000LL) return x*one; /* return +-0 */ - SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0); + INSERT_WORDS64 (xhi, i0&0x8000000000000000ULL); + x = xhi; return x; } else { /* fraction part in low x */ i = -1ULL>>(j0-52); if((i1&i)==0) { /* x is integral */ *iptr = x; /* return +-0 */ - SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0); + INSERT_WORDS64 (xhi, i0&0x8000000000000000ULL); + x = xhi; return x; } else { - SET_LDOUBLE_WORDS64(*iptr,i0,i1&(~i)); + INSERT_WORDS64 (xhi, i0); + INSERT_WORDS64 (xlo, i1&(~i)); + *iptr = ldbl_pack (xhi, xlo); return x - *iptr; } } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nearbyintl.c b/sysdeps/ieee754/ldbl-128ibm/s_nearbyintl.c index bfcd11044d..92ced5218b 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_nearbyintl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_nearbyintl.c @@ -34,11 +34,11 @@ __nearbyintl (long double x) fenv_t env; static const long double TWO52 = 4503599627370496.0L; union ibm_extended_long_double u; - u.d = x; + u.ld = x; - if (fabs (u.dd[0]) < TWO52) + if (fabs (u.d[0].d) < TWO52) { - double high = u.dd[0]; + double high = u.d[0].d; feholdexcept (&env); if (high > 0.0) { @@ -52,13 +52,13 @@ __nearbyintl (long double x) high += TWO52; if (high == 0.0) high = -0.0; } - u.dd[0] = high; - u.dd[1] = 0.0; - math_force_eval (u.dd[0]); - math_force_eval (u.dd[1]); + u.d[0].d = high; + u.d[1].d = 0.0; + math_force_eval (u.d[0]); + math_force_eval (u.d[1]); fesetenv (&env); } - else if (fabs (u.dd[1]) < TWO52 && u.dd[1] != 0.0) + else if (fabs (u.d[1].d) < TWO52 && u.d[1].d != 0.0) { double high, low, tau; /* In this case we have to round the low double and handle any @@ -67,57 +67,57 @@ __nearbyintl (long double x) may already be rounded and the low double may have the opposite sign to compensate. */ feholdexcept (&env); - if (u.dd[0] > 0.0) + if (u.d[0].d > 0.0) { - if (u.dd[1] > 0.0) + if (u.d[1].d > 0.0) { /* If the high/low doubles are the same sign then simply round the low double. */ - high = u.dd[0]; - low = u.dd[1]; + high = u.d[0].d; + low = u.d[1].d; } - else if (u.dd[1] < 0.0) + else if (u.d[1].d < 0.0) { /* Else the high double is pre rounded and we need to adjust for that. */ - tau = __nextafter (u.dd[0], 0.0); - tau = (u.dd[0] - tau) * 2.0; - high = u.dd[0] - tau; - low = u.dd[1] + tau; + tau = __nextafter (u.d[0].d, 0.0); + tau = (u.d[0].d - tau) * 2.0; + high = u.d[0].d - tau; + low = u.d[1].d + tau; } low += TWO52; low -= TWO52; } - else if (u.dd[0] < 0.0) + else if (u.d[0].d < 0.0) { - if (u.dd[1] < 0.0) + if (u.d[1].d < 0.0) { /* If the high/low doubles are the same sign then simply round the low double. */ - high = u.dd[0]; - low = u.dd[1]; + high = u.d[0].d; + low = u.d[1].d; } - else if (u.dd[1] > 0.0) + else if (u.d[1].d > 0.0) { /* Else the high double is pre rounded and we need to adjust for that. */ - tau = __nextafter (u.dd[0], 0.0); - tau = (u.dd[0] - tau) * 2.0; - high = u.dd[0] - tau; - low = u.dd[1] + tau; + tau = __nextafter (u.d[0].d, 0.0); + tau = (u.d[0].d - tau) * 2.0; + high = u.d[0].d - tau; + low = u.d[1].d + tau; } low = TWO52 - low; low = -(low - TWO52); } - u.dd[0] = high + low; - u.dd[1] = high - u.dd[0] + low; - math_force_eval (u.dd[0]); - math_force_eval (u.dd[1]); + u.d[0].d = high + low; + u.d[1].d = high - u.d[0].d + low; + math_force_eval (u.d[0]); + math_force_eval (u.d[1]); fesetenv (&env); } - return u.d; + return u.ld; } long_double_symbol (libm, __nearbyintl, nearbyintl); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c b/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c index 7e581274a3..c050944c0c 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c @@ -30,27 +30,28 @@ static char rcsid[] = "$NetBSD: $"; long double __nextafterl(long double x, long double y) { - int64_t hx,hy,ihx,ihy,ilx; - u_int64_t lx; - u_int64_t ly __attribute__ ((unused)); + int64_t hx,hy,ihx,ihy; + uint64_t lx; + double xhi, xlo, yhi; - GET_LDOUBLE_WORDS64(hx,lx,x); - GET_LDOUBLE_WORDS64(hy,ly,y); + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); + yhi = ldbl_high (y); + EXTRACT_WORDS64 (hy, yhi); ihx = hx&0x7fffffffffffffffLL; /* |hx| */ - ilx = lx&0x7fffffffffffffffLL; /* |lx| */ ihy = hy&0x7fffffffffffffffLL; /* |hy| */ - if((((ihx&0x7ff0000000000000LL)==0x7ff0000000000000LL)&& - ((ihx&0x000fffffffffffffLL)!=0)) || /* x is nan */ - (((ihy&0x7ff0000000000000LL)==0x7ff0000000000000LL)&& - ((ihy&0x000fffffffffffffLL)!=0))) /* y is nan */ + if((ihx>0x7ff0000000000000LL) || /* x is nan */ + (ihy>0x7ff0000000000000LL)) /* y is nan */ return x+y; /* signal the nan */ if(x==y) return y; /* x=y, return y */ - if(ihx == 0 && ilx == 0) { /* x == 0 */ - long double u; + if(ihx == 0) { /* x == 0 */ + long double u; /* return +-minsubnormal */ hy = (hy & 0x8000000000000000ULL) | 1; - SET_LDOUBLE_WORDS64(x,hy,0ULL);/* return +-minsubnormal */ + INSERT_WORDS64 (yhi, hy); + x = yhi; u = math_opt_barrier (x); u = u * u; math_force_eval (u); /* raise underflow flag */ @@ -59,10 +60,16 @@ long double __nextafterl(long double x, long double y) long double u; if(x > y) { /* x > y, x -= ulp */ + /* This isn't the largest magnitude correctly rounded + long double as you can see from the lowest mantissa + bit being zero. It is however the largest magnitude + long double with a 106 bit mantissa, and nextafterl + is insane with variable precision. So to make + nextafterl sane we assume 106 bit precision. */ if((hx==0xffefffffffffffffLL)&&(lx==0xfc8ffffffffffffeLL)) return x+x; /* overflow, return -inf */ if (hx >= 0x7ff0000000000000LL) { - SET_LDOUBLE_WORDS64(u,0x7fefffffffffffffLL,0x7c8ffffffffffffeLL); + u = 0x1.fffffffffffff7ffffffffffff8p+1023L; return u; } if(ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */ @@ -77,16 +84,19 @@ long double __nextafterl(long double x, long double y) return x; } if (ihx < 0x06a0000000000000LL) { /* ulp will denormal */ - SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL),0ULL); + INSERT_WORDS64 (yhi, hx & (0x7ffLL<<52)); + u = yhi; u *= 0x1.0000000000000p-105L; - } else - SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL)-0x0690000000000000LL,0ULL); + } else { + INSERT_WORDS64 (yhi, (hx & (0x7ffLL<<52))-(0x069LL<<52)); + u = yhi; + } return x - u; } else { /* x < y, x += ulp */ if((hx==0x7fefffffffffffffLL)&&(lx==0x7c8ffffffffffffeLL)) return x+x; /* overflow, return +inf */ - if ((u_int64_t) hx >= 0xfff0000000000000ULL) { - SET_LDOUBLE_WORDS64(u,0xffefffffffffffffLL,0xfc8ffffffffffffeLL); + if ((uint64_t) hx >= 0xfff0000000000000ULL) { + u = -0x1.fffffffffffff7ffffffffffff8p+1023L; return u; } if(ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */ @@ -103,10 +113,13 @@ long double __nextafterl(long double x, long double y) return x; } if (ihx < 0x06a0000000000000LL) { /* ulp will denormal */ - SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL),0ULL); + INSERT_WORDS64 (yhi, hx & (0x7ffLL<<52)); + u = yhi; u *= 0x1.0000000000000p-105L; - } else - SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL)-0x0690000000000000LL,0ULL); + } else { + INSERT_WORDS64 (yhi, (hx & (0x7ffLL<<52))-(0x069LL<<52)); + u = yhi; + } return x + u; } } diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c b/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c index 7e288a4368..b40cf167f3 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c @@ -34,23 +34,22 @@ double __nexttoward(double x, long double y) { int32_t hx,ix; int64_t hy,iy; - u_int32_t lx; - u_int64_t ly,uly; + uint32_t lx; + double yhi; EXTRACT_WORDS(hx,lx,x); - GET_LDOUBLE_WORDS64(hy,ly,y); + yhi = ldbl_high (y); + EXTRACT_WORDS64(hy,yhi); ix = hx&0x7fffffff; /* |x| */ iy = hy&0x7fffffffffffffffLL; /* |y| */ - uly = ly&0x7fffffffffffffffLL; /* |y| */ if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */ - ((iy>=0x7ff0000000000000LL)&&((iy-0x7ff0000000000000LL)|uly)!=0)) - /* y is nan */ + iy>0x7ff0000000000000LL) /* y is nan */ return x+y; if((long double) x==y) return y; /* x=y, return y */ if((ix|lx)==0) { /* x == 0 */ double u; - INSERT_WORDS(x,(u_int32_t)((hy>>32)&0x80000000),1);/* return +-minsub */ + INSERT_WORDS(x,(uint32_t)((hy>>32)&0x80000000),1);/* return +-minsub */ u = math_opt_barrier (x); u = u * u; math_force_eval (u); /* raise underflow flag */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c b/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c index b387a91192..19522f4762 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c @@ -27,16 +27,16 @@ float __nexttowardf(float x, long double y) { int32_t hx,ix; int64_t hy,iy; - u_int64_t ly, uly; + double yhi; GET_FLOAT_WORD(hx,x); - GET_LDOUBLE_WORDS64(hy,ly,y); + yhi = ldbl_high (y); + EXTRACT_WORDS64 (hy, yhi); ix = hx&0x7fffffff; /* |x| */ iy = hy&0x7fffffffffffffffLL; /* |y| */ - uly = ly&0x7fffffffffffffffLL; /* |y| */ if((ix>0x7f800000) || /* x is nan */ - ((iy>=0x7ff0000000000000LL)&&((iy-0x7ff0000000000000LL)|uly)!=0)) + (iy>0x7ff0000000000000LL)) /* y is nan */ return x+y; if((long double) x==y) return y; /* x=y, return y */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_remquol.c b/sysdeps/ieee754/ldbl-128ibm/s_remquol.c index f4777a0e1a..195e108ca9 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_remquol.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_remquol.c @@ -33,20 +33,24 @@ __remquol (long double x, long double y, int *quo) int64_t hx,hy; u_int64_t sx,lx,ly,qs; int cquo; - - GET_LDOUBLE_WORDS64 (hx, lx, x); - GET_LDOUBLE_WORDS64 (hy, ly, y); + double xhi, xlo, yhi, ylo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); + ldbl_unpack (y, &yhi, &ylo); + EXTRACT_WORDS64 (hy, yhi); + EXTRACT_WORDS64 (ly, ylo); sx = hx & 0x8000000000000000ULL; qs = sx ^ (hy & 0x8000000000000000ULL); hy &= 0x7fffffffffffffffLL; hx &= 0x7fffffffffffffffLL; /* Purge off exception values. */ - if ((hy | (ly & 0x7fffffffffffffff)) == 0) + if (hy == 0) return (x * y) / (x * y); /* y = 0 */ if ((hx >= 0x7ff0000000000000LL) /* x not finite */ - || ((hy >= 0x7ff0000000000000LL) /* y is NaN */ - && (((hy - 0x7ff0000000000000LL) | ly) != 0))) + || (hy > 0x7ff0000000000000LL)) /* y is NaN */ return (x * y) / (x * y); if (hy <= 0x7fbfffffffffffffLL) diff --git a/sysdeps/ieee754/ldbl-128ibm/s_rintl.c b/sysdeps/ieee754/ldbl-128ibm/s_rintl.c index 48dbe8569c..5fd6bb8702 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_rintl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_rintl.c @@ -40,7 +40,7 @@ __rintl (long double x) __builtin_inf ()), 1)) { double orig_xh; - int save_round = fegetround (); + int save_round = __fegetround (); /* Long double arithmetic, including the canonicalisation below, only works in round-to-nearest mode. */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_scalblnl.c b/sysdeps/ieee754/ldbl-128ibm/s_scalblnl.c index d752568885..03d4597271 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_scalblnl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_scalblnl.c @@ -41,11 +41,15 @@ long double __scalblnl (long double x, long int n) { int64_t k,l,hx,lx; union { int64_t i; double d; } u; - GET_LDOUBLE_WORDS64(hx,lx,x); + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); k = (hx>>52)&0x7ff; /* extract exponent */ l = (lx>>52)&0x7ff; if (k==0) { /* 0 or subnormal x */ - if (((hx|lx)&0x7fffffffffffffffULL)==0) return x; /* +-0 */ + if ((hx&0x7fffffffffffffffULL)==0) return x; /* +-0 */ u.i = hx; u.d *= two54; hx = u.i; @@ -61,7 +65,9 @@ long double __scalblnl (long double x, long int n) if (k > 0) { /* normal result */ hx = (hx&0x800fffffffffffffULL)|(k<<52); if ((lx & 0x7fffffffffffffffULL) == 0) { /* low part +-0 */ - SET_LDOUBLE_WORDS64(x,hx,lx); + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x; } if (l == 0) { /* low part subnormal */ @@ -81,14 +87,19 @@ long double __scalblnl (long double x, long int n) u.d *= twom54; lx = u.i; } - SET_LDOUBLE_WORDS64(x,hx,lx); + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x; } if (k <= -54) return tiny*__copysignl(tiny,x); /*underflow*/ k += 54; /* subnormal result */ lx &= 0x8000000000000000ULL; - SET_LDOUBLE_WORDS64(x,(hx&0x800fffffffffffffULL)|(k<<52),lx); + hx &= 0x800fffffffffffffULL; + INSERT_WORDS64 (xhi, hx|(k<<52)); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x*twolm54; } long_double_symbol (libm, __scalblnl, scalblnl); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c b/sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c index bcdb23bdaa..161172db6e 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c @@ -41,11 +41,15 @@ long double __scalbnl (long double x, int n) { int64_t k,l,hx,lx; union { int64_t i; double d; } u; - GET_LDOUBLE_WORDS64(hx,lx,x); + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); k = (hx>>52)&0x7ff; /* extract exponent */ l = (lx>>52)&0x7ff; if (k==0) { /* 0 or subnormal x */ - if (((hx|lx)&0x7fffffffffffffffULL)==0) return x; /* +-0 */ + if ((hx&0x7fffffffffffffffULL)==0) return x; /* +-0 */ u.i = hx; u.d *= two54; hx = u.i; @@ -61,7 +65,9 @@ long double __scalbnl (long double x, int n) if (k > 0) { /* normal result */ hx = (hx&0x800fffffffffffffULL)|(k<<52); if ((lx & 0x7fffffffffffffffULL) == 0) { /* low part +-0 */ - SET_LDOUBLE_WORDS64(x,hx,lx); + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x; } if (l == 0) { /* low part subnormal */ @@ -81,14 +87,19 @@ long double __scalbnl (long double x, int n) u.d *= twom54; lx = u.i; } - SET_LDOUBLE_WORDS64(x,hx,lx); + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x; } if (k <= -54) return tiny*__copysignl(tiny,x); /*underflow*/ k += 54; /* subnormal result */ lx &= 0x8000000000000000ULL; - SET_LDOUBLE_WORDS64(x,(hx&0x800fffffffffffffULL)|(k<<52),lx); + hx &= 0x800fffffffffffffULL; + INSERT_WORDS64 (xhi, hx|(k<<52)); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); return x*twolm54; } #ifdef IS_IN_libm diff --git a/sysdeps/ieee754/ldbl-128ibm/s_signbitl.c b/sysdeps/ieee754/ldbl-128ibm/s_signbitl.c index ee4aea6cf4..aecb1fd792 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_signbitl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_signbitl.c @@ -25,8 +25,10 @@ int ___signbitl (long double x) { int64_t e; + double xhi; - GET_LDOUBLE_MSW64 (e, x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (e, xhi); return e < 0; } #ifdef IS_IN_libm diff --git a/sysdeps/ieee754/ldbl-128ibm/s_sincosl.c b/sysdeps/ieee754/ldbl-128ibm/s_sincosl.c index 3b1e547bdc..a9e2f3d19a 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_sincosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_sincosl.c @@ -27,9 +27,11 @@ void __sincosl (long double x, long double *sinx, long double *cosx) { int64_t ix; + double xhi; /* High word of x. */ - GET_LDOUBLE_MSW64 (ix, x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (ix, xhi); /* |x| ~< pi/4 */ ix &= 0x7fffffffffffffffLL; diff --git a/sysdeps/ieee754/ldbl-128ibm/s_sinl.c b/sysdeps/ieee754/ldbl-128ibm/s_sinl.c index 6fec16f851..087921a913 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_sinl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_sinl.c @@ -53,9 +53,11 @@ long double __sinl(long double x) { long double y[2],z=0.0L; int64_t n, ix; + double xhi; /* High word of x. */ - GET_LDOUBLE_MSW64(ix,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (ix, xhi); /* |x| ~< pi/4 */ ix &= 0x7fffffffffffffffLL; diff --git a/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c b/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c index 138b63cd1a..c63e25345d 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_tanhl.c @@ -47,10 +47,12 @@ static const long double one=1.0L, two=2.0L, tiny = 1.0e-300L; long double __tanhl(long double x) { long double t,z; - int64_t jx,ix,lx; + int64_t jx,ix; + double xhi; /* High word of |x|. */ - GET_LDOUBLE_WORDS64(jx,lx,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (jx, xhi); ix = jx&0x7fffffffffffffffLL; /* x is INF or NaN */ @@ -61,7 +63,7 @@ long double __tanhl(long double x) /* |x| < 22 */ if (ix < 0x4036000000000000LL) { /* |x|<22 */ - if ((ix | (lx&0x7fffffffffffffffLL)) == 0) + if (ix == 0) return x; /* x == +-0 */ if (ix<0x3c60000000000000LL) /* |x|<2**-57 */ return x*(one+x); /* tanh(small) = small */ diff --git a/sysdeps/ieee754/ldbl-128ibm/s_tanl.c b/sysdeps/ieee754/ldbl-128ibm/s_tanl.c index 9967d0c206..66b8a0621e 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_tanl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_tanl.c @@ -53,9 +53,11 @@ long double __tanl(long double x) { long double y[2],z=0.0L; int64_t n, ix; + double xhi; /* High word of x. */ - GET_LDOUBLE_MSW64(ix,x); + xhi = ldbl_high (x); + EXTRACT_WORDS64 (ix, xhi); /* |x| ~< pi/4 */ ix &= 0x7fffffffffffffffLL; diff --git a/sysdeps/ieee754/ldbl-128ibm/strtold_l.c b/sysdeps/ieee754/ldbl-128ibm/strtold_l.c index 04e3288571..93a80c5eec 100644 --- a/sysdeps/ieee754/ldbl-128ibm/strtold_l.c +++ b/sysdeps/ieee754/ldbl-128ibm/strtold_l.c @@ -43,11 +43,11 @@ libc_hidden_proto (STRTOF) #define FLOAT_HUGE_VAL HUGE_VALL # define SET_MANTISSA(flt, mant) \ do { union ibm_extended_long_double u; \ - u.d = (flt); \ - u.ieee_nan.mantissa0 = (mant) >> 32; \ - u.ieee_nan.mantissa1 = (mant); \ - if ((u.ieee.mantissa0 | u.ieee.mantissa1) != 0) \ - (flt) = u.d; \ + u.ld = (flt); \ + u.d[0].ieee_nan.mantissa0 = (mant) >> 32; \ + u.d[0].ieee_nan.mantissa1 = (mant); \ + if ((u.d[0].ieee.mantissa0 | u.d[0].ieee.mantissa1) != 0) \ + (flt) = u.ld; \ } while (0) #include <strtod_l.c> diff --git a/sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c b/sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c index ed0d4a5067..06dcf02ffe 100644 --- a/sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c +++ b/sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c @@ -89,23 +89,23 @@ __x2y2m1l (long double x, long double y) double vals[12]; SET_RESTORE_ROUND (FE_TONEAREST); union ibm_extended_long_double xu, yu; - xu.d = x; - yu.d = y; - if (fabs (xu.dd[1]) < 0x1p-500) - xu.dd[1] = 0.0; - if (fabs (yu.dd[1]) < 0x1p-500) - yu.dd[1] = 0.0; - mul_split (&vals[1], &vals[0], xu.dd[0], xu.dd[0]); - mul_split (&vals[3], &vals[2], xu.dd[0], xu.dd[1]); + xu.ld = x; + yu.ld = y; + if (fabs (xu.d[1].d) < 0x1p-500) + xu.d[1].d = 0.0; + if (fabs (yu.d[1].d) < 0x1p-500) + yu.d[1].d = 0.0; + mul_split (&vals[1], &vals[0], xu.d[0].d, xu.d[0].d); + mul_split (&vals[3], &vals[2], xu.d[0].d, xu.d[1].d); vals[2] *= 2.0; vals[3] *= 2.0; - mul_split (&vals[5], &vals[4], xu.dd[1], xu.dd[1]); - mul_split (&vals[7], &vals[6], yu.dd[0], yu.dd[0]); - mul_split (&vals[9], &vals[8], yu.dd[0], yu.dd[1]); + mul_split (&vals[5], &vals[4], xu.d[1].d, xu.d[1].d); + mul_split (&vals[7], &vals[6], yu.d[0].d, yu.d[0].d); + mul_split (&vals[9], &vals[8], yu.d[0].d, yu.d[1].d); vals[8] *= 2.0; vals[9] *= 2.0; - mul_split (&vals[11], &vals[10], yu.dd[1], yu.dd[1]); - if (xu.dd[0] >= 0.75) + mul_split (&vals[11], &vals[10], yu.d[1].d, yu.d[1].d); + if (xu.d[0].d >= 0.75) vals[1] -= 1.0; else { |