Re: [musl] casinh function accuracy problem

[Rich Felker <dalias@libc.org>]

On Wed, Feb 12, 2020 at 03:07:17PM +1100, Damian McGuckin wrote:
> On Tue, 11 Feb 2020, Rich Felker wrote:
> 
> >My minimal criterion for large-scale improvements of src/complex would
> >be fixing any remaining cases where inf/nan behavior is badly wrong or
> >there's catastrophic error (>2^52 ulp, or even just >2^20 ulp or so).
> >Beyond that, I think "reducing ulp error" would be nice but hard to
> >quantify and make a goal without having an idea how bad it is now, not
> >to mention without having rigorous error bounds on the real math
> >library functions.
> 
> I think INF/NaN behaviour at the fundamental level is flawed.
> 
> This initialization:
> 
> 	double complex x = 1.0e+200 + INFINITY * I;
> 
> on every compiler I try, yields an 'x' of
> 
> 	NaN + INFINITY i

This is expected unless the implementation supports pure-imaginary
types and a _Imaginary_I (as opposed to _Complex_I). Note that Annex G
requires this, but GCC treats it as optional and does not support it
despite claiming Annex G conformance.

In this particular kind of case, the right solution is to use the
CMLPX macro instead of the * operator with I:

	double complex x = CMPLX(1.0e+200, INFINITY);

> Whereas if I compute
> 
> 	double complex a = 2.0 + 1.0e200 * I;
> 	double complex b = 1.0e200 + 1.0 * I;
> 	double x = a * b;
> 
> then 'x' prints correctly as
> 
> 	1.0e+200 + INF * I;

Indeed the compiler is required to handle these overflow cases unless
#pragma STDC CX_LIMITED_RANGE ON is in effect. I don't know if GCC
honors the pragma but it has -fcx-limited-range to set it globally.
Note that without CX_LIMITED_RANGE, complex math performance is
atrociously bad due to extensive overflow checking.

Rich