Small integer powers of real±0j are invalid

[skirpichev]

For example,

>>> z1 = complex(1,+0.0)
>>> z1**2
(1+0j)
>>> z1*z1
(1+0j)
>>> z2 = complex(1,-0.0)
>>> z2**2
(1+0j)
>>> z2*z2
(1-0j)

"Big enough" or non-integer powers are unaffected:

>>> z2**2000
(1-0j)
>>> z2**2.1
(1-0j)

It seems, this is related to

cpython/Objects/complexobject.c

Line 164 in c179c0e

r = c_1;

Instead, imaginary component of r should copy the sign from x.imag.

Linked PRs

• gh-117999: fixed small nonnegative integer powers of complex numbers #118000
• gh-117999: drop special case for small integer exponents in complex_pow() #123283
• gh-117999: Fix small integer powers of complex numbers #124243

[gaogaotiantian]

Is (1-0j)**2 == 1+0j a bug? In another way, is 1+0j different than 1-0j?

[skirpichev]

Is (1-0j)**2 == 1+0j a bug?

Huh. There is a bug, but a different one from this bugreport: 1-0j is not a complex(1,-0.0). See this thread. In short, a±bj is not an equivalent of rectangular notation for complex numbers if components are special floating point numbers (signed zeros, nans, infinities).

What is a bug here: complex(1,-0.0)**2 and complex(1,-0.0)*complex(1,-0.0) have different imaginary components:

>>> complex(1,-0.0)**2
(1+0j)
>>> complex(1,-0.0)*complex(1,-0.0)
(1-0j)

They are equal in sense of == (just as 0.0==-0.0), but you can see the difference in the repr() output (or you could use math.copysign() to test this).

[skirpichev]

The _Py_c_quot() is affected as well, not sure if this worth a separate issue:

>>> 1/complex(1, 0.)  # should be (1-0j)
(1+0j)
>>> 1/complex(1, -0.)
(1+0j)
>>> complex(1, 0.)**-1.00001
(1-0j)
>>> complex(1, 0.)**-0.99999
(1-0j)
>>> complex(1, -0.)**-1.00001
(1+0j)
>>> complex(1, -0.)**-0.99999
(1+0j)

[nineteendo]

Wait until a decision has been made for the existing pull request. That avoids unnecessary review churn.

[serhiy-storchaka]

Consider also cases when the real part is -1.

Consider also cases when the real part is ±0 and the imaginary part is ±1.

And there are also interesting cases when one of components is ±inf and other is ±1, or when both components are ±inf.

[skirpichev]

Consider also cases when the real part is -1.

Thanks, that was missed.

Consider also cases when the real part is ±0 and the imaginary part is ±1. [...] And there are also interesting cases when one of components is ±inf and other is ±1, or when both components are ±inf.

In short, all special numbers (i.e. when any component is infinity, signed zero or nan) require a special treatment. Iff... this is a bug.

@serhiy-storchaka, should I fill a separate report for _Py_c_quot()? Given

>>> complex(1, 0.0)*complex(1, 0.0)
(1+0j)
>>> complex(1, 0.0)*complex(1, -0.0)
(1+0j)

[serhiy-storchaka]

It is a different case. 1/complex(1, 0.) is currently complex(1.0, 0.0)/complex(1.0, 0.0). The sign of the imaginary part of the result depends on signs of imaginary parts of the nominator and the denominator. Since 0.0 - 0.0 and 0.0 + -0.0 are 0.0, I think that the current result is correct. When (if) we implement mixed real-complex operations separately, the right result of 1/complex(1, 0.) will be complex(1.0, -0.0). But it is a part of much larger issue, all arithmetic operations will be affected.

[skirpichev]

Then we have discontinuity, e.g.

>>> complex(1,-0.0)**1.0000001
(1-0j)
>>> complex(1,-0.0)**1  # oops
(1+0j)
>>> complex(1,-0.0)**0.9999999
(1-0j)

I think that the current result is correct.

I was thinking about definition of the inversion z from z*q = 1 == 1+0j. For z=1+0j current arithmetic rules allow both q=z and q=complex(1,-0.0). So, we should use something else (e.g continuity argument) to choose the answer.

I agree, with the mixed-mode arithmetic this decision is trivial: "just use 1/z". Perhaps, this is one argument for adding imaginary type, that was not mentioned in the d.p.o thread, unfortunately: i.e. preservation of analytic identities. For example:

>>> z = 2+1e-10j
>>> cmath.asin(z)
(1.5707963267371616+1.3169578969248166j)
>>> -1j*cmath.log(1j*z + cmath.sqrt(1-z*z))
(1.5707963267371616+1.3169578969248164j)
>>> z = 2+0j
>>> cmath.asin(z)
(1.5707963267948966+1.3169578969248166j)
>>> -1j*cmath.log(1j*z + cmath.sqrt(1-z*z))  # oops
(1.5707963267948966-1.3169578969248166j)
>>> z = 2-1e-10j
>>> cmath.asin(z)
(1.5707963267371616-1.3169578969248166j)
>>> -1j*cmath.log(1j*z + cmath.sqrt(1-z*z))
(1.5707963267371616-1.3169578969248166j)
>>> z = complex(2,-0.0)
>>> cmath.asin(z)
(1.5707963267948966-1.3169578969248166j)
>>> -1j*cmath.log(1j*z + cmath.sqrt(1-z*z))
(1.5707963267948966-1.3169578969248166j)

