diff options
| -rw-r--r-- | ChangeLog | 27 | ||||
| -rw-r--r-- | math/libm-test.inc | 19 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/e_fmodl.c | 125 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/e_hypotl.c | 88 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/e_remainderl.c | 18 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_erfl.c | 55 |
6 files changed, 198 insertions, 134 deletions
diff --git a/ChangeLog b/ChangeLog index 829d4334a8..e3647f33ef 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,32 @@ 2013-10-04 Alan Modra <amodra@gmail.com> + [BZ #15734], [BZ #15735] + * sysdeps/ieee754/ldbl-128ibm/e_fmodl.c (__ieee754_fmodl): Rewrite + all uses of ieee875 long double macros and unions. Simplify test + for 0.0L. Correct |x|<|y| and |x|=|y| test. Use + ldbl_extract_mantissa value for ix,iy exponents. Properly + normalize after ldbl_extract_mantissa, and don't add hidden bit + already handled. Don't treat low word of ieee854 mantissa like + low word of IBM long double and mask off bit when testing for + zero. + * sysdeps/ieee754/ldbl-128ibm/e_hypotl.c (__ieee754_hypotl): Rewrite + all uses of ieee875 long double macros and unions. Simplify tests + for 0.0L and inf. Correct double adjustment of k. Delete dead code + adjusting ha,hb. Simplify code setting kld. Delete two600 and + two1022, instead use their values. Recognise that tests for large + "a" and small "b" are mutually exclusive. Rename vars. Comment. + * sysdeps/ieee754/ldbl-128ibm/e_remainderl.c (__ieee754_remainderl): + Rewrite all uses of ieee875 long double macros and unions. Simplify + test for 0.0L and nan. Correct negation. + * sysdeps/ieee754/ldbl-128ibm/s_erfl.c (__erfl): Rewrite all uses of + ieee875 long double macros and unions. Correct output for large + magnitude x. Correct absolute value calculation. + (__erfcl): Likewise. + * math/libm-test.inc: Add tests for errors discovered in IBM long + double versions of fmodl, remainderl, erfl and erfcl. + +2013-10-04 Alan Modra <amodra@gmail.com> + * sysdeps/ieee754/ldbl-128ibm/e_atan2l.c (__ieee754_atan2l): Rewrite all uses of ieee854 long double macros and unions. Simplify tests for long doubles that are fully specified by the high double. diff --git a/math/libm-test.inc b/math/libm-test.inc index 2cb3d2c614..9677052e00 100644 --- a/math/libm-test.inc +++ b/math/libm-test.inc @@ -7829,6 +7829,11 @@ static const struct test_f_f_data erf_test_data[] = TEST_f_f (erf, 2.0L, 0.995322265018952734162069256367252929L), TEST_f_f (erf, 4.125L, 0.999999994576599200434933994687765914L), TEST_f_f (erf, 27.0L, 1.0L), + TEST_f_f (erf, -27.0L, -1.0L), +#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 54 + /* The input is not exactly representable as a double. */ + TEST_f_f (erf, -0x1.fffffffffffff8p-2L, -0.5204998778130465132916303345518417673509L), +#endif }; static void @@ -7857,6 +7862,10 @@ static const struct test_f_f_data erfc_test_data[] = TEST_f_f (erfc, 0x1.ffa002p+2L, 1.233585992097580296336099501489175967033e-29L), TEST_f_f (erfc, 0x1.ffffc8p+2L, 1.122671365033056305522366683719541099329e-29L), #ifdef TEST_LDOUBLE +# if LDBL_MANT_DIG >= 54 + /* The input is not exactly representable as a double. */ + TEST_f_f (erfc, -0x1.fffffffffffff8p-2L, 1.52049987781304651329163033455184176735L), +# endif /* The result can only be represented in long double. */ # if LDBL_MIN_10_EXP < -319 TEST_f_f (erfc, 27.0L, 0.523704892378925568501606768284954709e-318L), @@ -9355,6 +9364,13 @@ static const struct test_ff_f_data fmod_test_data[] = #if defined TEST_LDOUBLE && LDBL_MIN_EXP <= -16381 TEST_ff_f (fmod, 0x0.fffffffffffffffep-16382L, 0x1p-16445L, plus_zero, NO_INEXACT_EXCEPTION), #endif +#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 56 + TEST_ff_f (fmod, -0x1.00000000000004p+0L, 0x1.fffffffffffff8p-1L, -0x1p-53L, NO_INEXACT_EXCEPTION), + TEST_ff_f (fmod, 0x1.fffffffffffffap-1L, 0x1.fffffffffffff8p-1L, 0x1p-56L, NO_INEXACT_EXCEPTION), + TEST_ff_f (fmod, -0x1.fffffffffffffap-1L, 0x1.fffffffffffff8p-1L, -0x1p-56L, NO_INEXACT_EXCEPTION), + TEST_ff_f (fmod, 0x1.fffffffffffffap-1L, -0x1.fffffffffffff8p-1L, 0x1p-56L, NO_INEXACT_EXCEPTION), + TEST_ff_f (fmod, -0x1.fffffffffffffap-1L, -0x1.fffffffffffff8p-1L, -0x1p-56L, NO_INEXACT_EXCEPTION), +#endif }; static void @@ -12316,6 +12332,9 @@ static const struct test_ff_f_data remainder_test_data[] = TEST_ff_f (remainder, -1.625, -1.0, 0.375, NO_INEXACT_EXCEPTION), TEST_ff_f (remainder, 5.0, 2.0, 1.0, NO_INEXACT_EXCEPTION), TEST_ff_f (remainder, 3.0, 2.0, -1.0, NO_INEXACT_EXCEPTION), +#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 56 + TEST_ff_f (remainder, -0x1.80000000000002p1L, 2.0, 0x1.fffffffffffff8p-1L, NO_INEXACT_EXCEPTION), +#endif }; static void 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_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_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/s_erfl.c b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c index 6a4475ed6b..c861c65ccf 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c @@ -760,16 +760,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 +778,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 +789,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 +815,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 +825,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 +847,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 +893,7 @@ __erfcl (long double x) y += C20a; break; } - if (sign & 0x80000000) + if (hx & 0x80000000) y = 2.0L - y; return y; } @@ -899,10 +901,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,21 +936,17 @@ __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) + if ((hx & 0x80000000) == 0) return r / x; else return two - r / x; } else { - if ((sign & 0x80000000) == 0) + if ((hx & 0x80000000) == 0) return tiny * tiny; else return two - tiny; |