Re: Patch for bz#11589 - jn calculation

[Andreas Schwab <schwab at redhat dot com>]

Andreas Jaeger <aj@suse.de> writes:

> diff --git a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
> index 0eea745..372f942 100644
> --- a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
> +++ b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
> @@ -286,7 +286,16 @@ __ieee754_jnl (n, x)
>  		    }
>  		}
>  	    }
> -	  b = (t * __ieee754_j0l (x) / b);
> +	  /* j0() and j1() suffer enormous loss of precision at and
> +	   * near zero; however, we know that their zero points never
> +	   * coincide, so just choose the one further away from zero.
> +	   */
> +	  z = __ieee754_j0l (x);
> +	  w = __ieee754_j1l (x);
> +	  if (fabsl (z) >= fabsl (w))
> +	    b = (t * z / b);
> +	  else
> +	    b = (t * w / a);
>  	}
>      }
>    if (sgn == 1)

Something's fishy here.  I'm getting huge ulps both on powerpc and s390
for all the new tests.  For example:

testing long double (without inline functions)
Failure: Test: jn (2, 2.4048255576957729) == 0.43175480701968038399746111312430703
Result:
 is:          4.31754807019680381803e-01   0x1.ba1deea029493e99a8fcp-2
 should be:   4.31754807019680383997e-01   0x1.ba1deea029493f3b959ep-2
 difference:  2.19448415868800798091e-18   0x1.43d9447b3c18e0000000p-59
 ulp       :  356075412601880.8750
 max.ulp   :  0.0000

Andreas.

-- 
Andreas Schwab, schwab@redhat.com
GPG Key fingerprint = D4E8 DBE3 3813 BB5D FA84  5EC7 45C6 250E 6F00 984E
"And now for something completely different."