Re: [PATCH] Fix spurious underflow exceptions for Bessel functions for double(from glibc bug 14155)

Corinna Vinschen via Newlib <[email protected]> Thu, 26 Mar 2020 12:27:17 +0100
Newsgroups	gmane.comp.lib.newlib
Message-ID	<[email protected]>
On Mar 25 11:18, Keith Packard via Newlib wrote:
> From: "Joseph S. Myers" <[email protected]>
> 
> 	This fix comes from glibc, from files which originated from
> 	the same place as the newlib files. Those files in glibc carry
> 	the same license as the newlib files.
> 
> Bug 14155 is spurious underflow exceptions from Bessel functions for
> large arguments.  (The correct results for large x are roughly
> constant * sin or cos (x + constant) / sqrt (x), so no underflow
> exceptions should occur based on the final result.)
> 
> There are various places underflows may occur in the intermediate
> calculations that cause the failures listed in that bug.  This patch
> fixes problems for the double version where underflows occur in
> calculating the intermediate functions P and Q (in particular, x**-12
> gets computed while calculating Q).  Appropriate approximations are
> used for P and Q for arguments at least 0x1p28 and above to avoid the
> underflows.
> 
> For sufficiently large x - 0x1p129 and above - the code already has a
> cut-off to avoid calculating P and Q at all, which means the
> approximations -0.125 / x and 0.375 / x can't themselves cause
> underflows calculating Q.  This cut-off is heuristically reasonable
> for the point beyond which Q can be neglected (based on expecting
> around 0x1p-64 to be the least absolute value of sin or cos for large
> arguments representable in double).
> 
> The float versions use a cut-off 0x1p17, which is less heuristically
> justifiable but should still only affect values near zeroes of the
> Bessel functions where these implementations are intrinsically
> inaccurate anyway (bugs 14469-14472), and should serve to avoid
> underflows (the float underflow for jn in bug 14155 probably comes
> from the recurrence to compute jn).  ldbl-96 uses 0x1p129, which may
> not really be enough heuristically (0x1p143 or so might be safer - 143
> = 64 + 79, number of mantissa bits plus total number of significant
> bits in representation) but again should avoid underflows and only
> affect values where the code is substantially inaccurate anyway.
> ldbl-128 and ldbl-128ibm share a completely different implementation
> with no such cut-off, which I propose to fix separately.
> 
> Signed-off-by: Keith Packard <[email protected]>
> ---
>  newlib/libm/math/e_j0.c  | 6 ++++--
>  newlib/libm/math/e_j1.c  | 6 ++++--
>  newlib/libm/math/ef_j0.c | 6 +++---
>  newlib/libm/math/ef_j1.c | 4 ++--
>  4 files changed, 13 insertions(+), 9 deletions(-)

Pushed.


Thanks,
Corinna

-- 
Corinna Vinschen
Cygwin Maintainer
Red Hat
signature.asc (application/pgp-signature, 833 B)