This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: PowerPC floating point little-endian [4 of 15]
- From: Alan Modra <amodra at gmail dot com>
- To: "Joseph S. Myers" <joseph at codesourcery dot com>
- Cc: libc-alpha at sourceware dot org
- Date: Sat, 13 Jul 2013 01:18:48 +0930
- Subject: Re: PowerPC floating point little-endian [4 of 15]
- References: <20130710012435 dot GN2602 at bubble dot grove dot modra dot org> <20130710012622 dot GQ2602 at bubble dot grove dot modra dot org> <Pine dot LNX dot 4 dot 64 dot 1307101718380 dot 13074 at digraph dot polyomino dot org dot uk>
On Wed, Jul 10, 2013 at 05:26:49PM +0000, Joseph S. Myers wrote:
> When fixing user-visible bugs, please file them in Bugzilla, and if
> possible include testcases for them (that fail before, pass after the
> patch) in libm-test.inc.
OK, I have some libm-test.inc additions.
> My guess is that fmod/remainder tests would be conditional on "#if defined
> TEST_LDOUBLE && LDBL_MANT_DIG >= 106"; I'm less sure about erf/erfc, but
> presumably by knowing where the error is, you should be able to make the
> incorrect negation give errors of up to around 2^52 ulps, and cases with
> large errors make suitable testcases. The erf issue with large negative
> input is certainly easy to test - libm-test.inc is currently lacking any
> coverage of finite negative values.
>
> > if(__builtin_expect(hx<=hy,0)) {
> > - if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
> > + if (hx < hy
> > + || (int64_t) (lx ^ sx) < (int64_t) (ly ^ sy))
> > + return x; /* |x|<|y| return x */
> > if(lx==ly)
> > - return Zero[(u_int64_t)sx>>63]; /* |x|=|y| return x*0*/
> > + return Zero[(uint64_t)sx>>63]; /* |x|=|y| return x*0*/
>
> It looks to me like the (lx==ly) conditional needs to check (lx^sx) and
> (ly^sy) instead to be correct, further illustrating the value of a range
> of new testcases.
What's more, the new testcases revealed that my change wasn't even
correct for the |x|<|y| test. The following is correct I think, and
may even be optimal.
if (hx < hy
|| (((ly ^ sy) & 0x8000000000000000LL) == 0
&& (int64_t) (lx ^ sx) < (int64_t) (ly ^ sy))
|| (((lx ^ sx) & 0x8000000000000000LL) != 0
&& (int64_t) (lx ^ sx) > (int64_t) (ly ^ sy)))
return x; /* |x|<|y| return x */
if ((lx ^ sx) == (ly ^ sy))
return Zero[(uint64_t)sx>>63]; /* |x|=|y| return x*0*/
However, I found that fmodl still didn't give correct results due to
using the exponents before ldbl_extract_mantissa makes an adjustment
for the "borrow from the hidden bit" case. I also found it was
possible to make fmodl hang due to a bogus mask of a significant bit
before the final nomalization.
* math/libm-test.inc: Add erf, erfc, fmod and remainder tests.
* 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.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 2e5237b..fb68e1e 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -7771,6 +7771,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
@@ -7799,6 +7804,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),
@@ -9293,6 +9302,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
@@ -12193,6 +12209,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 a60963c..ffb3bec 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
@@ -27,76 +27,65 @@ 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 (hx < hy
+ || (((ly ^ sy) & 0x8000000000000000LL) == 0
+ && (int64_t) (lx ^ sx) < (int64_t) (ly ^ sy))
+ || (((lx ^ sx) & 0x8000000000000000LL) != 0
+ && (int64_t) (lx ^ sx) > (int64_t) (ly ^ sy)))
+ return x; /* |x|<|y| return x */
+ if ((lx ^ sx) == (ly ^ sy))
+ return Zero[(uint64_t)sx>>63]; /* |x|=|y| return x*0*/
}
- /* 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 +93,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 +102,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;
--
Alan Modra
Australia Development Lab, IBM