diff options
| -rw-r--r-- | Lib/test/output/test_complex | 1 | ||||
| -rw-r--r-- | Lib/test/test_complex.py | 65 | ||||
| -rw-r--r-- | Objects/complexobject.c | 66 |
3 files changed, 124 insertions, 8 deletions
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 @@ 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; |