mp.asin(z).imag and fp.asin(z).imag use inconsistent sign choices when z.real > 1 and z.imag==0 about mpmath HOT 23 CLOSED

Considering all data points, I think the expected result of asin(2) is 1.570796+1.316958i, especially because this is the result that the evaluation of the mathematical definition of asin(z) = -i log(i z + sqrt(1 - z*z)) would give with mixed-mode arithmetic.

I would interpret the identity exactly rather than with mixed mode arithmetic:

asin(2) = - i log(i * 2 + sqrt(1 - 2^2))
    = - i log (2 i + sqrt(3) i )
    = - i (log(| 2 + sqrt(3) |) + i pi / 2)
    = pi/2 - i log(2 + sqrt(3))
    ~= 1.57 - 1.31 i

With SymPy that is:

In [1]: asin(2)
Out[1]: asin(2)

In [2]: asin(2).rewrite(log)
Out[2]: 
   ⎛              ⅈ⋅π⎞
-ⅈ⋅⎜log(√3 + 2) + ───⎟
   ⎝               2 ⎠

In [3]: asin(2).rewrite(log).expand()
Out[3]: 
π                
─ - ⅈ⋅log(√3 + 2)
2                

In [4]: asin(2).rewrite(log).expand().n()
Out[4]: 1.5707963267949 - 1.31695789692482⋅ⅈ

In [5]: asin(z).rewrite(log) # holds for all z
Out[5]: 
      ⎛         ________⎞
      ⎜        ╱      2 ⎟
-ⅈ⋅log⎝ⅈ⋅z + ╲╱  1 - z  ⎠

I don't think that signed zeros are really relevant unless you want to think about what to do with 2 - 0.0 i. In general signed zeros are (IMO) a mistake or at least to be implemented properly you need three cases: -0, 0 and +0. Without all three cases you cannot use signed zeros to resolve branch cuts in general.

The decision about which way to resolve the branch cut is always arbitrary but there is at least a principle here that a useful identity holds with a consistent choice of branch for each function.

Besides consistency what would be an a priori reason to prefer the alternative?

If consistency is the goal then I would put mpmath/arb consistency as most important.

from mpmath.

I wonder if that's something, that must be documented in https://mpmath.readthedocs.io/en/latest/contexts.html#contexts

from mpmath.

This coming from:

acos() case seems wrong too:

>>> import mpmath
>>> mpmath.mp.acos(2)
mpc(real='0.0', imag='1.3169578969248168')
>>> mpmath.fp.acos(2)
-1.3169578969248166j

from mpmath.

#788

from mpmath.

It seems, mpmath behaves here as Mathematica:

$ wolframscript 
Wolfram Language 14.0.0 Engine for Linux x86 (64-bit)
Copyright 1988-2023 Wolfram Research, Inc.

In[1]:= ArcSin[2]//N

Out[1]= 1.5708 - 1.31696 I

Looks weird, that different contexts use different definitions of elementary functions on the complex plane. So, I would prefer to change this. On another hand, it will be a compatibility break.

from mpmath.

This is what arb gives:

In [9]: flint.acb(2).asin()
Out[9]: [1.57079632679490 +/- 3.56e-15] + [-1.31695789692482 +/- 3.71e-15]j

In [10]: flint.acb(2+.001j).asin()
Out[10]: [1.57021897662193 +/- 2.11e-15] + [1.31695808937485 +/- 4.93e-15]j

In [11]: flint.acb(2-.001j).asin()
Out[11]: [1.57021897662193 +/- 2.11e-15] + [-1.31695808937485 +/- 3.07e-15]j

Both arb and mpmath document the identity

asin(z) = -i log(i z + sqrt(1 - z^2))

The mpmath docs explicitly say that the principal branches of log and sqrt are intended in this identity:

In [15]: -1j*cmath.log(1j*2 + cmath.sqrt(1-2**2))
Out[15]: (1.5707963267948966-1.3169578969248166j)

