Fixed inaccuracy of j0f (BZ #28185)
Checks
Context |
Check |
Description |
dj/TryBot-apply_patch |
success
|
Patch applied to master at the time it was sent
|
dj/TryBot-32bit |
success
|
Build for i686
|
Commit Message
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".
---
sysdeps/ieee754/flt-32/e_j0f.c | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
Comments
* 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
@@ -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