Skip to content

Commit 97358c7

Browse files
skirpichevskirpichev
committed
pythongh-120010: fix invalid (nan+nanj) results in _Py_c_prod()
In some cases, previosly computed as (nan+nanj), we could recover meaningful component values in the result, see e.g. the C11, Annex G.5.2, routine _Cmuld(): >>> z = 1e300+1j >>> z*(inf+infj) # was (nan+nanj) (nan+infj) That also fix some complex powers for small integer exponents, computed by optimised algorithm (by squaring): >>> z**5 # was (nan+nanj) Traceback (most recent call last): File "<python-input-1>", line 1, in <module> z**5 ~^^~ OverflowError: complex exponentiation
1 parent aa9fe98 commit 97358c7

File tree

3 files changed

+85
-3
lines changed

3 files changed

+85
-3
lines changed

‎Lib/test/test_complex.py

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -94,6 +94,10 @@ def assertFloatsAreIdentical(self, x, y):
9494
msg += ': zeros have different signs'
9595
self.fail(msg.format(x, y))
9696

97+
def assertComplexesAreIdentical(self, x, y):
98+
self.assertFloatsAreIdentical(x.real, y.real)
99+
self.assertFloatsAreIdentical(x.imag, y.imag)
100+
97101
def assertClose(self, x, y, eps=1e-9):
98102
"""Return true iff complexes x and y "are close"."""
99103
self.assertCloseAbs(x.real, y.real, eps)
@@ -227,6 +231,43 @@ def test_mul(self):
227231
self.assertRaises(TypeError, operator.mul, 1j, None)
228232
self.assertRaises(TypeError, operator.mul, None, 1j)
229233

234+
self.assertComplexesAreIdentical((1e300+1j) * complex(INF, INF),
235+
complex(NAN, INF))
236+
self.assertComplexesAreIdentical(complex(INF, INF) * (1e300+1j),
237+
complex(NAN, INF))
238+
self.assertComplexesAreIdentical((1e300+1j) * complex(NAN, INF),
239+
complex(-INF, INF))
240+
self.assertComplexesAreIdentical(complex(NAN, INF) * (1e300+1j),
241+
complex(-INF, INF))
242+
self.assertComplexesAreIdentical((1e300+1j) * complex(INF, NAN),
243+
complex(INF, INF))
244+
self.assertComplexesAreIdentical(complex(INF, NAN) * (1e300+1j),
245+
complex(INF, INF))
246+
self.assertComplexesAreIdentical(complex(INF, 1) * complex(NAN, INF),
247+
complex(NAN, INF))
248+
self.assertComplexesAreIdentical(complex(INF, 1) * complex(INF, NAN),
249+
complex(INF, NAN))
250+
self.assertComplexesAreIdentical(complex(NAN, INF) * complex(INF, 1),
251+
complex(NAN, INF))
252+
self.assertComplexesAreIdentical(complex(INF, NAN) * complex(INF, 1),
253+
complex(INF, NAN))
254+
self.assertComplexesAreIdentical(complex(NAN, 1) * complex(1, INF),
255+
complex(-INF, NAN))
256+
self.assertComplexesAreIdentical(complex(1, NAN) * complex(1, INF),
257+
complex(NAN, INF))
258+
259+
self.assertComplexesAreIdentical(complex(1e200, NAN) * complex(1e200, NAN),
260+
complex(INF, NAN))
261+
self.assertComplexesAreIdentical(complex(1e200, NAN) * complex(NAN, 1e200),
262+
complex(NAN, INF))
263+
self.assertComplexesAreIdentical(complex(NAN, 1e200) * complex(1e200, NAN),
264+
complex(NAN, INF))
265+
self.assertComplexesAreIdentical(complex(NAN, 1e200) * complex(NAN, 1e200),
266+
complex(-INF, NAN))
267+
268+
self.assertComplexesAreIdentical(complex(NAN, NAN) * complex(NAN, NAN),
269+
complex(NAN, NAN))
270+
230271
def test_mod(self):
231272
# % is no longer supported on complex numbers
232273
with self.assertRaises(TypeError):
@@ -268,6 +309,7 @@ def test_pow(self):
268309
self.assertAlmostEqual(pow(1j, 200), 1)
269310
self.assertRaises(ValueError, pow, 1+1j, 1+1j, 1+1j)
270311
self.assertRaises(OverflowError, pow, 1e200+1j, 1e200+1j)
312+
self.assertRaises(OverflowError, pow, 1e200+1j, 5)
271313
self.assertRaises(TypeError, pow, 1j, None)
272314
self.assertRaises(TypeError, pow, None, 1j)
273315
self.assertAlmostEqual(pow(1j, 0.5), 0.7071067811865476+0.7071067811865475j)

‎Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
Correct invalid corner cases (resulted in ``(nan+nanj)`` output) in complex
2+
multiplication, e.g. ``(1e300+1j)*(inf+infj)``. Patch by Sergey B
3+
Kirpichev.

‎Objects/complexobject.c

Lines changed: 40 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -53,11 +53,48 @@ _Py_c_neg(Py_complex a)
5353
}
5454

5555
Py_complex
56-
_Py_c_prod(Py_complex a, Py_complex b)
56+
_Py_c_prod(Py_complex z, Py_complex w)
5757
{
5858
Py_complex r;
59-
r.real = a.real*b.real - a.imag*b.imag;
60-
r.imag = a.real*b.imag + a.imag*b.real;
59+
double a = z.real, b = z.imag, c = w.real, d = w.imag;
60+
double ac = a*c, bd = b*d, ad = a*d, bc = b*c;
61+
62+
r.real = ac - bd;
63+
r.imag = ad + bc;
64+
65+
if (isnan(r.real) && isnan(r.imag)) {
66+
/* Recover infinities that computed as (nan+nanj) */
67+
int recalc = 0;
68+
if (isinf(a) || isinf(b)) { /* z is infinite */
69+
/* "Box" the infinity and change nans in the other factor to 0 */
70+
a = copysign(isinf(a) ? 1.0 : 0.0, a);
71+
b = copysign(isinf(b) ? 1.0 : 0.0, b);
72+
if (isnan(c)) c = copysign(0.0, c);
73+
if (isnan(d)) d = copysign(0.0, d);
74+
recalc = 1;
75+
}
76+
if (isinf(c) || isinf(d)) { /* w is infinite */
77+
/* "Box" the infinity and change nans in the other factor to 0 */
78+
c = copysign(isinf(c) ? 1.0 : 0.0, c);
79+
d = copysign(isinf(d) ? 1.0 : 0.0, d);
80+
if (isnan(a)) a = copysign(0.0, a);
81+
if (isnan(b)) b = copysign(0.0, b);
82+
recalc = 1;
83+
}
84+
if (!recalc && (isinf(ac) || isinf(bd) || isinf(ad) || isinf(bc))) {
85+
/* Recover infinities from overflow by changing nans to 0 */
86+
if (isnan(a)) a = copysign(0.0, a);
87+
if (isnan(b)) b = copysign(0.0, b);
88+
if (isnan(c)) c = copysign(0.0, c);
89+
if (isnan(d)) d = copysign(0.0, d);
90+
recalc = 1;
91+
}
92+
if (recalc) {
93+
r.real = Py_INFINITY*(a*c - b*d);
94+
r.imag = Py_INFINITY*(a*d + b*c);
95+
}
96+
}
97+
6198
return r;
6299
}
63100

0 commit comments

Comments
 (0)