fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny

fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny

[Paul Zimmermann]

       Hi,

before the patch below, the maximum ulp error for j0 in the whole binary32
range is 6177902 ulps (for x = 3.153646966e+38).

After this patch, it is 900691 ulps (for x = 2.404825449e+00).

The patch fixes the case where x >= 2^127 and tiny sin(x)+cos(x).

Large remaining errors are due to a cancellation in another branch of the code.

Paul

PS: the same method can be applied to j1 and y1.
PS2: this can wait for 2.33 of course.

From 6b731f36b1a5badf4704645d0dda40957cedd0db Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Mon, 27 Jul 2020 19:01:18 +0200
Subject: [PATCH] fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is
 tiny

---
 sysdeps/ieee754/flt-32/e_j0f.c | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index c89b9f2688..f85d8a59e0 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -56,6 +56,22 @@ __ieee754_j0f(float x)
 		    if ((s*c)<zero) cc = z/ss;
 		    else	    ss = z/cc;
 		}
+                else {
+                  /* we subtract (exactly) a value x0 such that cos(x0)+sin(x0)
+                     is very near from 0, and use the identity
+                     sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
+                     sin(x) + cos(x) with extra accuracy */
+                   float x0 = 3.153646966e+38f;
+                   float y = x - x0; /* exact */
+                   /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
+                   z = __sinf (y);
+                   float eps = 8.17583368e-8f;
+                   /* cos(x0) ~ -sin(x0) + eps */
+                   z += eps * __cosf (x);
+                   /* now z ~ (sin(x)-cos(x))*cos(x0) */
+                   float cosx0 = -0.707106740f;
+                   cc = z / cosx0;
+                }
 	/*
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
-- 
2.27.0

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny

[Joseph Myers]

On Mon, 27 Jul 2020, Paul Zimmermann wrote:

> +                   float x0 = 3.153646966e+38f;
> +                   float y = x - x0; /* exact */
> +                   /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
> +                   z = __sinf (y);
> +                   float eps = 8.17583368e-8f;
> +                   /* cos(x0) ~ -sin(x0) + eps */
> +                   z += eps * __cosf (x);
> +                   /* now z ~ (sin(x)-cos(x))*cos(x0) */
> +                   float cosx0 = -0.707106740f;

In new code we generally prefer to use hex float constants in such cases 
where a specific floating-point value is wanted.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny

[Paul Zimmermann]

> In new code we generally prefer to use hex float constants in such cases 
> where a specific floating-point value is wanted.

thank you Joseph. Here is a new version. The maximal error for x >= 2^127
is now 4 ulps (attained for x=1.740713465e+38).

Total: errors=4220511 (0.10%) errors2=393216 maxerr=4 ulp(s)

Paul

From 6b731f36b1a5badf4704645d0dda40957cedd0db Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Mon, 27 Jul 2020 19:01:18 +0200
Subject: [PATCH 1/2] fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x)
 is tiny

---
 sysdeps/ieee754/flt-32/e_j0f.c | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index c89b9f2688..f85d8a59e0 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -56,6 +56,22 @@ __ieee754_j0f(float x)
 		    if ((s*c)<zero) cc = z/ss;
 		    else	    ss = z/cc;
 		}
+                else {
+                  /* we subtract (exactly) a value x0 such that cos(x0)+sin(x0)
+                     is very near from 0, and use the identity
+                     sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
+                     sin(x) + cos(x) with extra accuracy */
+                   float x0 = 3.153646966e+38f;
+                   float y = x - x0; /* exact */
+                   /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
+                   z = __sinf (y);
+                   float eps = 8.17583368e-8f;
+                   /* cos(x0) ~ -sin(x0) + eps */
+                   z += eps * __cosf (x);
+                   /* now z ~ (sin(x)-cos(x))*cos(x0) */
+                   float cosx0 = -0.707106740f;
+                   cc = z / cosx0;
+                }
 	/*
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
-- 
2.27.0

From 44124c42fe519c7dcac829160181ba0bb6c8751c Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Tue, 28 Jul 2020 10:05:38 +0200
Subject: [PATCH 2/2] use hex float constants as advised by Joseph Myers

---
 sysdeps/ieee754/flt-32/e_j0f.c | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index f85d8a59e0..883a1a4d13 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -61,15 +61,15 @@ __ieee754_j0f(float x)
                      is very near from 0, and use the identity
                      sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
                      sin(x) + cos(x) with extra accuracy */
