Skip to content

gh-117999: fixed small integer powers of complex numbers - #118000

Open
skirpichev wants to merge 16 commits into
python:mainfrom
skirpichev:fix-c_powi-117999
Open

gh-117999: fixed small integer powers of complex numbers#118000
skirpichev wants to merge 16 commits into
python:mainfrom
skirpichev:fix-c_powi-117999

Conversation

@skirpichev

@skirpichev skirpichev commented Apr 17, 2024 •
edited
Loading

Copy link
Copy Markdown
Member

Before, handling of numbers with special values in components (infinities, nans, signed zero) was invalid. Simple example:

    >>> z = complex(1, -0.0)
    >>> z*z
    (1-0j)
    >>> z**2
    (1+0j)
    >>> complex('inf')**-1
    0j

Now:

    >>> z**2
    (1-0j)
    >>> complex('inf')**-1
    -0j

Also, positive powers of infinite complex not raise anymore OverflowError's:

   >>> complex('inf')**1
   (inf+0j)
   >>> complex('inf')**0.5
   (inf+nanj)

@bedevere-app bedevere-app Bot mentioned this pull request Apr 17, 2024
Open
@skirpichev
skirpichev force-pushed the fix-c_powi-117999 branch 2 times, most recently from 2a6babf to 5ee43ee Compare April 18, 2024 08:49
@skirpichev skirpichev changed the title gh-117999: fixed invalid small integer powers of real±0j gh-117999: fixed small nonnegative integer powers of complex numbers with special components Apr 24, 2024
@skirpichev
skirpichev force-pushed the fix-c_powi-117999 branch 2 times, most recently from 51c6cad to cd3e11e Compare May 12, 2024 17:00
@skirpichev skirpichev changed the title gh-117999: fixed small nonnegative integer powers of complex numbers with special components gh-117999: fixed small nonnegative integer powers of complex numbers May 12, 2024
@skirpichev skirpichev mentioned this pull request May 29, 2024
Merged
skirpichev added 2 commits May 30, 2024 07:31
@skirpichev
pythongh-117999: fixed small nonnegative integer powers of complex nu…
70d21c7 
…mbers

Before, handling of numbers with special values in components
(infinities, nans, signed zero) was invalid.  Simple example:

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

Now:

    >>> z**2
    (1-0j)
@skirpichev
Polishing of c_powu(): exclude unreachible cases and magic numbers
8e0c482
@skirpichev
skirpichev force-pushed the fix-c_powi-117999 branch from 58c2fb1 to 8e0c482 Compare May 30, 2024 04:31
@skirpichev skirpichev mentioned this pull request Jul 8, 2024
Open
@skirpichev

skirpichev commented Aug 18, 2024

Copy link
Copy Markdown
Member Author

@picnixz, I would appreciate your review on this pr. Or your opinion in the issue thread.

@picnixz

picnixz commented Aug 18, 2024

Copy link
Copy Markdown
Member

I'll do it tomorrow! (Monday, Paris time)

@picnixz
picnixz self-requested a review August 18, 2024 13:54
picnixz
picnixz reviewed Aug 19, 2024
View reviewed changes
Comment thread Lib/test/test_complex.py Outdated
Comment thread Lib/test/test_complex.py
Comment thread Lib/test/test_complex.py Outdated
Comment thread Lib/test/test_complex.py Outdated
Comment thread Lib/test/test_complex.py Outdated
Comment thread Objects/complexobject.c Outdated
Comment thread Lib/test/test_complex.py Outdated
@picnixz

picnixz commented Aug 19, 2024

Copy link
Copy Markdown
Member

I need to think a bit more on the issue. I'll try to have something by the end of the day or tomorrow. Ideally, I would like to have no inconsistency between the generic algorithm and the non-generic one (namely, the result should be as if we were using the generic algorithm).

@skirpichev

skirpichev commented Aug 19, 2024 •
edited
Loading

Copy link
Copy Markdown
Member Author

I would like to have no inconsistency between the generic algorithm and the non-generic one

I'm not sure if it's possible without too much code, that affects performance severely. BTW, I think that numpy code has no such special version for integer exponents. I'll double check.

Edit: Ah, no. numpy mimics CPython here, at least in npy_math_complex.c.src.

Edit2:

JFR, some simple benchmarks.

