SF bug [ #409448 ] Complex division is braindead

http://sourceforge.net/tracker/?func=detail&aid=409448&group_id=5470&atid=105470
Now less braindead.  Also added test_complex.py, which doesn't test much, but
fails without this patch.

diff --git a/Lib/test/output/test_complex b/Lib/test/output/test_complex
new file mode 100644 (file)
index 0000000..b5224a9
--- /dev/null
+++ b/Lib/test/output/test_complex
@@ -0,0 +1 @@
+test_complex

diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py
new file mode 100644 (file)
index 0000000..ef44fbc
--- /dev/null
+++ b/Lib/test/test_complex.py
@@ -0,0 +1,65 @@
+from test_support import TestFailed
+from random import random
+
+# XXX need many, many more tests here.
+
+nerrors = 0
+
+def check_close_real(x, y, eps=1e-12):
+    """Return true iff floats x and y "are close\""""
+    # put the one with larger magnitude second
+    if abs(x) > abs(y):
+        x, y = y, x
+    if y == 0:
+        return abs(x) < eps
+    if x == 0:
+        return abs(y) < eps
+    # check that relative difference < eps
+    return abs((x-y)/y) < eps
+
+def check_close(x, y, eps=1e-12):
+    """Return true iff complexes x and y "are close\""""
+    return check_close_real(x.real, y.real, eps) and \
+           check_close_real(x.imag, y.imag, eps)
+
+def test_div(x, y):
+    """Compute complex z=x*y, and check that z/x==y and z/y==x."""
+    global nerrors
+    z = x * y
+    if x != 0:
+        q = z / x
+        if not check_close(q, y):
+            nerrors += 1
+            print `z`, "/", `x`, "==", `q`, "but expected", `y`
+    if y != 0:
+        q = z / y
+        if not check_close(q, x):
+            nerrors += 1
+            print `z`, "/", `y`, "==", `q`, "but expected", `x`
+
+simple_real = [float(i) for i in range(-5, 6)]
+simple_complex = [complex(x, y) for x in simple_real for y in simple_real]
+for x in simple_complex:
+    for y in simple_complex:
+        test_div(x, y)
+
+# A naive complex division algorithm (such as in 2.0) is very prone to
+# nonsense errors for these (overflows and underflows).
+test_div(complex(1e200, 1e200), 1+0j)
+test_div(complex(1e-200, 1e-200), 1+0j)
+
+# Just for fun.
+for i in range(100):
+    test_div(complex(random(), random()),
+             complex(random(), random()))
+
+try:
+    z = 1.0 / (0+0j)
+except ZeroDivisionError:
+    pass
+else:
+    nerrors += 1
+    raise TestFailed("Division by complex 0 didn't raise ZeroDivisionError")
+
+if nerrors:
+    raise TestFailed("%d tests failed" % nerrors)

diff --git a/Objects/complexobject.c b/Objects/complexobject.c
index 3c1301c7a0c175d2ae4c928d44f61357e37ee36b..34dbab0becaaaf98a6689ff1d177fad4dd202228 100644 (file)
--- a/Objects/complexobject.c
+++ b/Objects/complexobject.c
@@ -29,7 +29,8 @@
 
 static Py_complex c_1 = {1., 0.};
 
-Py_complex c_sum(Py_complex a, Py_complex b)
+Py_complex
+c_sum(Py_complex a, Py_complex b)
 {
        Py_complex r;
        r.real = a.real + b.real;
@@ -37,7 +38,8 @@ Py_complex c_sum(Py_complex a, Py_complex b)
        return r;
 }
 
-Py_complex c_diff(Py_complex a, Py_complex b)
+Py_complex
+c_diff(Py_complex a, Py_complex b)
 {
        Py_complex r;
        r.real = a.real - b.real;
@@ -45,7 +47,8 @@ Py_complex c_diff(Py_complex a, Py_complex b)
        return r;
 }
 
-Py_complex c_neg(Py_complex a)
+Py_complex
+c_neg(Py_complex a)
 {
        Py_complex r;
        r.real = -a.real;
@@ -53,7 +56,8 @@ Py_complex c_neg(Py_complex a)
        return r;
 }
 
-Py_complex c_prod(Py_complex a, Py_complex b)
+Py_complex
+c_prod(Py_complex a, Py_complex b)
 {
        Py_complex r;
        r.real = a.real*b.real - a.imag*b.imag;
@@ -61,8 +65,16 @@ Py_complex c_prod(Py_complex a, Py_complex b)
        return r;
 }
 
-Py_complex c_quot(Py_complex a, Py_complex b)
+Py_complex
+c_quot(Py_complex a, Py_complex b)
 {
+       /******************************************************************
+       This was the original algorithm.  It's grossly prone to spurious
+       overflow and underflow errors.  It also merrily divides by 0 despite
+       checking for that(!).  The code still serves a doc purpose here, as
+       the algorithm following is a simple by-cases transformation of this
+       one:
+
        Py_complex r;
        double d = b.real*b.real + b.imag*b.imag;
        if (d == 0.)
@@ -70,9 +82,45 @@ Py_complex c_quot(Py_complex a, Py_complex b)
        r.real = (a.real*b.real + a.imag*b.imag)/d;
        r.imag = (a.imag*b.real - a.real*b.imag)/d;
        return r;
+       ******************************************************************/
+
+       /* This algorithm is better, and is pretty obvious:  first divide the
+        * numerators and denominator by whichever of {b.real, b.imag} has
+        * larger magnitude.  The earliest reference I found was to CACM
+        * Algorithm 116 (Complex Division, Robert L. Smith, Stanford
+        * University).  As usual, though, we're still ignoring all IEEE
+        * endcases.
+        */
+        Py_complex r;  /* the result */
+        const double abs_breal = b.real < 0 ? -b.real : b.real;
+        const double abs_bimag = b.imag < 0 ? -b.imag : b.imag;
+
+        if (abs_breal >= abs_bimag) {
+               /* divide tops and bottom by b.real */
+               if (abs_breal == 0.0) {
+                       errno = EDOM;
+                       r.real = r.imag = 0.0;
+               }
+               else {
+                       const double ratio = b.imag / b.real;
+                       const double denom = b.real + b.imag * ratio;
+                       r.real = (a.real + a.imag * ratio) / denom;
+                       r.imag = (a.imag - a.real * ratio) / denom;
+               }
+       }
+       else {
+               /* divide tops and bottom by b.imag */
+               const double ratio = b.real / b.imag;
+               const double denom = b.real * ratio + b.imag;
+               assert(b.imag != 0.0);
+               r.real = (a.real * ratio + a.imag) / denom;
+               r.imag = (a.imag * ratio - a.real) / denom;
+       }
+       return r;
 }
 
-Py_complex c_pow(Py_complex a, Py_complex b)
+Py_complex
+c_pow(Py_complex a, Py_complex b)
 {
        Py_complex r;
        double vabs,len,at,phase;
@@ -101,7 +149,8 @@ Py_complex c_pow(Py_complex a, Py_complex b)
        return r;
 }
 
-static Py_complex c_powu(Py_complex x, long n)
+static Py_complex
+c_powu(Py_complex x, long n)
 {
        Py_complex r, p;
        long mask = 1;
@@ -116,7 +165,8 @@ static Py_complex c_powu(Py_complex x, long n)
        return r;
 }
 
-static Py_complex c_powi(Py_complex x, long n)
+static Py_complex
+c_powi(Py_complex x, long n)
 {
        Py_complex cn;