-                   float x0 = 3.153646966e+38f;
+                   float x0 = 0xe.d4108p+124f;
                    float y = x - x0; /* exact */
                    /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
                    z = __sinf (y);
-                   float eps = 8.17583368e-8f;
+                   float eps = 0x1.5f263ep-24f;
                    /* cos(x0) ~ -sin(x0) + eps */
                    z += eps * __cosf (x);
                    /* now z ~ (sin(x)-cos(x))*cos(x0) */
-                   float cosx0 = -0.707106740f;
+                   float cosx0 = -0xb.504f3p-4f;
                    cc = z / cosx0;
                 }
 	/*
-- 
2.27.0

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny

[Andreas Schwab]

On Jul 28 2020, Paul Zimmermann wrote:

> +                  /* we subtract (exactly) a value x0 such that cos(x0)+sin(x0)
> +                     is very near from 0, and use the identity

Did you mean "near to"?

Andreas.

-- 
Andreas Schwab, schwab@linux-m68k.org
GPG Key fingerprint = 7578 EB47 D4E5 4D69 2510  2552 DF73 E780 A9DA AEC1
"And now for something completely different."

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny

[Paul Zimmermann]

       Dear Andreas,

yes thanks. Sorry my english is not perfect.

Paul

From 6b731f36b1a5badf4704645d0dda40957cedd0db Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Mon, 27 Jul 2020 19:01:18 +0200
Subject: [PATCH 1/3] fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x)
 is tiny

---
 sysdeps/ieee754/flt-32/e_j0f.c | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index c89b9f2688..f85d8a59e0 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -56,6 +56,22 @@ __ieee754_j0f(float x)
 		    if ((s*c)<zero) cc = z/ss;
 		    else	    ss = z/cc;
 		}
+                else {
+                  /* we subtract (exactly) a value x0 such that cos(x0)+sin(x0)
+                     is very near from 0, and use the identity
+                     sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
+                     sin(x) + cos(x) with extra accuracy */
+                   float x0 = 3.153646966e+38f;
+                   float y = x - x0; /* exact */
+                   /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
+                   z = __sinf (y);
+                   float eps = 8.17583368e-8f;
+                   /* cos(x0) ~ -sin(x0) + eps */
+                   z += eps * __cosf (x);
+                   /* now z ~ (sin(x)-cos(x))*cos(x0) */
+                   float cosx0 = -0.707106740f;
+                   cc = z / cosx0;
+                }
 	/*
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
-- 
2.27.0

From 44124c42fe519c7dcac829160181ba0bb6c8751c Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Tue, 28 Jul 2020 10:05:38 +0200
Subject: [PATCH 2/3] use hex float constants as advised by Joseph Myers

---
 sysdeps/ieee754/flt-32/e_j0f.c | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index f85d8a59e0..883a1a4d13 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -61,15 +61,15 @@ __ieee754_j0f(float x)
                      is very near from 0, and use the identity
                      sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
                      sin(x) + cos(x) with extra accuracy */
-                   float x0 = 3.153646966e+38f;
+                   float x0 = 0xe.d4108p+124f;
                    float y = x - x0; /* exact */
                    /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
                    z = __sinf (y);
-                   float eps = 8.17583368e-8f;
+                   float eps = 0x1.5f263ep-24f;
                    /* cos(x0) ~ -sin(x0) + eps */
                    z += eps * __cosf (x);
                    /* now z ~ (sin(x)-cos(x))*cos(x0) */
-                   float cosx0 = -0.707106740f;
+                   float cosx0 = -0xb.504f3p-4f;
                    cc = z / cosx0;
                 }
 	/*
-- 
2.27.0

From 409781259c5d150ec62079e608ea94b520eafa58 Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Tue, 28 Jul 2020 12:47:31 +0200
Subject: [PATCH 3/3] fixed typo (thanks Andreas Schwab)

---
 sysdeps/ieee754/flt-32/e_j0f.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index 883a1a4d13..66249d2812 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -58,7 +58,7 @@ __ieee754_j0f(float x)
 		}
                 else {
                   /* we subtract (exactly) a value x0 such that cos(x0)+sin(x0)
-                     is very near from 0, and use the identity
+                     is very near to 0, and use the identity
                      sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
                      sin(x) + cos(x) with extra accuracy */
                    float x0 = 0xe.d4108p+124f;
-- 
2.27.0

