[<prev] [next>] [<thread-prev] [day] [month] [year] [list]
```Date: Thu, 17 May 2018 18:11:53 +0200
From: Szabolcs Nagy <nsz@...t70.net>
To: musl@...ts.openwall.com
Subject: Re: What should be the result of CACOSH(F)(CCOSH(F)(-2 + 1j))?

* Paolo Mantegazza <paolo.mantegazza@...imi.it> [2018-05-17 14:32:40 +0200]:
> Hi,
>
> calling: j = sqrt(-1),
>
> MUSL        answer is: -0.200000e+1 + 0.100000e+1j;
>
> GLIBC       answer is: +0.20000e+1  -  0.100000e+1j;
>
> SCILAB     answer is: +0.20000e+1  -  0.100000e+1j;
>
> MATLAB   answer is: +0.20000e+1  -  0.100000e+1j;
>
> UCLIB-NG answer is: +0.20000e+1  -  0.100000e+1j;
>
> Math is not democracy so maybe MUSL's answer is the right one. In fact, with
> infinite precision at least, one should expect that, by applying the inverse
> of a function to a function, the result should be the used function
> argument.
>
> So, does it either show a missed correct principal value or that MUSL is the
> smartest one?
>

musl does not try very hard to provide correct complex/
so it's easily possible that the answer is wrong.

we can try to get these edge cases right, but there are
lot of large ulp error cases and intermediate value
overflows cases that are not trivial to fix.
(that's why we use the simplest expression in most cases
which is wrong but in an obvious way)

> In any case, following:
> http://mathworld.wolfram.com/InverseHyperbolicCosine.html, a way to have
> MUSL match GLIBC:SCILAB:MATLAB:UCLIB-NG is to change the two code lines in
> MUSL ./src/complex/cacosh.c
>
>      z = cacos(z);
>      return CMPLX(-cimag(z), creal(z));  // AKA j*cacos(z)
>
> into
>
>     return = clog(z + csqrt(z + 1) * csqrt(z - 1));  // AKA a definition of
> cacosh
>
> As a further info, NEWLIB cacosh.c (not tested here) recites:
>
> #if 0 /* does not give the principal value */
>         w = I * cacos(z);
> #else
>         w = clog(z + csqrt(z + 1) * csqrt(z - 1));
> #endif
>         return w;
>
> it should be remarked that, to provide the correct principal value using
> just cacos(z), the above mentioned link addresses the need of testing the
> cacosh argument in order to appropriately use either  j*cacos(z) or
> -j*cacos(z). It is therefore likely that the fix to be chosen, if any,
> should be based on computational efficiency.
>
> Once more math is not democracy, so the answer must be left to the math
> savviest.
>
> Regards, Paolo Mantegazza.
```