In [16]: cmath.asin(2)
Out[16]: (1.5707963267948966+1.3169578969248166j)

from mpmath.

More data points:

Numpy, PyTorch, Boost, cmath: asin(2) -> 1.5608 + 1.3170 j
Mathematica, Maple, sympy, matlab, jax: asin(2) -> 1.5608 - 1.3170 j

from mpmath.

FWIW, one more data point:

CLN: asin(2.0) -> 1.5707964-1.3169578969248166d0i

This is consistent with its documentation for cln::asin, which says

cl_N asin (const cl_N& z)

Returns arcsin(z). This is defined as arcsin(z) = log(iz+sqrt(1-z^2))/i and satisfies arcsin(-z) = -arcsin(z). The range of the result is the strip in the complex domain -pi/2 <= realpart(arcsin(z)) <= pi/2, excluding the numbers with realpart = -pi/2 and imagpart < 0 and the numbers with realpart = pi/2 and imagpart > 0.

from mpmath.

https://dl.acm.org/doi/10.1145/275323.275324 defines

from mpmath.

In [15]: -1jcmath.log(1j2 + cmath.sqrt(1-2**2))
Out[15]: (1.5707963267948966-1.3169578969248166j)

Hmm, I'm not sure if this is not a CPython bug (see below C version)...

Consider:

>>> cmath.asin(complex(2,+0.0001))
(1.5707385917680738+1.3169578988493174j)
>>> -1j*cmath.log(1j*complex(2,0.0001) + cmath.sqrt(1-complex(2,0.0001)**2))
(1.5707385917680738+1.3169578988493182j)
>>> cmath.asin(complex(2,-0.0001))
(1.5707385917680738-1.3169578988493174j)
>>> -1j*cmath.log(1j*complex(2,-0.0001) + cmath.sqrt(1-complex(2,-0.0001)**2))
(1.5707385917680738-1.3169578988493176j)

Also,

$ cat a.c 
#include <complex.h>
#include <stdio.h>
void
main()
{
    double complex z;
    z = CMPLX(2, 0.0); z = casin(z);
    printf("casin(+2+0i) = %f%+fi\n", creal(z), cimag(z));
    z = CMPLX(2, 0.0); z = -I*clog(I*z + csqrt(1 - z*z));
    printf("casin(+2+0i) = %f%+fi\n", creal(z), cimag(z));
    z = CMPLX(2, -0.0); z = casin(z);
    printf("casin(+2-0i) = %f%+fi\n", creal(z), cimag(z));
    z = CMPLX(2, -0.0); z = -I*clog(I*z + csqrt(1 - z*z));
    printf("casin(+2-0i) = %f%+fi\n", creal(z), cimag(z));
}
$ cc a.c -lm && ./a.out
casin(+2+0i) = 1.570796+1.316958i
casin(+2+0i) = 1.570796+1.316958i
casin(+2-0i) = 1.570796-1.316958i
casin(+2-0i) = 1.570796-1.316958i

from mpmath.

I am not sure what do you consider as a bug here... I think it is just unfortunate that there is a split between asin implementations and apparently the split is defined by either fixing a formula for asin function or choosing a branch that corresponds to a PV of arcsin. Of course, from the mathematical point of view, both selections give a correct answer in the sense that sin(asin(2)) == 2 holds in both cases.

On the other hand, I'd like to consider mpmath as a reference library to elementary functions and use it for testing the correctness of functions from other numerical libraries. Since different libraries may implement mathematically equivalent but numerically different algorithms, as evident from this asin case, one approach of resolving this dilemma is to introduce configuration or environment variable that defines the behavior of asin at the positive real line. Say, have:

$ MPMATH_ASIN_USE_PV_BRANCH=0 python -c 'import mpmath; print(mpmath.asin(2))'
(1.5707963267949 - 1.31695789692482j)
$ MPMATH_ASIN_USE_PV_BRANCH=1 python -c 'import mpmath; print(mpmath.asin(2))'
(1.5707963267949 + 1.31695789692482j)