> From: Andreas Schwab <schwab@linux-m68k.org>
> Date: Tue, 28 Jul 2020 11:19:28 +0200
> 
> On Jul 28 2020, Paul Zimmermann wrote:
> 
> > +                  /* we subtract (exactly) a value x0 such that cos(x0)+sin(x0)
> > +                     is very near from 0, and use the identity
> 
> Did you mean "near to"?
> 
> Andreas.
> 
> -- 
> Andreas Schwab, schwab@linux-m68k.org
> GPG Key fingerprint = 7578 EB47 D4E5 4D69 2510  2552 DF73 E780 A9DA AEC1
> "And now for something completely different."

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny

[Adhemerval Zanella via Libc-alpha]

On 28/07/2020 07:50, Paul Zimmermann wrote:
>        Dear Andreas,
> 
> yes thanks. Sorry my english is not perfect.
> 
> Paul

Could you send v2 patch with all the fixes indicated by Joseph and Andreas
(this change from a change format is confusing)?  Also please fix the
indentation issue and the open brackets on next line. 

I also think this fix should also add an entry on math/auto-libm-test-out-y0
that exercises this code path and with a check if the ULPs file require
some adjustments as well. 

> 
> From 6b731f36b1a5badf4704645d0dda40957cedd0db Mon Sep 17 00:00:00 2001
> From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
> Date: Mon, 27 Jul 2020 19:01:18 +0200
> Subject: [PATCH 1/3] fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x)
>  is tiny
> 
> ---
>  sysdeps/ieee754/flt-32/e_j0f.c | 16 ++++++++++++++++
>  1 file changed, 16 insertions(+)
> 
> diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
> index c89b9f2688..f85d8a59e0 100644
> --- a/sysdeps/ieee754/flt-32/e_j0f.c
> +++ b/sysdeps/ieee754/flt-32/e_j0f.c
> @@ -56,6 +56,22 @@ __ieee754_j0f(float x)
>  		    if ((s*c)<zero) cc = z/ss;
>  		    else	    ss = z/cc;
>  		}
> +                else {
> +                  /* we subtract (exactly) a value x0 such that cos(x0)+sin(x0)
> +                     is very near from 0, and use the identity
> +                     sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
> +                     sin(x) + cos(x) with extra accuracy */
> +                   float x0 = 3.153646966e+38f;
> +                   float y = x - x0; /* exact */
> +                   /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
> +                   z = __sinf (y);
> +                   float eps = 8.17583368e-8f;
> +                   /* cos(x0) ~ -sin(x0) + eps */
> +                   z += eps * __cosf (x);
> +                   /* now z ~ (sin(x)-cos(x))*cos(x0) */
> +                   float cosx0 = -0.707106740f;
> +                   cc = z / cosx0;
> +                }
>  	/*
>  	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
>  	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
> 

fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny (v2)

[Paul Zimmermann]

       Dear Adhemerval,

thank you for your review. Here is a v2 with all fixes. I did fix the
indentation issue, but did not see the "open brackets on next line" issue:
in other places of this file, we have "else {" with the open brackets on the
same line. I also did add an entry in math/auto-libm-test-in that corresponds
to the larger error for the new code path. No adjustment is needed however,
since the new code is more accurate.

Best regards,
Paul

From a4fe8c3e6c172c7eea198f3225efb05b6afe0f65 Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Wed, 29 Jul 2020 08:59:12 +0200
Subject: [PATCH] fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is
 tiny (v2)

---
 math/auto-libm-test-in         |  2 ++
 sysdeps/ieee754/flt-32/e_j0f.c | 16 ++++++++++++++++
 2 files changed, 18 insertions(+)

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 4414e54d93..5d488a8711 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -5748,6 +5748,8 @@ j0 0x1p16382
 j0 0x1p16383
 # the next value generates larger error bounds on x86_64 (binary32)
 j0 0x2.602774p+0 xfail-rounding:ibm128-libgcc
+# the next value exercises the flt-32 code path for x >= 2^127
+j0 0x8.2f4ecp+124
 
 j1 -1.0
 j1 0.0
diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index c89b9f2688..8540d00b7b 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -56,6 +56,22 @@ __ieee754_j0f(float x)
 		    if ((s*c)<zero) cc = z/ss;
 		    else	    ss = z/cc;
 		}
+                else {
+                  /* we subtract (exactly) a value x0 such that cos(x0)+sin(x0)
+                     is very near to 0, and use the identity
+                     sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
+                     sin(x) + cos(x) with extra accuracy */
+                  float x0 = 0xe.d4108p+124f;
+                  float y = x - x0; /* exact */
+                  /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
+                  z = __sinf (y);
+                  float eps = 0x1.5f263ep-24f;
+                  /* cos(x0) ~ -sin(x0) + eps */
+                  z += eps * __cosf (x);
+                  /* now z ~ (sin(x)-cos(x))*cos(x0) */
+                  float cosx0 = -0xb.504f3p-4f;
+                  cc = z / cosx0;
+                }
 	/*
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
-- 
2.27.0

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny (v2)

[Adhemerval Zanella via Libc-alpha]

On 29/07/2020 04:03, Paul Zimmermann wrote:
>        Dear Adhemerval,
> 
> thank you for your review. Here is a v2 with all fixes. I did fix the
> indentation issue, but did not see the "open brackets on next line" issue:
> in other places of this file, we have "else {" with the open brackets on the
> same line. I also did add an entry in math/auto-libm-test-in that corresponds
> to the larger error for the new code path. No adjustment is needed however,
> since the new code is more accurate.

Some code were imported from other projects and does not follow the glibc
code guidelines.  I am not sure which is the correct strategy to follow
in such cases, but even this patch does not really follow the file 
indentation (for instance, 'else' are added after the close bracket, so
no new line in this case).

In any case, following the already set file style should be fine.

> 
> Best regards,
> Paul
> 
> From a4fe8c3e6c172c7eea198f3225efb05b6afe0f65 Mon Sep 17 00:00:00 2001
> From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
> Date: Wed, 29 Jul 2020 08:59:12 +0200
> Subject: [PATCH] fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is
>  tiny (v2)
> 
> ---
>  math/auto-libm-test-in         |  2 ++
>  sysdeps/ieee754/flt-32/e_j0f.c | 16 ++++++++++++++++
>  2 files changed, 18 insertions(+)
> 
> diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
> index 4414e54d93..5d488a8711 100644
> --- a/math/auto-libm-test-in
> +++ b/math/auto-libm-test-in
> @@ -5748,6 +5748,8 @@ j0 0x1p16382
>  j0 0x1p16383
>  # the next value generates larger error bounds on x86_64 (binary32)
>  j0 0x2.602774p+0 xfail-rounding:ibm128-libgcc
> +# the next value exercises the flt-32 code path for x >= 2^127
> +j0 0x8.2f4ecp+124
>  
>  j1 -1.0
>  j1 0.0

Does it require any ULP adjustment (at least on the architecture you
tested it)?

> diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
> index c89b9f2688..8540d00b7b 100644
> --- a/sysdeps/ieee754/flt-32/e_j0f.c
> +++ b/sysdeps/ieee754/flt-32/e_j0f.c
> @@ -56,6 +56,22 @@ __ieee754_j0f(float x)
>  		    if ((s*c)<zero) cc = z/ss;
>  		    else	    ss = z/cc;
>  		}
> +                else {

No new line for the else ( '} else {' ).

> +                  /* we subtract (exactly) a value x0 such that cos(x0)+sin(x0)

s/we/We.

> +                     is very near to 0, and use the identity
> +                     sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
> +                     sin(x) + cos(x) with extra accuracy */

