GeezerVison is showing github.com — Small integer powers of real±0j are invalid · Issue #117999 · python/cpython
Small integer powers of real±0j are invalid #117999Labels interpreter-core(Objects, Python, Grammar, and Parser dirs)(Objects, Python, Grammar, and Parser dirs)type-bugAn unexpected behavior, bug, or errorAn unexpected behavior, bug, or error DescriptionFor 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
Instead, imaginary component of r should copy the sign from x.imag.
Linked PRs
Activitygaogaotiantian commented on Apr 17, 2024 Is (1-0j)**2 == 1+0j a bug? In another way, is 1+0j different than 1-0j? skirpichev commented on Apr 17, 2024
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 commented on Apr 18, 2024 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 commented on Apr 18, 2024 Wait until a decision has been made for the existing pull request. That avoids unnecessary review churn. serhiy-storchaka commented on Apr 23, 2024 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 commented on Apr 23, 2024
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 commented on Apr 23, 2024 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 commented on Apr 24, 2024 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 commented on Apr 24, 2024 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. serhiy-storchaka commented on Aug 25, 2024 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 commented on Aug 26, 2024 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 commented on Sep 19, 2024 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). skirpichev commented on Aug 18, 2026
It guarantees that c**1 is c for any complex number c. c**2 is now equivalent to c*c. etc
Unfortunately, this is not true for any c (regardless on overflows). That fail for values like shown above:
>>> x = -1-1j
>>> x*x*x*x*x*x
-8j
>>> x**6 # wrong
(-0-8j)
>>> x**8 # wrong
(16-0j)
>>> x*x*x*x*x*x*x*x
(16+0j) Note, that this also holds for the MPC, you can use gmpy2 to see it:
>>> from gmpy2 import *
>>> x = mpc('-1-1j')
>>> x*x*x*x*x*x
mpc('0.0-8.0j')
>>> x**6 # correct
mpc('0.0-8.0j')
>>> x*x*x*x*x*x*x*x
mpc('16.0+0.0j')
>>> x**8 # wrong
mpc('16.0-0.0j')
>>> (x*x*x*x)*(x*x*x*x)
mpc('16.0-0.0j')I don't see an easy way to solve this. (Above cancellation is just a simplest example, in one binop.)
So, I think it's fine for now. @serhiy-storchaka, lets merge either #124243 or #118000. serhiy-storchaka commented on Sep 3, 2026 x**6 is calculated as (x*x)*(x*x)*(x*x), and x**8 as ((x*x)*(x*x))*((x*x)*(x*x)). This is more accurate even for floats.
serhiy-storchaka commented on Sep 3, 2026 Your MPC example is version dependent: x**8 is mpc('16.0+0.0j') with MPC 1.3.1, but mpc('16.0-0.0j'), as in #124243, with MPC 1.4.0. From the MPC 1.4.0 NEWS: "mpc_pow: Agree on and implement the sign of the imaginary part when both inputs are real".
MPC does not implement complex(x, y)**0 = complex(1.0, copysign(0.0, atan2(y, x))) either: mpc(1.0,-0.0)**0 is mpc('1.0-0.0j'), but mpc(2.0,-0.0)**0 is mpc('1.0+0.0j') — the sign depends on the modulus, not on the argument. Let us keep 1+0j and consider the zero exponent in gh-137746. skirpichev commented on Sep 3, 2026
MetadataMetadatainterpreter-core(Objects, Python, Grammar, and Parser dirs)(Objects, Python, Grammar, and Parser dirs)type-bugAn unexpected behavior, bug, or errorAn unexpected behavior, bug, or error DevelopmentNo branches or pull requestsIssue actions
|