Re: [musl] casinh function accuracy problem

[Rich Felker <dalias@libc.org>]

On Tue, Feb 11, 2020 at 08:30:59PM +0100, Szabolcs Nagy wrote:
> * mg1633068 <songyue@smail.nju.edu.cn> [2020-02-11 23:51:51 +0800]:
> > Hi everyone,
> > I'm writing testcases for complex math function. Considering the following simple code
> > 
> > #include <stdio.h>
> > #include <complex.h>
> > int main(int argc, char *argv[])
> > {
> >     double complex d = 3.0+6.6*I;
> >     double complex ret = casinh(d);
> >     printf("casinh(3.0+6.6*I)=%.15f+%.15f*I\n", creal(ret), cimag(ret));
> > 
> >     return 0;
> > }
> > 
> > With musl libc, the result is:
> >     casinh(3.0+6.6*I)=2.671002221994648+1.140551372972568*I
> > but with glibc, the result is:
> >     casinh(3.0+6.6*I)=2.671002221994652+1.140551372972565*I
> > 
> > We can see that musl is less accurate. I'm trying to solve this problem.
> > With little knowledge of numerical computing, any comment is appreciated!
> 
> do you mean you are trying to fix the code in musl?
> 
> that's welcome, but i think it will be hard without
> numerical computing knowledge.
> 
> several complex math functions in musl are not
> correct (implemented in a very naive way), but
> fixing them is significant effort.
> 
> in this particular case the 8 and 12 ulp errors
> on the real and imaginary parts are still
> considered small errors (glibc has 1 and 2 ulp
> errors compared to the correctly rounded result).

Indeed, doing much better than musl's complex functions do now looks
like it would require a good deal of numerical analysis. The current
ones are based pretty much entirely on identities that are true on the
complex numbers, but that introduce additional error at each step when
used on floating point representations. I'm kinda surprised they do as
well as they do now.

On targets with a long double type having more precision than double,
you might be able to cheat and call the long double versions instead,
then drop to double in the result, to get a few more places of
precision. But that won't help on targets where ld==double.

High-quality complex math functions are a long-term wishlist item for
musl but nobody has stepped up to do them and I don't really feel like
doing it, at least not over other improvements I could be working on.
This might be an area well-served by sponsored enhancement if there's
a user who needs them improved with resources to pay someone to do it.

Rich