BUG: arcsin evaluation on branch cuts deviates from Python array API standard

[pearu]

Describe the issue:

According to Python array API standard, arcsin(z) is defined as

-1j * ln(z*1j + sqrt(1 - z**2))

As demonstrated below, the evaluation of numpy.arcsin and using the definition above lead to different results on the real line when abs(z.real) > 1.

Reproduce the code example:

import numpy
z = complex(2, 0)
print(numpy.arcsin(z))                            # -> (1.5707963267948966+1.3169578969248166j)
print(-1j*numpy.log(z*1j + numpy.sqrt(1 - z*z)))  # -> (1.5707963267948966-1.3169578969248166j)

Error message:

No response

Python and NumPy Versions:

2.0.0rc1
3.11.8 | packaged by conda-forge | (main, Feb 16 2024, 20:53:32) [GCC 12.3.0]

Runtime Environment:

No response

Context for the issue:

See also mpmath/mpmath#786 that contains some data points on how other software evaluate arcsin on its branch cuts.

[skirpichev]

I doubt it's a numpy issue. Rather it's a known issue of the CPython's complex arithmetic, see this thread. The C11 standard (annex G) has special rules for mixed-mode arithmetic (real+complex and imaginary+complex; most compilers implement at least the first), while the CPython's complex numbers do type casts (real -> complex(real, 0), etc). Your example could be corrected:

>>> import numpy
>>> z = complex(2, 0)
>>> numpy.arcsin(z)
(1.5707963267948966+1.3169578969248166j)
>>> -1j*numpy.log(z*1j + numpy.sqrt(complex(1 - (z*z).real, -(z*z).imag)))
(1.5707963267948966+1.3169578969248164j)

[pearu]

It is a numpy issue in the sense that the actual behavior of numpy.arcsin does not match with the (target) promise that NumPy is fully Python array API standard compliant. Although, such issue could be fixed at the documentation level.

However, it is not solely of NumPy or CPython issue. The same behavior is present in Boost asin implementation (because the same Hull et al algorithm is used there), for instance.

The corrected example is not sufficient to resolve this issue as users are not expected to manipulate their code in the demonstrated way.

[skirpichev]

the actual behavior of numpy.arcsin does not match with the (target) promise that NumPy is fully Python array API standard compliant.

I doubt. The documentation mention mathematical expression, \[\operatorname{asin}(z) = -j\ \ln(zj + \sqrt{1-z^2})\] (not it's Python version). The actual behaviour of arcsin() correctly match this. The -1j * ln(z*1j + sqrt(1 - z**2)) is an incorrect translation of the above expression to Python, it fails to handle special numbers (like complex(2,+0.0)) properly.

Same issue affects other elementary inverse functions, e.g. arcsinh():

>>> numpy.arcsinh(complex(0.0,2))
(1.3169578969248166+1.5707963267948966j)
>>> numpy.arcsinh(complex(-0.0,2))
(-1.3169578969248166+1.5707963267948966j)
>>> z1 = complex(0.0,2)
>>> z=z1; numpy.log(z + numpy.sqrt(1+z*z))
(1.3169578969248166+1.5707963267948966j)
>>> z2 = complex(-0.0,2)
>>> z=z2; numpy.log(z + numpy.sqrt(1+z*z))
(1.3169578969248166+1.5707963267948966j)

Although, such issue could be fixed at the documentation level.

Not sure how it's possible.

Proper fix would be adding an imaginary type in CPython and implementing all mixed-mode arithmetic rules (like it's specified in the c11's annex G). This will require a PEP. From the quoted thread, I doubt someone from CPython core devs will sponsor this...

However, it is not solely of NumPy or CPython issue. The same behavior is present in Boost asin implementation

Probably, it depends on the complex number backend (mpc, for example). Not all implement special rules for mixed-mode arithmetic. This is valid for the gmpy2:

>>> import gmpy2
>>> z = gmpy2.mpc(2, 0)
>>> gmpy2.mpc(0,-1) * gmpy2.log(z*gmpy2.mpc(0,1) + gmpy2.sqrt(1 - z*z))
mpc('1.5707963267948966-1.3169578969248166j')

[pearu]

... it fails to handle special numbers (like complex(2,+0.0)) properly.

OK, I see. So, this issue is a consequence of having defects in the basic complex arithmetic of Python or numpy objects. Say, take

>>> z = numpy.array([2, 2], dtype=numpy.complex128)
>>> z.imag[0] *= -1
>>> z
array([2.-0.j, 2.+0.j])
>>> I = numpy.array([0], dtype=numpy.complex128)
>>> I.imag[0] = 1
>>> I
array([0.+1.j])

then the sign bit gets easily lost even in non-mixed-mode arithmetic:

>>> z * I
array([0.+2.j, 0.+2.j])

Although, such issue could be fixed at the documentation level.
Not sure how it's possible.

One possibility is to adjust the Python array API standard to use a different definition for arcsin, say, the one that Hull et al paper uses. This definition ensures that sign(arcsin(z).imag) == sign(z.imag) on the arcsin branch cut.

Another possibility is to mention in the Python array API standard that on the branch cut, the arcsin expression is not guaranteed to evaluate to the same value with the arcsin return value when expression evaluation is signed-zero ignorant.

Same issue affects other elementary inverse functions, e.g. arcsinh

Absolutely. The reason is that in the set of functions {arcsinh, arccosh, arcsin, arccos}, one often implements only one of them, say, arcsin, and the rest of the functions can be expressed in terms of arcsin.

[skirpichev]

then the sign bit gets easily lost even in non-mixed-mode arithmetic z * I

In fact, this is also case of mixed-mode arithmetic:) The C11 has special imaginary types (the macro I expands to this) and special rules for imaginary x complex operations, namely:

