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
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
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 3c1301c..34dbab0 100644
--- 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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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;