gh-117999: Fix small integer powers of complex numbers #124243
Conversation
* Fix the sign of zero components in the result. E.g. complex(1,-0.0)**2 now evaluates to complex(1,-0.0) instead of complex(1,-0.0). * Fix negative small integer powers of infinite complex numbers. E.g. complex(inf)**-1 now evaluates to complex(0,-0.0) instead of complex(nan,nan). * Powers of infinite numbers no longer raise OverflowError. E.g. complex(inf)**1 now evaluates to complex(inf) and complex(inf)**0.5 now evaluates to complex(inf,nan).
|
I'll take look, but probably already tomorrow :( Meanwhile, I'm not sure if fix for division should be included here. That does make sense for a separate pr, that will provide mixed-mode arithmetic for complex. I was thinking (In fact, I did, patch just lacks more tests: skirpichev#4) on splitting out from skirpichev#1 just cases for real op complex, as specified by C99+ (note, that there is no special case for |
There was a problem hiding this comment.
That looks close to #118000. I doubt there is an easy way to fix powi() algorithm to match generic one (same signs, etc).
|
In mathematics, the non-integer power of a complex number is a multi-value function. We can choose the main value, but this is arbitrary. The integer power of a complex number is always defined unambiguously. It is also more useful. So it is expected that the integer power can be calculated with higher precision and produce less exceptions or NaNs that the general real or complex power. |
No, this will introduce a huge incompatibility, not just with the C standard.
I doubt you can avoid nan component for |
There was a problem hiding this comment.
It turns to be same as #118000, modulo changes in _Py_ADJUST_ERANGE2. I stuck with examples like z=1-1j, but probably we can't do anything in such cases.
Now it looks good for me.
Few remarks:
- See suggestion.
- _Py_ADJUST_ERANGE2 change looks incomplete for me. I think proper fix should include special case for real/integer exponents, just like we did for arithmetic (I'll open an issue.) But if you wish to keep this, lets add tests.
| assert(n > 0); | ||
| while ((n & 1) == 0) { | ||
| x = _Py_c_prod(x, x); | ||
| n >>= 1; | ||
| } | ||
| Py_complex r = x; | ||
| n >>= 1; | ||
| while (n) { | ||
| x = _Py_c_prod(x, x); | ||
| if (n & 1) { | ||
| r = _Py_c_prod(r, x); | ||
| } | ||
| n >>= 1; |
There was a problem hiding this comment.
For me, your tests works fine with my patch from #118000. Here it is, with slightly adjusted formatting (and using new _Py_rc_quot):
diff --git a/Objects/complexobject.c b/Objects/complexobject.c
index b66ebe131a..9dc5d22b37 100644
--- a/Objects/complexobject.c
+++ b/Objects/complexobject.c
@@ -333,19 +333,26 @@ _Py_c_pow(Py_complex a, Py_complex b)
r.real = len*cos(phase);
r.imag = len*sin(phase);
- _Py_ADJUST_ERANGE2(r.real, r.imag);
+ if (isfinite(a.real) && isfinite(a.imag)
+ && isfinite(b.real) && isfinite(b.imag))
+ {
+ _Py_ADJUST_ERANGE2(r.real, r.imag);
+ }
}
return r;
}
+#define INT_EXP_CUTOFF 100
+
static Py_complex
c_powu(Py_complex x, long n)
{
- Py_complex r, p;
+ Py_complex p = x, r = n-- ? p : c_1;
long mask = 1;
- r = c_1;
- p = x;
- while (mask > 0 && n >= mask) {
+ assert(-1 <= n);
+ assert(n < INT_EXP_CUTOFF);
+ while (n >= mask) {
+ assert(mask > 0);
if (n & mask)
r = _Py_c_prod(r,p);
mask <<= 1;
@@ -357,11 +364,10 @@ c_powu(Py_complex x, long n)
static Py_complex
c_powi(Py_complex x, long n)
{
- if (n > 0)
- return c_powu(x,n);
+ if (n >= 0)
+ return c_powu(x, n);
else
- return _Py_c_quot(c_1, c_powu(x,-n));
-
+ return _Py_rc_quot(1.0, c_powu(x, -n));
}
double
@@ -751,9 +757,13 @@ complex_pow(PyObject *v, PyObject *w, PyObject *z)
errno = 0;
// Check whether the exponent has a small integer value, and if so use
// a faster and more accurate algorithm.
- if (b.imag == 0.0 && b.real == floor(b.real) && fabs(b.real) <= 100.0) {
+ if (b.imag == 0.0 && b.real == floor(b.real)
+ && fabs(b.real) <= INT_EXP_CUTOFF)
+ {
p = c_powi(a, (long)b.real);
- _Py_ADJUST_ERANGE2(p.real, p.imag);
+ if (isfinite(a.real) && isfinite(a.imag)) {
+ _Py_ADJUST_ERANGE2(p.real, p.imag);
+ }
}
else {
p = _Py_c_pow(a, b);(asserts are outcome of review, I'm not sure they are useful here.)
What do you think?
|
CC @mdickinson |
Uh oh!
There was an error while loading. Please reload this page.