Message ID | 20210804094609.245719-1-Paul.Zimmermann@inria.fr |
---|---|
State | New |
Headers | show |
Series | Fixed inaccuracy of j0f (BZ #28185) | expand |
* Paul Zimmermann: > The largest errors over the full binary32 range are after this > patch (on x86_64): > > RNDN: libm wrong by up to 9.00e+00 ulp(s) [9] for x=0x1.04c39cp+6 > RNDZ: libm wrong by up to 9.00e+00 ulp(s) [9] for x=0x1.04c39cp+6 > RNDU: libm wrong by up to 9.00e+00 ulp(s) [9] for x=0x1.04c39cp+6 > RNDD: libm wrong by up to 8.98e+00 ulp(s) [9] for x=0x1.4b7066p+7 > > Some constants used in the code are also changed to "static const". Would it make sense to add a few previously broken arguments to the test cases? We currently do not see any test failures. Thanks, Florian
Dear Florian, > From: Florian Weimer <fweimer@redhat.com> > Cc: libc-alpha@sourceware.org > Date: Wed, 04 Aug 2021 12:46:00 +0200 > User-Agent: Gnus/5.13 (Gnus v5.13) Emacs/27.2 (gnu/linux) > > * Paul Zimmermann: > > > The largest errors over the full binary32 range are after this > > patch (on x86_64): > > > > RNDN: libm wrong by up to 9.00e+00 ulp(s) [9] for x=0x1.04c39cp+6 > > RNDZ: libm wrong by up to 9.00e+00 ulp(s) [9] for x=0x1.04c39cp+6 > > RNDU: libm wrong by up to 9.00e+00 ulp(s) [9] for x=0x1.04c39cp+6 > > RNDD: libm wrong by up to 8.98e+00 ulp(s) [9] for x=0x1.4b7066p+7 > > > > Some constants used in the code are also changed to "static const". > > Would it make sense to add a few previously broken arguments to the test > cases? We currently do not see any test failures. I could, but I fear it will yield plenty of errors for other formats, since the "broken arguments" are near roots of j0, where glibc code behaves quite poorly. For example if I add in auto-libm-test-in 44 inputs near the first root of j0 I get: $ make -j4 check subdirs=math ... FAIL: math/test-double-j0 FAIL: math/test-float128-j0 FAIL: math/test-float32x-j0 FAIL: math/test-float64-j0 FAIL: math/test-float64x-j0 FAIL: math/test-ldouble-j0 Would it make sense to add just one entry (the one with the largest error) with xfail-rounding:double,float128,float32x,float64,float64x,ldouble? Paul
On Wed, 4 Aug 2021, Paul Zimmermann wrote: > Would it make sense to add just one entry (the one with the largest error) > with xfail-rounding:double,float128,float32x,float64,float64x,ldouble? xfail-rounding is only if the problems are only in non-default rounding modes, if there are large errors in round-to-nearest you should use xfail (and note that the arguments are the format names as used by gen-auto-libm-tests, e.g. binary64, not the names used in the testcase names).
thank you Joseph. I will submit a new patch with lines like: j0 0x1.31ec02p+1 xfail:binary64 xfail:intel96 xfail:binary128 Paul > Date: Wed, 4 Aug 2021 17:36:45 +0000 > From: Joseph Myers <joseph@codesourcery.com> > > On Wed, 4 Aug 2021, Paul Zimmermann wrote: > > > Would it make sense to add just one entry (the one with the largest error) > > with xfail-rounding:double,float128,float32x,float64,float64x,ldouble? > > xfail-rounding is only if the problems are only in non-default rounding > modes, if there are large errors in round-to-nearest you should use xfail > (and note that the arguments are the format names as used by > gen-auto-libm-tests, e.g. binary64, not the names used in the testcase > names). > > -- > Joseph S. Myers > joseph@codesourcery.com
diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c index 462518c365..080395eab5 100644 --- a/sysdeps/ieee754/flt-32/e_j0f.c +++ b/sysdeps/ieee754/flt-32/e_j0f.c @@ -40,7 +40,7 @@ S04 = 1.1661400734e-09; /* 0x30a045e8 */ static const float zero = 0.0; /* This is the nearest approximation of the first zero of j0. */ -#define FIRST_ZERO_J0 0xf.26247p-28f +#define FIRST_ZERO_J0 0x1.33d152p+1f #define SMALL_SIZE 64 @@ -212,7 +212,7 @@ j0f_asympt (float x) /* Now x - pi/4 - alpha0 = h + n*pi/2 mod (2*pi). */ float xr = (float) h; n = n & 3; - float cst = 0xc.c422ap-4; /* sqrt(2/pi) rounded to nearest */ + static const float cst = 0xc.c422ap-4f; /* sqrt(2/pi) rounded to nearest */ float t = cst / sqrtf (x) * (float) beta0; if (n == 0) return t * __cosf (xr); @@ -286,7 +286,7 @@ __ieee754_j0f(float x) /* The following threshold is optimal: for x=0x1.3b58dep+1 and rounding upwards, |cc|=0x1.579b26p-4 and z is 10 ulps far from the correctly rounded value. */ - float threshold = 0x1.579b26p-4; + static const float threshold = 0x1.579b26p-4f; if (fabsf (cc) > threshold) return z; else @@ -493,7 +493,7 @@ y0f_asympt (float x) /* Now x - pi/4 - alpha0 = h + n*pi/2 mod (2*pi). */ float xr = (float) h; n = n & 3; - float cst = 0xc.c422ap-4; /* sqrt(2/pi) rounded to nearest */ + static const float cst = 0xc.c422ap-4; /* sqrt(2/pi) rounded to nearest */ float t = cst / sqrtf (x) * (float) beta0; if (n == 0) return t * __sinf (xr); @@ -584,7 +584,7 @@ __ieee754_y0f(float x) z = (invsqrtpi*ss)/sqrtf(x); /* The following threshold is optimal (determined on aarch64-linux-gnu). */ - float threshold = 0x1.be585ap-4; + static const float threshold = 0x1.be585ap-4; if (fabsf (ss) > threshold) return z; else