The C variant (on gcc+glibc) has no such issue (yet cpow() looks buggy just as CPython's pow() in this issue, hence z*z), because mixed-mode arithmetic rules at least for complex and real pairs, see this.

[skirpichev]

Ok, I did update of the pr to fix only nonnegative powers. I think, it works, here are few "differences" (which are due to different arrangements of parenthesis in the c_powu(), see below):

$ ./python a.py 
(1-1j)**4:  (-4-0j)  product:  (-4+0j)
(1-1j)**6:  (-0+8j)  product:  8j

script code

# a.py
from functools import reduce
from itertools import combinations
from math import inf, nan

cases = [complex(*x) for x in combinations([1, -1, 0.0, -0.0, 2,
                                            inf, -inf, nan], 2)]
powers = [0, 1, 2, 3, 4, 5, 6]

for c in cases:
    for p in powers:
        try:
            r_pow = repr(c**p)
        except OverflowError:
            continue
        r_pro = repr(reduce(lambda x, y: x*y, [c]*p) if p else 1+0j)
        if r_pow != r_pro:
            print(f'{c}**{p}: ', r_pow, " product: ", r_pro)

For illustration of a failing example:

>>> z = 1-1j
>>> z*z
-2j
>>> z*z*z*z
(-4+0j)
>>> (1+0j)*(z*z)*(z*z)  # z**4 implementation
(-4-0j)

As associativity is broken already for floats - I think it's fine.

[skirpichev]

Hmm, it seems that similar issue for negative (small) integer exponents can't be solved without mixed-mode arithmetic in some parts.

Example in the main:

>>> z = 1+0j
>>> z**-111.0000001
(1-0j)
>>> z**-111  # big enough to use generic power algorithm
(1-0j)
>>> z**-110.9999999
(1-0j)
>>> z**-1.0000001
(1-0j)
>>> z**-1  # oops
(1+0j)
>>> c_pow(z, -1)  # with generic algorithm in python
(1-0j)
>>> z**-0.9999999
(1-0j)

Same happens for other special components. This seems to be an issue, regardless on the question of how we define e.g. 1/(1+0j).

The problem part in this case is this line:

cpython/Objects/complexobject.c

Line 181 in 3c890b5

return _Py_c_quot(c_1, c_powu(x,-n));

Assuming that c_powu() is fixed, true division on this line should be implemented using mixed-mode arithmetic rules (real / complex), i.e. 1/z=complex(z.real,-z.imag)/abs(z)**2.

pure-python version of generic pow() algorithm

def c_pow(a, b):
    from math import atan2, cos, exp, hypot, log, sin
    if b.real == b.imag == 0:
        return 1+0j
    elif a.real == a.imag == 0:
        if b.imag or b.real < 0:
            raise ZeroDivisionError
        return 0j

    vabs = hypot(a.real, a.imag)
    len = pow(vabs, b.real)
    at = atan2(a.imag, a.real)
    phase = at*b.real
    if b.imag:
        len /= exp(at*b.imag)
        phase += b.imag*log(vabs)

    return complex(len*cos(phase), len*sin(phase))

[skirpichev]

See also #60200 (comment) example from SO (last case lacks OverflowError in the current main).

[serhiy-storchaka]

The problem in the current algorithm for small integer degree is that it assumes that c**0 is c_1 = complex(1.0, 0.0) for every complex number c. But it should be complex(1.0, -0.0) for complex(1.0, -0.0). There is a similar error in calculation for the negative integer degree.

For $c = r e^{i\phi}$, $\lim_{\varepsilon \to 0}{c^\varepsilon} = \lim_{\varepsilon \to 0}{e^{\varepsilon(\ln{r}+i\phi)}} = \lim_{\varepsilon \to 0}{e^{\varepsilon\ln{r}}} \cdot \lim_{\varepsilon \to 0}{e^{i\varepsilon\phi}} = \lim_{\varepsilon \to 0}{e^{i\varepsilon\phi}} = 1 + i\lim_{\varepsilon \to 0}{\varepsilon\phi}$.

So, complex(x, y)**0 should be complex(1.0, copysign(0.0, atan2(y, x)).

[skirpichev]

Unfortunately, I doubt that binary exponentiation can be easily (if at all) fixed to address this. Another example from:

>>> import _testcapi
>>> z = -1-1j
>>> _testcapi._py_c_pow(z, 6)[0]
(4.408728476930473e-15-8.000000000000004j)
>>> z**6
(-0-8j)
>>> z*z*z*z*z*z
-8j

[serhiy-storchaka]

See also #124243 which fixes calculation of small integer powers of complex numbers. It guarantees that c**1 is c for any complex number c. c**2 is now equivalent to c*c. etc (if no OverflowError was raised).