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

[Adhemerval Zanella]

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