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

[Adhemerval Zanella via Libc-alpha <libc-alpha@sourceware.org>]

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