gh-113841: fix possible undefined division by 0 in _Py_c_pow()

[skirpichev]

• Issue: Possible undefined behavior division by zero in complex's _Py_c_pow() #113841

[skirpichev]

CC @gpshead

[gpshead]

cute math trick & thanks for the detailed analysis on the issue.

[miss-islington-app]

Thanks @skirpichev for the PR, and @gpshead for merging it 🌮🎉.. I'm working now to backport this PR to: 3.13.
🐍🍒⛏🤖

[bedevere-app]

GH-127216 is a backport of this pull request to the 3.13 branch.

[skirpichev]

Sorry, @gpshead, your description is slightly wrong;)

Of course, x**-y != x**y in general (for floating-point numbers), even for real. E.g.:

>>> f = 0.992856470905735
>>> math.exp(f)
2.6989328952459926
>>> 1/math.exp(-f)
2.698932895245993

Perhaps, I should mention that in the pr description. It was assumed as an ordinary truth:(

Though, the difference should be ~1ULP for good libm implementations (like in above example). So, I don't expect regressions in accuracy. My tests seems to confirm this.

But rigor error analysis is hard. In fact, I'm not sure why division was used from beginning. Closest algorithm I've found (no error checks, of course;)) - it's Algorithm 190 from ACM: https://dl.acm.org/doi/10.1145/366663.366679. It uses subtraction in the exponent.

Edit:

Here simple tests for a kind of random input.

>>> def c_pow(a, b, v):
...     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:
...         if v == 1:
...             len /= exp(at*b.imag)
...         else:
...             len *= exp(-at*b.imag)
...         phase += b.imag*log(vabs)
... 
...     return complex(len*cos(phase), len*sin(phase))
...
>>> m = 0.0
... for i in range(10000):
...     z1 = random.random() + random.random()*1j
...     z2 = random.random() + random.random()*1j
...     v1 = c_pow(z1, z2, v=1)
...     v2 = c_pow(z1, z2, v=2)
...     d = max(map(abs, [v1, v2]))
...     va = max(abs(v1), abs(v2)) 
...     m2 = d/va
...     if m2>1: 
...         print(i, z1, z2)
...         break 
...     if m2 > m:
...         m = m2
...         
>>> m
1.0

[gpshead]

LOL I can fix the commit description in the backport's squash merge. :)

[skirpichev]

Oh, yes, please. Next time I'll try not forgot to fill pr description.