[PATCH] Fixed inaccuracy of j0f (BZ #28185)

[PATCH] Fixed inaccuracy of j0f (BZ #28185)

[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".
---
 sysdeps/ieee754/flt-32/e_j0f.c | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

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
-- 
2.30.2

Re: [PATCH] Fixed inaccuracy of j0f (BZ #28185)

[Florian Weimer]

* 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

Re: [PATCH] Fixed inaccuracy of j0f (BZ #28185)

[Paul Zimmermann]

       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

Re: [PATCH] Fixed inaccuracy of j0f (BZ #28185)

[Joseph Myers]

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

Re: [PATCH] Fixed inaccuracy of j0f (BZ #28185)

[Paul Zimmermann]

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