gh-156695: Improve accuracy of complex powers with small negative integer exponents

Author: Aniketsy

Objects/complexobject.c

@@ -357,11 +357,29 @@ c_powu(Py_complex x, long n)
 static Py_complex
 c_powi(Py_complex x, long n)
 {
-    if (n >= 0)
+    if (n > 0)
         return c_powu(x,n);
-    else
-        return c_powu(_Py_c_quot(c_1, x), -n);
 
+    Py_complex r = _Py_c_quot(c_1, c_powu(x, -n));
+    if (errno == EDOM
+        || (isfinite(r.real) && isfinite(r.imag)
+            && (r.real != 0.0 || r.imag != 0.0)))
+        return r;
+
+    /* gh-156695: x**|n| left the exponent range although the result is
+       representable.  Redo it with x scaled to exponent zero; both the
+       scaling and its undoing are exact. */
+    double m = fabs(x.real) > fabs(x.imag) ? fabs(x.real) : fabs(x.imag);
+    if (m == 0.0 || !isfinite(m))
+        return r;
+
+    int e;
+    frexp(m, &e);
+    Py_complex w = {ldexp(x.real, -e), ldexp(x.imag, -e)};
+    r = _Py_c_quot(c_1, c_powu(w, -n));
+    r.real = ldexp(r.real, (int)(e * n));
+    r.imag = ldexp(r.imag, (int)(e * n));
+    return r;
 }
 
 double