Both cmath and numpy give "incorrect" value of asin(10)

Question

I recently needed to quickly find the arcsine of 10. I decided to use python to calculate it for me:

cmath.asin(10)

Based on experience I had expected a result in the 4th quadrant (positive real(pi/2) and negative imaginary). Surprise... it returned the 1st quadrant result. I tried numpy.arcsin as well... same result. While it is true that the sine of the returned value is 10, I don't believe this is the standard principal value for the arcsine function.

>>> import math
>>> import numpy as np
>>> z=cmath.asin(10)
>>> z
(1.5707963267948966+2.993222846126381j)
>>> cmath.sin(z)
(9.999999999999998+6.092540900222251e-16j)
>>> z2=np.arcsin(10+0j)
>>> np.sin(z2)
(10+6.092540900222253e-16j)

I found this same result using numpy (as seen above).

Is there a python module from which I can expect to get the standard principal values (i.e. follows the principal branch cuts) for complex valued functions? Or, is the notion of standard principal value too fluid at this point to expect consistency?

---

• Pushing through the discussion on the Wikipedia page (https://en.wikipedia.org/wiki/Inverse_trigonometric_functions) section on complex values leads to the 4th quadrant solution.

• This is also true for an article (http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.413.5729&rep=rep1&type=pdf) published by a PSU mathematician on complex inverse trig functions.

• Wolfram Alpha gives the 4th quadrant solution.

Answer 1

The cmath behaviour is somewhat standard, in that it's not limited to just cmath and NumPy: it also matches the behaviour recommended in Annex G of the C standards (at least, C99 and later), as well as the definitions laid out by William Kahan in his paper "Branch Cuts for Complex Elementary Functions", subtitled "Much Ado About Nothing's Sign Bit".

But what we're really seeing here is yet another divergence between the world of pure mathematics and the world of floating-point arithmetic.

The "standard" behaviour above is specific to mathematics performed using floating-point arithmetic, and in particular using floating-point arithmetic formats where there's a "negative zero" value that's distinct (equal, but distinct) from "positive zero". This includes the by now almost ubiquitous IEEE 754 floating-point standard.

The branch cuts in cmath match the "standard" mathematics ones (e.g., for asin, we make the cuts along the positive real axis from 1 to infinity, and along the negative real axis from -1 to negative infinity), and as usual, the values on the subinterval [-1, 1] of the real line match the usual real asin function in both standard mathematics and the cmath module. Given that, continuity away from the branch cuts forces cmath.asin to agree with the standard mathematical definition everywhere except possibly on the branch cuts.

Mathematically, to extend asin to the branch cuts you need to choose whether to be "continuous from above" or "continuous from below" on each cut, and the usual choice for asin is to be continuous from below on [1, inf) and continuous from above on (-inf, -1], which would give you the fourth-quadrant result for asin(10) that you expect. But if you're working with IEEE 754 floating-point, another option emerges: on the branch cut for asin, the imaginary part of the argument is always zero. You can now use the sign of the zero to determine which side of the branch cut you interpret the argument as lying on. So we get for example:

>>> from cmath import asin
>>> asin(complex(10.0, 0.0))  # 'top' of the branch cut
(1.5707963267948966+2.993222846126381j)
>>> asin(complex(10.0, -0.0))  # 'bottom' of the branch cut
(1.5707963267948966-2.993222846126381j)

This is similar to the way that the real-valued atan2 function works in most languages: usually, atan2(0.0, -1.0) is defined to be pi, while atan2(-0.0, -1.0) is defined to be -pi; the sign of the zero is used to distinguish. Mathematically, it's a bit of a cheat, but it has some nice properties in floating-point land. For example, we get that asin(z.conjugate()) is interchangeable with asin(z).conjugate() for all z, including all the floating-point special cases. And behaviour on quadrants is well-defined (if again you determine membership of each quadrant using the signs of the zeros where relevant).

As to your question about a Python module that gives the "standard" values, I'm not aware of one, though @hpaulj mentions SymPy in the comments. Or you could fudge the zero signs to force values to be interpreted the right way.