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.

Lib/test/output/test_complex

@@ -0,0 +1 @@
+test_complex

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)

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;