(x + iy) * (iv) = (-y*v) + i(x*v)

Though, neither gcc nor clang support this.

[seberg]

There doesn't seem to be anything actionable to me here in NumPy, so going to close it. The "array api" even says:

Note: branch cuts follow C99 and have provisional status

basically with the assumption that C99 is a reasonable expectation. If NumPy fails that, it would most likely (admittedly not necessarily) mean that compilers/math libs used by NumPy ignore C99.
And I would probably need a far better reason than a note in the Array API that was written with the assumption that C99 is widely adopted. For one, musl complex functions were pretty imprecise IIRC, and for the most part we just live with it (at least the longdouble versions).

---

I am certain there are some interesting things here, like whether complex_arr * real should actually be implemented as complex_arr * complex(real, -0.0) or need a dedicated loop to get -0 right?
But even if we want to do something about it, it seems like a different issue.

[skirpichev]

whether complex_arr * real should actually be implemented as complex_arr * complex(real, -0.0)

(1) complex_arr*complex(real, 0.0) actually (as in CPython, Go, Julia, etc).

These coercion rules break using usual rectangular notation for complex numbers (e.g. in CPython: 1-0.0j is 1+0.0j). But also they break common mathematical identities like above.

or need a dedicated loop to get -0 right

(2) Or follow C99, C11, etc. I.e. (a+ib)*r = a*r+i(b*r).

it seems like a different issue.

On another hand, as it seems from the linked thread, CPython core devs are mostly happy with the present state of art (1) and I think, that any change (i.e. adoption of c99+, (2)) here could come only from numeric libraries, outside of the stdlib.

[seberg]

here could come only from numeric libraries, outside of the stdlib

Right, the scientific python community could certainly push for such changes and may succeed with sufficient motivation.
However, the Array API assuming that following C99 makes sense (when things seem trickier on multiple levels) isn't a motivation/reason for it.
But the Python discussion and you have a point: The only true fix, apparently requires the existence of imaginary-only scalars which adds a lot of complexity... but for whose gain exactly?

The main thing is: While I agree in principle with all of these subtleties, the only actual issues I remember seeing are from those explicitly writing tests for them with no motivation beyond testing standards like C99 compliance.
Although maybe you know about some actual issues biting users...

(This is not isolated to imaginary numbers, similar things exist with -0.0 handling for float in some places.)

[rgommers]

I agree with @seberg here. Branch cuts are nice to have clearly documented when every array/tensor library is already compliant. The requirements were checked during writing of the spec standard, but there are so many that surely a few mismatches were missed. In such cases, it seems perfectly fine to skip the relevant test and document what the mismatch is - it can be revisited if there's a more compelling reason to put actual work into it.

[pearu]