With specialized code (main):

$ python -m pyperf timeit -q -s 'z=1+1j;e=20' 'pow(z,e)'
Mean +- std dev: 237 ns +- 5 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=90' 'pow(z,e)'
Mean +- std dev: 259 ns +- 4 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=110' 'pow(z,e)'
Mean +- std dev: 579 ns +- 5 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=200' 'pow(z,e)'
Mean +- std dev: 571 ns +- 3 ns

Without:

$ python -m pyperf timeit -q -s 'z=1+1j;e=2' 'pow(z,e)'
Mean +- std dev: 562 ns +- 3 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=20' 'pow(z,e)'
Mean +- std dev: 576 ns +- 2 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=90' 'pow(z,e)'
Mean +- std dev: 576 ns +- 3 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=110' 'pow(z,e)'
Mean +- std dev: 578 ns +- 4 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=200' 'pow(z,e)'
Mean +- std dev: 573 ns +- 3 ns

skirpichev
skirpichev commented Aug 19, 2024
View reviewed changes
Comment thread Lib/test/test_complex.py Outdated
skirpichev
skirpichev commented Aug 19, 2024
View reviewed changes
Comment thread Objects/complexobject.c Outdated
skirpichev
skirpichev commented Aug 19, 2024
View reviewed changes
Comment thread Lib/test/test_complex.py Outdated
skirpichev
skirpichev commented Aug 19, 2024
View reviewed changes
Comment thread Lib/test/test_complex.py Outdated
skirpichev and others added 3 commits August 19, 2024 11:55
@skirpichev @picnixz
Apply suggestions from code review
389441d 
Co-authored-by: Bénédikt Tran <10796600+picnixz@users.noreply.github.com>
@skirpichev @picnixz
Apply suggestions from code review
9f24ad1 
Co-authored-by: Bénédikt Tran <10796600+picnixz@users.noreply.github.com>
@skirpichev
Address review: subTest and asserts
807bcc6
@skirpichev skirpichev closed this Aug 20, 2024
@skirpichev
skirpichev deleted the fix-c_powi-117999 branch August 20, 2024 03:29
@skirpichev
skirpichev restored the fix-c_powi-117999 branch September 17, 2025 04:29
@skirpichev skirpichev reopened this Sep 17, 2025
@skirpichev
skirpichev marked this pull request as draft September 17, 2025 04:29
@github-actions

github-actions Bot commented Apr 16, 2026

Copy link
Copy Markdown

This PR is stale because it has been open for 30 days with no activity.

@github-actions github-actions Bot added the stale Stale PR or inactive for long period of time. label Apr 16, 2026
@skirpichev skirpichev self-assigned this May 24, 2026
@github-actions github-actions Bot removed the stale Stale PR or inactive for long period of time. label May 30, 2026
@skirpichev
skirpichev marked this pull request as ready for review August 18, 2026 08:28
@skirpichev

skirpichev commented Aug 18, 2026

Copy link
Copy Markdown
Member Author

@picnixz, this is ready for review again.

See #117999 (comment) - we still can't avoid differences wrt generic algorithm for all values. But probably it's fine.

@skirpichev

skirpichev commented Aug 19, 2026 •
edited
Loading

Copy link
Copy Markdown
Member Author