What do you think?

from mpmath.

I am not sure what do you consider as a bug here...

Well, we have an identity asin(z) = -i log(i z + sqrt(1 - z^2)). As you can see, the C version preserves this for z=2±0i, but not the cmath's version. One reason for this is python/cpython#117999.

one approach of resolving this dilemma is to introduce configuration or environment variable that defines the behavior of asin at the positive real line

Does make sense for me, if there will be just one configuration variable. In general, each boolean variable will double test cases...

https://dl.acm.org/doi/10.1145/275323.275324 defines

Hmm, it corresponds to the mpmath docs and to the C stdlib (z=2+0i).

BTW, note that so far the mp context has no signed zero (technically, it's here, but not supported even in basic arithmetic or I/O). The fp context has it.

from mpmath.

Yet one point (gmpy2):

>>> import gmpy2
>>> c = gmpy2.get_context()
>>> c.allow_complex=True
>>> gmpy2.asin(2)
1.5707963267948966+1.3169578969248168j

In short, it seems all libraries, that support IEEE-style arithmetic - use above convention. Not sure about benefits of the Mathematica approach (except for compatibility). @fredrik-johansson?

BTW, Oscar's example -1j*cmath.log(1j*2 + cmath.sqrt(1-2**2)) suffers from float -> complex(float, 0.0) type-cast of the CPython complex arithmetic (yet another instance of this issue). It can be corrected:

>>> -1j*cmath.log(1j*z + cmath.sqrt(complex(1-(z*z).real, -(z*z).imag)))
(1.5707963267948966+1.3169578969248164j)

The C stdlib's complex numbers have special rules for mixed-mode arithmetic (real+complex and imaginary+complex; however most implementations support only real+complex specialization). That's why C version does work without workarounds:

$ cat a.c
#include <complex.h>
#include <stdio.h>
void
main()
{
    double complex z;
    z = CMPLX(2, 0.0); z = casin(z);
    printf("%f%+fi\n", creal(z), cimag(z));
    z = CMPLX(2, 0.0); z = -I*clog(I*z + csqrt(1 - z*z));
    printf("%f%+fi\n", creal(z), cimag(z));
}
$ cc a.c -lm && ./a.out 
1.570796+1.316958i
1.570796+1.316958i

"Generally, mixed-mode arithmetic combining real and complex variables should be performed directly, not by first coercing the real to complex, lest the sign of zero be rendered uninformative;the same goes for combinations of pure imaginary quantities with complex variables." (c) Kahan, W: Branch cuts for complex elementary functions

from mpmath.

Considering all data points, I think the expected result of asin(2) is 1.570796+1.316958i, especially because this is the result that the evaluation of the mathematical definition of asin(z) = -i log(i z + sqrt(1 - z*z)) would give with mixed-mode arithmetic.

When a software returns 1.570796-1.316958i (notice the negative imaginary part) then it is very likely that it (carelessly) uses the mathematical definition of asin using complex arithmetic that makes no difference between +0 and -0.

BC issue with respect to Mathematica/Maple and implications to sympy is an annoying problem but I think the correctness over weights BC, not to mention the fact that a majority of math libraries already implement the correct behavior w/ this issue.

What do you think?

from mpmath.

using complex arithmetic that makes no difference between +0 and -0.

That was my guess, but it seems Maple supports signed zeros.

the fact that a majority of math libraries already implement the correct behavior w/ this issue.

I doubt the Mathematica behaviour is "incorrect". Yes, it makes more sense in floating point arithmetic without support for signed zeros; see above Kahan reference and CCC (Counter-Clockwise Continuity) convention. To illustrate one for re<0:

In[9]:= ArcSin[-2.+I*0.0001]
Out[9]= -1.57074 + 1.31696 I
In[10]:= ArcSin[-2.]
Out[10]= -1.5708 + 1.31696 I

BTW, yet another data point, Maxima:

(%i5) asin(2.0);
(%o5)              1.570796326794897 - 1.316957896924817 %i

PS: I would worry about flint's convention also as it's a possible backend for mpmath, see #733.

from mpmath.

Considering all data points, I think the expected result of asin(2) is 1.570796+1.316958i, especially because this is the result that the evaluation of the mathematical definition of asin(z) = -i log(i z + sqrt(1 - z*z)) would give with mixed-mode arithmetic.

I would interpret the identity exactly rather than with mixed mode arithmetic:

asin(2) = - i log(i * 2 + sqrt(1 - 2^2))

Considering that sqrt(...) denotes +-sqrt(...), you could write this also as

asin(2) = - i log(i * 2 - sqrt(1 - 2^2))

that will lead to the result with positive imaginary part. See, for instance, https://proofwiki.org/wiki/Complex_Inverse_Sine/Examples/2

It looks like symbolic and numeric libraries consistently pick up different branches when evaluating arcsin at branch cuts and since mpmath is relevant for both universes then letting the user to choose the branch cut as suggested in #786 (comment) still makes sense.

from mpmath.

Considering that sqrt(...) denotes +-sqrt(...)

The point is that the formula holds with a consistent choice of branch cuts for each function so sqrt means +sqrt here:

from mpmath.

Interestingly, the paper On quality of implementation of Fortran 2008
complex intrinsic functions on branch cuts finds that

arcsin(a+i0) == pi/2 + b i
b = log(sqrt(a**2 - 1) + a) >= 0
a >= 1

(see A.3, p 23).

from mpmath.

Ok, @oscarbenjamin convinced me that the current convention does make sense in floating point arithmetic without signed zero. And, in turn, that it doesn't make sense with signed zero.

Given mp's mpf type has no signed zero ("officially" (tm)), we should keep current convention, unless mpf got support for signed zero per #167. I think we could switch convention in the experimental pr #768.

I'm going to close this issue in few days, unless someone introduce new arguments on the table.

from mpmath.

Just to confirm: when mpmath.mp and mpmath.fp functions return different results, it can be considered acceptable?

from mpmath.

Yes, for branch cuts: mp and fp are different contexts. In the fp we can't say "we are on the branch cut", but "we are on the left side or on the right side"/"above or below". Both contexts should be consistent, however, for arbitrary small (but nonzero!) deviations from the branch cut. This is the case for asin as well:

>>> fp.asin(2+1e-10j)
(1.5707963267371616+1.3169578969248166j)
>>> mp.asin(2+1e-10j)
mpc(real='1.5707963267371616', imag='1.3169578969248168')
>>> fp.asin(2-1e-10j)
(1.5707963267371616-1.3169578969248166j)
>>> mp.asin(2-1e-10j)
mpc(real='1.5707963267371616', imag='-1.3169578969248168')

On another hand, I think this will require that functions should be implemented slightly differently for fp and mp contexts. When fp context fallback to generic code (which is the case for most special functions) - we will have invalid results... An example:

>>> complex(mp.besselk(1, complex(-1, +1e-10)))
(-0.6019072304174309-1.775499689314474j)
>>> complex(fp.besselk(1, complex(-1, +1e-10)))
(-0.6019200766099881-1.7753783061922508j)
>>> complex(fp.besselk(1, complex(-1, +0.)))
(-0.6019200763910812-1.775378306087049j)
>>> complex(mp.besselk(1, complex(-1, -1e-10)))
(-0.6019072304174309+1.775499689314474j)
>>> complex(fp.besselk(1, complex(-1, -1e-10)))
(-0.6019200766099881+1.7753783061922508j)
>>> complex(fp.besselk(1, complex(-1, -0.)))  # oops
(0.6019200763910814+1.775378306087049j)

from mpmath.

#795

from mpmath.

#797

from mpmath.