The main thing is: While I agree in principle with all of these subtleties, the only actual issues I remember seeing are from those explicitly writing tests for them with no motivation beyond testing standards like C99 compliance. Although maybe you know about some actual issues biting users...

If I understand your message correctly, this is a disappointing conclusion for building software quality upon testing (where standards can be used as a reference). In this particular case, the advice seems to be to use

expected = complex(1.5707963267948966, 1.3169578969248166)  # or use arcsin from a reference library such as numpy or other
assert arcsin(complex(2, 0)) in {expected, expected.conjugate()}

approach which is more annoying compared to agreeing what is the expected result and implementing it accordingly, especially when relevant standards already exists. I realize that this is not an easy problem to resolve, say, within Python array API standard, considering that there exists two groups of major software that define the arcsin value at its branch cut differently for various reasons:

• in symbolic software like Mathematica, Maple, Maxima, CLN, SymPy, Matlab, mpmath.mp asin(2) returns 1.5608 - 1.3170j
• in floating-point software like NumPy, cmath, PyTorch, Boost, mpmath.fp, and gmpy2, asin(2) returns 1.5608 + 1.3170j

When implementing arcsin for a third software, a decision must be made to which group to belong which is the actual issue that I am facing atm.

[seberg]

Well, I don't mind if other NumPy devs take up this. But basically, similar libraries seem to agree with NumPy and many of them have far more experience and even if we agree it would make more sense to follow mpmath (which anyone is free to chime in there).

More importantly: I also think there is also way too much half knowledge here and that explicitly includes the NumPy implementation and the Array API.

So if you want to make the point, then you really need to dig deeper, and the only thing I can think of is the C99 Annex G (short of some old Kahan paper or so).
The only thing that I am able to find there that seems applicable is (I may have missed something), is that the complex arccos is specified fully w.r.t. to branch cuts and it includes, this as the first bullet point:

• cacos(conj(z)) = conj(cacos(z))

Which, translated, means that:

np.arccos(complex(2, -0.)) == np.arccos(complex(2, 0.)).conjugate()

and thus the branch cut must differentiate between -0.0 and NumPy is correct.
And honestly, that has also the very valuable property that if you underflow the floating point value while convergin towards 0 you will not get a discontinuity!

EDIT: That Kahan reference, which to me would be authoritative is: "Branch Cuts for Complex Elementary Functions or Much Ado About Nothing's Sign Bit", which I think agrees, but I don't want to read it that in depth...

[skirpichev]

The only true fix, apparently requires the existence of imaginary-only scalars which adds a lot of complexity...

That's only part of the fix. Another one is special mixed-mode rules for arithmetic, that include "complex op real" cases and "complex op imaginary" cases. Exactly as specified by Annex G of the C99.

but for whose gain exactly?

The main thing is: While I agree in principle with all of these
subtleties, the only actual issues I remember seeing are from those
explicitly writing tests for them with no motivation beyond testing
standards like C99 compliance.

Well, motivation is mostly same as for the C standard.

First, we could keep ordinary rectangular notation for complex numbers, without corner cases for special numbers (signed zero, infinities, nans). E.g. to input complex(1,-0.0) you should use something like -(1+0j).

Second, we could preserve common identities like asin(z)=-1j*log(1j*z+sqrt(1-z**2). Without this, every identity, including real and imaginary literals could be easily broken for special numbers, like z=2+0j. "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

Sadly, the second argument was missing in the mentioned thread.

Although maybe you know about some actual issues biting users...

Well, many people assume that usual rectangular notation is valid for CPython's complex numbers. (BTW, this is claimed in the CPython docs! See python/cpython#109218) The quoted thread links many CPython's issues, here is another recent example in the mpmath: mpmath/mpmath#774 (comment) (this time related to infinities, not signed zeros).

(This is not isolated to imaginary numbers, similar things exist with -0.0
handling for float in some places.)

Floats in CPython mostly are just IEEE 754 floats, I doubt there are much place for confusion.

[seberg]

Fair, unlike the branch-cuts (which everything I saw until now make me think they are right as they are), there are issues around mixed mode operations.

I am not sure I want to follow up on it in NumPy (I think this is something where it would be nice to have Python/others take the lead before taking any real action).
But maybe it is an issue, although asking on the mailing list or SciPy discussion is probably a better venue to spark interest in the Python discussion.