I run on this following script, comparing pow() output wrt GNU MPC. This PR (and #124243 as well) introduce additional differences (32 broken test):

--- 1   2026-08-20 03:33:59.886826488 +0300
+++ 2   2026-08-20 03:58:01.747376340 +0300
@@ -6,6 +6,8 @@
 (-inf-infj)  (-0.5+nanj)  0j           (nan+nanj)  
 (-inf-infj)  (-0-0j)      (1-0j)       (1+0j)      
 (-inf-infj)  -0j          (1-0j)       (1+0j)      
+(-inf-infj)  (0.5-0j)     (inf+nanj)   (inf-infj)  
+(-inf-infj)  (0.5+0j)     (inf+nanj)   (inf-infj)  
 (-inf-infj)  (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (-inf-infj)  (inf+0.5j)   (inf+nanj)   (nan+nanj)  
 (-inf-infj)  (inf+infj)   (inf+nanj)   (nan+nanj)  
@@ -16,6 +18,8 @@
 (-inf-0.5j)  (-0.5+nanj)  0j           (nan+nanj)  
 (-inf-0.5j)  (-0-0j)      (1-0j)       (1+0j)      
 (-inf-0.5j)  -0j          (1-0j)       (1+0j)      
+(-inf-0.5j)  (0.5-0j)     (inf+nanj)   (inf-infj)  
+(-inf-0.5j)  (0.5+0j)     (inf+nanj)   (inf-infj)  
 (-inf-0.5j)  (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (-inf-0.5j)  (inf+0.5j)   (inf+nanj)   (nan+nanj)  
 (-inf-0.5j)  (inf+infj)   (inf+nanj)   (nan+nanj)  
@@ -26,6 +30,8 @@
 (-inf-0j)    (-0.5+nanj)  0j           (nan+nanj)  
 (-inf-0j)    (-0-0j)      (1-0j)       (1+0j)      
 (-inf-0j)    -0j          (1-0j)       (1+0j)      
+(-inf-0j)    (0.5-0j)     (inf+nanj)   (inf-infj)  
+(-inf-0j)    (0.5+0j)     (inf+nanj)   (inf-infj)  
 (-inf-0j)    (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (-inf-0j)    (inf+0.5j)   (inf+nanj)   (nan+nanj)  
 (-inf-0j)    (inf+infj)   (inf+nanj)   (nan+nanj)  
@@ -38,6 +44,8 @@
 (-inf+0j)    (-0.5+nanj)  0j           (nan+nanj)  
 (-inf+0j)    (-0-0j)      (1-0j)       (1+0j)      
 (-inf+0j)    -0j          (1-0j)       (1+0j)      
+(-inf+0j)    (0.5-0j)     (inf+nanj)   (inf+infj)  
+(-inf+0j)    (0.5+0j)     (inf+nanj)   (inf+infj)  
 (-inf+0j)    (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (-inf+0j)    (inf-infj)   (inf+nanj)   (nan+nanj)  
 (-inf+0j)    (inf-0.5j)   (inf+nanj)   (nan+nanj)  
@@ -50,6 +58,8 @@
 (-inf+0.5j)  (-0.5+nanj)  0j           (nan+nanj)  
 (-inf+0.5j)  (-0-0j)      (1-0j)       (1+0j)      
 (-inf+0.5j)  -0j          (1-0j)       (1+0j)      
+(-inf+0.5j)  (0.5-0j)     (inf+nanj)   (inf+infj)  
+(-inf+0.5j)  (0.5+0j)     (inf+nanj)   (inf+infj)  
 (-inf+0.5j)  (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (-inf+0.5j)  (inf-infj)   (inf+nanj)   (nan+nanj)  
 (-inf+0.5j)  (inf-0.5j)   (inf+nanj)   (nan+nanj)  
@@ -62,6 +72,8 @@
 (-inf+infj)  (-0.5+nanj)  0j           (nan+nanj)  
 (-inf+infj)  (-0-0j)      (1-0j)       (1+0j)      
 (-inf+infj)  -0j          (1-0j)       (1+0j)      
+(-inf+infj)  (0.5-0j)     (inf+nanj)   (inf+infj)  
+(-inf+infj)  (0.5+0j)     (inf+nanj)   (inf+infj)  
 (-inf+infj)  (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (-inf+infj)  (inf-infj)   (inf+nanj)   (nan+nanj)  
 (-inf+infj)  (inf-0.5j)   (inf+nanj)   (nan+nanj)  
@@ -98,6 +110,8 @@
 (-0.5-infj)  (-0.5+nanj)  0j           (nan+nanj)  
 (-0.5-infj)  (-0-0j)      (1-0j)       (1+0j)      
 (-0.5-infj)  -0j          (1-0j)       (1+0j)      
+(-0.5-infj)  (0.5-0j)     (inf+nanj)   (inf-infj)  
+(-0.5-infj)  (0.5+0j)     (inf+nanj)   (inf-infj)  
 (-0.5-infj)  (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (-0.5-infj)  (inf+0.5j)   (inf+nanj)   (nan+nanj)  
 (-0.5-infj)  (inf+infj)   (inf+nanj)   (nan+nanj)  
@@ -142,6 +156,8 @@
 (-0.5+infj)  (-0.5+nanj)  0j           (nan+nanj)  
 (-0.5+infj)  (-0-0j)      (1-0j)       (1+0j)      
 (-0.5+infj)  -0j          (1-0j)       (1+0j)      
+(-0.5+infj)  (0.5-0j)     (inf+nanj)   (inf+infj)  
+(-0.5+infj)  (0.5+0j)     (inf+nanj)   (inf+infj)  
 (-0.5+infj)  (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (-0.5+infj)  (inf-infj)   (inf+nanj)   (nan+nanj)  
 (-0.5+infj)  (inf-0.5j)   (inf+nanj)   (nan+nanj)  
@@ -154,6 +170,8 @@
 (-0-infj)    (-0.5+nanj)  0j           (nan+nanj)  
 (-0-infj)    (-0-0j)      (1-0j)       (1+0j)      
 (-0-infj)    -0j          (1-0j)       (1+0j)      
+(-0-infj)    (0.5-0j)     (inf+nanj)   (inf-infj)  
+(-0-infj)    (0.5+0j)     (inf+nanj)   (inf-infj)  
 (-0-infj)    (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (-0-infj)    (inf+0.5j)   (inf+nanj)   (nan+nanj)  
 (-0-infj)    (inf+infj)   (inf+nanj)   (nan+nanj)  
@@ -186,6 +204,8 @@
 (-0+infj)    (-0.5+nanj)  0j           (nan+nanj)  
 (-0+infj)    (-0-0j)      (1-0j)       (1+0j)      
 (-0+infj)    -0j          (1-0j)       (1+0j)      
+(-0+infj)    (0.5-0j)     (inf+nanj)   (inf+infj)  
+(-0+infj)    (0.5+0j)     (inf+nanj)   (inf+infj)  
 (-0+infj)    (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (-0+infj)    (inf-infj)   (inf+nanj)   (nan+nanj)  
 (-0+infj)    (inf-0.5j)   (inf+nanj)   (nan+nanj)  
@@ -198,6 +218,8 @@
 -infj        (-0.5+nanj)  0j           (nan+nanj)  
 -infj        (-0-0j)      (1-0j)       (1+0j)      
 -infj        -0j          (1-0j)       (1+0j)      
+-infj        (0.5-0j)     (inf+nanj)   (inf-infj)  
+-infj        (0.5+0j)     (inf+nanj)   (inf-infj)  
 -infj        (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 -infj        (inf+0.5j)   (inf+nanj)   (nan+nanj)  
 -infj        (inf+infj)   (inf+nanj)   (nan+nanj)  
@@ -230,6 +252,8 @@
 infj         (-0.5+nanj)  0j           (nan+nanj)  
 infj         (-0-0j)      (1-0j)       (1+0j)      
 infj         -0j          (1-0j)       (1+0j)      
+infj         (0.5-0j)     (inf+nanj)   (inf+infj)  
+infj         (0.5+0j)     (inf+nanj)   (inf+infj)  
 infj         (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 infj         (inf-infj)   (inf+nanj)   (nan+nanj)  
 infj         (inf-0.5j)   (inf+nanj)   (nan+nanj)  
@@ -242,6 +266,8 @@
 (0.5-infj)   (-0.5+nanj)  0j           (nan+nanj)  
 (0.5-infj)   (-0-0j)      (1-0j)       (1+0j)      
 (0.5-infj)   -0j          (1-0j)       (1+0j)      
+(0.5-infj)   (0.5-0j)     (inf+nanj)   (inf-infj)  
+(0.5-infj)   (0.5+0j)     (inf+nanj)   (inf-infj)  
 (0.5-infj)   (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (0.5-infj)   (inf+0.5j)   (inf+nanj)   (nan+nanj)  
 (0.5-infj)   (inf+infj)   (inf+nanj)   (nan+nanj)  
@@ -304,6 +330,8 @@
 (0.5+infj)   (-0.5+nanj)  0j           (nan+nanj)  
 (0.5+infj)   (-0-0j)      (1-0j)       (1+0j)      
 (0.5+infj)   -0j          (1-0j)       (1+0j)      
+(0.5+infj)   (0.5-0j)     (inf+nanj)   (inf+infj)  
+(0.5+infj)   (0.5+0j)     (inf+nanj)   (inf+infj)  
 (0.5+infj)   (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (0.5+infj)   (inf-infj)   (inf+nanj)   (nan+nanj)  
 (0.5+infj)   (inf-0.5j)   (inf+nanj)   (nan+nanj)  
@@ -316,6 +344,8 @@
 (inf-infj)   (-0.5+nanj)  0j           (nan+nanj)  
 (inf-infj)   (-0-0j)      (1-0j)       (1+0j)      
 (inf-infj)   -0j          (1-0j)       (1+0j)      
+(inf-infj)   (0.5-0j)     (inf+nanj)   (inf-infj)  
+(inf-infj)   (0.5+0j)     (inf+nanj)   (inf-infj)  
 (inf-infj)   (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (inf-infj)   (inf+0.5j)   (inf+nanj)   (nan+nanj)  
 (inf-infj)   (inf+infj)   (inf+nanj)   (nan+nanj)  
@@ -404,6 +434,8 @@
 (inf+infj)   (-0.5+nanj)  0j           (nan+nanj)  
 (inf+infj)   (-0-0j)      (1-0j)       (1+0j)      
 (inf+infj)   -0j          (1-0j)       (1+0j)      
+(inf+infj)   (0.5-0j)     (inf+nanj)   (inf+infj)  
+(inf+infj)   (0.5+0j)     (inf+nanj)   (inf+infj)  
 (inf+infj)   (0.5+nanj)   (inf+nanj)   (nan+nanj)  
 (inf+infj)   (inf-infj)   (inf+nanj)   (nan+nanj)  
 (inf+infj)   (inf-0.5j)   (inf+nanj)   (nan+nanj)  
@@ -497,4 +529,4 @@
 (nan+nanj)   (-0-0j)      (1-0j)       (1+0j)      
 (nan+nanj)   (-0+0j)      (1-0j)       (1+0j)      
 ====================================================
-Failures: 496 (21%)
+Failures: 528 (22%)

I think it's fine, assuming that on above examples output must match cmath.sqrt(base).

@skirpichev skirpichev mentioned this pull request Aug 19, 2026
Open
@skirpichev skirpichev changed the title gh-117999: fixed small nonnegative integer powers of complex numbers gh-117999: fixed small integer powers of complex numbers Aug 20, 2026
@skirpichev skirpichev removed their assignment Aug 24, 2026
hpkfft
hpkfft reviewed Aug 30, 2026
View reviewed changes

@hpkfft hpkfft left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This PR looks good to me, and I would approve it if I had the necessary access to do so.

Without this PR, (inf+0j)**(1+0j) raises an OverflowError, which I think is broken. It's bad to raise an error, since it can crash an entire Python program if the programmer is not expecting it to be possible. Note that for the real-valued case:

>>> math.inf**1
inf

Comment thread Lib/test/test_complex.py
continue
for n in range(1, 9):
with self.subTest(exponent=-n):
self.assertComplexesAreIdentical(c**-n, 1/(c**n))

@hpkfft hpkfft Aug 30, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This assertion is an implementation detail. In the future, for example, one might want to use the C math library's cpow for both positive and negative small integers, in which case the two expressions may not be identical.
If this test is seen as documenting an eternal, unbreakable contract with users, then I don't like it. On the other hand, the understanding might be that Python can change the implementation and then simply change the test as well. I'm perfectly fine with that. (I'm not sure of the development culture here.)

I guess the comment about identical also applies to c**2 == c*c and maybe even c**1 == c. In MPC, (inf+1j)**(1+0j) gives (inf+nanj)....

@skirpichev

skirpichev commented Aug 31, 2026

Copy link
Copy Markdown
Member Author

Without this PR, (inf+0j)**(1+0j) raises an OverflowError, which I think is broken.

Perhaps, we should factor out this fix, as it looks less controversial. Problem noticed e.g. here: #60200 (comment)

@skirpichev skirpichev mentioned this pull request Aug 31, 2026
Open
@skirpichev skirpichev self-assigned this Aug 31, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

@picnixz picnixz Awaiting requested review from picnixz
@serhiy-storchaka serhiy-storchaka Awaiting requested review from serhiy-storchaka
1 more reviewer
@hpkfft hpkfft hpkfft left review comments
Reviewers whose approvals may not affect merge requirements

Assignees

@skirpichev skirpichev

Labels

awaiting review

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

3 participants

@skirpichev @picnixz @hpkfft