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: use generic algorithm in complex_pow() if base has special components #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.