Period and double space prior comment close ('[...] accuracy.  */').

> +                  float x0 = 0xe.d4108p+124f;
> +                  float y = x - x0; /* exact */
> +                  /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
> +                  z = __sinf (y);
> +                  float eps = 0x1.5f263ep-24f;
> +                  /* cos(x0) ~ -sin(x0) + eps */
> +                  z += eps * __cosf (x);
> +                  /* now z ~ (sin(x)-cos(x))*cos(x0) */
> +                  float cosx0 = -0xb.504f3p-4f;
> +                  cc = z / cosx0;
> +                }
>  	/*
>  	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
>  	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
> 

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny (v3)

[Paul Zimmermann]

       Dear Adhemerval,

thank you again for your useful comments.

> Does it require any ULP adjustment (at least on the architecture you
> tested it)?

I reset to 0 the j0/float values in sysdeps/x86_64/fpu/libm-test-ulps,
then ran "make regen-ulps", and got:

$ diff ../sysdeps/x86_64/fpu/libm-test-ulps /localdisk/zimmerma/glibc/build/math/NewUlps
1315c1315
< float: 0
---
> float: 8
1321c1321
< float: 0
---
> float: 4
1327c1327
< float: 0
---
> float: 5
1333c1333
< float: 0
---
> float: 5

which were the original values except for j0_towardzero (it was 6).
However the same applies to "master", so this is independent from my change.
Anyway, I decreased the value from 6 to 5.

Below is a v3.

Paul

From 34d336cf2e8495bba21254157be13e02f52d934a Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Thu, 30 Jul 2020 09:16:36 +0200
Subject: [PATCH] fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is
 tiny (v3)

---
 math/auto-libm-test-in            |  2 ++
 sysdeps/ieee754/flt-32/e_j0f.c    | 17 ++++++++++++++++-
 sysdeps/x86_64/fpu/libm-test-ulps |  2 +-
 3 files changed, 19 insertions(+), 2 deletions(-)

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 4414e54d93..5d488a8711 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -5748,6 +5748,8 @@ j0 0x1p16382
 j0 0x1p16383
 # the next value generates larger error bounds on x86_64 (binary32)
 j0 0x2.602774p+0 xfail-rounding:ibm128-libgcc
+# the next value exercises the flt-32 code path for x >= 2^127
+j0 0x8.2f4ecp+124
 
 j1 -1.0
 j1 0.0
diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index c89b9f2688..91e8de8fe3 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -55,7 +55,22 @@ __ieee754_j0f(float x)
 		    z = -__cosf(x+x);
 		    if ((s*c)<zero) cc = z/ss;
 		    else	    ss = z/cc;
-		}
+		} else {
+                  /* We subtract (exactly) a value x0 such that cos(x0)+sin(x0)
+                     is very near to 0, and use the identity
+                     sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
+                     sin(x) + cos(x) with extra accuracy. */
+                  float x0 = 0xe.d4108p+124f;
+                  float y = x - x0; /* exact */
+                  /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
+                  z = __sinf (y);
+                  float eps = 0x1.5f263ep-24f;
+                  /* cos(x0) ~ -sin(x0) + eps */
+                  z += eps * __cosf (x);
+                  /* now z ~ (sin(x)-cos(x))*cos(x0) */
+                  float cosx0 = -0xb.504f3p-4f;
+                  cc = z / cosx0;
+                }
 	/*
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index d71d86d961..2561243fb7 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1324,7 +1324,7 @@ ldouble: 4
 
 Function: "j0_towardzero":
 double: 5
-float: 6
+float: 5
 float128: 2
 ldouble: 5
 
-- 
2.27.0

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny (v3)

[Florian Weimer via Libc-alpha]

* Paul Zimmermann:

>        Dear Adhemerval,
>
> thank you again for your useful comments.
>
>> Does it require any ULP adjustment (at least on the architecture you
>> tested it)?
>
> I reset to 0 the j0/float values in sysdeps/x86_64/fpu/libm-test-ulps,
> then ran "make regen-ulps", and got:
>
> $ diff ../sysdeps/x86_64/fpu/libm-test-ulps /localdisk/zimmerma/glibc/build/math/NewUlps
> 1315c1315
> < float: 0
> ---
>> float: 8
> 1321c1321
> < float: 0
> ---
>> float: 4
> 1327c1327
> < float: 0
> ---
>> float: 5
> 1333c1333
> < float: 0
> ---
>> float: 5
>
> which were the original values except for j0_towardzero (it was 6).
> However the same applies to "master", so this is independent from my change.
> Anyway, I decreased the value from 6 to 5.

I suggest to leave it at 6, other CPU variants may still need the 6
there.

Thanks,
Florian

fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny (v4)

[Paul Zimmermann]

       Dear Florian,

> I suggest to leave it at 6, other CPU variants may still need the 6
> there.

here is a new version.

Paul

From 971c832c2087e9463951d1b07b48d4d9a998e8c0 Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Mon, 3 Aug 2020 15:16:39 +0200
Subject: [PATCH] fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is
 tiny (v4)

---
 math/auto-libm-test-in         |  2 ++
 sysdeps/ieee754/flt-32/e_j0f.c | 17 ++++++++++++++++-
 2 files changed, 18 insertions(+), 1 deletion(-)

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 4414e54d93..5d488a8711 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -5748,6 +5748,8 @@ j0 0x1p16382
 j0 0x1p16383
 # the next value generates larger error bounds on x86_64 (binary32)
 j0 0x2.602774p+0 xfail-rounding:ibm128-libgcc
+# the next value exercises the flt-32 code path for x >= 2^127
+j0 0x8.2f4ecp+124
 
 j1 -1.0
 j1 0.0
diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index c89b9f2688..91e8de8fe3 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -55,7 +55,22 @@ __ieee754_j0f(float x)
 		    z = -__cosf(x+x);
 		    if ((s*c)<zero) cc = z/ss;
 		    else	    ss = z/cc;
-		}
+		} else {
+                  /* We subtract (exactly) a value x0 such that cos(x0)+sin(x0)
+                     is very near to 0, and use the identity
+                     sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
+                     sin(x) + cos(x) with extra accuracy. */
+                  float x0 = 0xe.d4108p+124f;
+                  float y = x - x0; /* exact */
+                  /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
+                  z = __sinf (y);
+                  float eps = 0x1.5f263ep-24f;
+                  /* cos(x0) ~ -sin(x0) + eps */
+                  z += eps * __cosf (x);
+                  /* now z ~ (sin(x)-cos(x))*cos(x0) */
+                  float cosx0 = -0xb.504f3p-4f;
+                  cc = z / cosx0;
+                }
 	/*
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
-- 
2.27.0

fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny (v4)

[Paul Zimmermann]

ping: https://sourceware.org/pipermail/libc-alpha/2020-August/116806.html

Have a nice week-end,
Paul

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny (v4)

[Adhemerval Zanella via Libc-alpha]

On 07/08/2020 08:25, Paul Zimmermann wrote:
> ping: https://sourceware.org/pipermail/libc-alpha/2020-August/116806.html
> 
> Have a nice week-end,
> Paul
> 

LGTM and I pushed it upstream with some minor fixes (I fixed the indentation
to match the rest of the file, added double space prior comment endding and
adjust one line that surpassed the 78 column limit).

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny (v4)

[Joseph Myers]

On Fri, 7 Aug 2020, Adhemerval Zanella via Libc-alpha wrote:

> On 07/08/2020 08:25, Paul Zimmermann wrote:
> > ping: https://sourceware.org/pipermail/libc-alpha/2020-August/116806.html
> > 
> > Have a nice week-end,
> > Paul
> > 
> 
> LGTM and I pushed it upstream with some minor fixes (I fixed the indentation
> to match the rest of the file, added double space prior comment endding and
> adjust one line that surpassed the 78 column limit).

That commit is missing the regeneration of auto-libm-test-out-j0.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny (v4)

[Adhemerval Zanella via Libc-alpha]

On 07/08/2020 17:30, Joseph Myers wrote:
> On Fri, 7 Aug 2020, Adhemerval Zanella via Libc-alpha wrote:
> 
>> On 07/08/2020 08:25, Paul Zimmermann wrote:
>>> ping: https://sourceware.org/pipermail/libc-alpha/2020-August/116806.html
>>>
>>> Have a nice week-end,
>>> Paul
>>>
>>
>> LGTM and I pushed it upstream with some minor fixes (I fixed the indentation
>> to match the rest of the file, added double space prior comment endding and
>> adjust one line that surpassed the 78 column limit).
> 
> That commit is missing the regeneration of auto-libm-test-out-j0.
> 

Done, I also updated both x86_64 and i686 ulp files that resulted from the
the new test.

Re: fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny (v4)

[Paul Zimmermann]

       Dear Joseph, Dear Adhemerval,

sorry for that, I'll try to remember next time to regenerate the
auto-libm-test-out-xxx file and to check the ulp files.

Paul

> From: Adhemerval Zanella <adhemerval.zanella@linaro.org>
> Date: Sat, 8 Aug 2020 16:53:45 -0300
> 
> On 07/08/2020 17:30, Joseph Myers wrote:
> > On Fri, 7 Aug 2020, Adhemerval Zanella via Libc-alpha wrote:
> > 
> >> On 07/08/2020 08:25, Paul Zimmermann wrote:
> >>> ping: https://sourceware.org/pipermail/libc-alpha/2020-August/116806.html
> >>>
> >>> Have a nice week-end,
> >>> Paul
> >>>
> >>
> >> LGTM and I pushed it upstream with some minor fixes (I fixed the indentation
> >> to match the rest of the file, added double space prior comment endding and
> >> adjust one line that surpassed the 78 column limit).
> > 
> > That commit is missing the regeneration of auto-libm-test-out-j0.
> > 
> 
> Done, I also updated both x86_64 and i686 ulp files that resulted from the
> the new test.
>