[Numpy-discussion] complex division
Tim Hochberg
tim.hochberg at cox.net
Sun Feb 19 14:34:02 CST 2006
While rummaging around Python's complexobject.c looking for code to
steal for complex power, I came across the following comment relating to
complex division:
/******************************************************************
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.)
errno = EDOM;
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.
*/
The algorithm shown, and maligned, in this comment is pretty much
exactly what is done in numpy at present. The function goes on to use
the improved algorithm, which I will include at the bottom of the post.
It seems nearly certain that using this algorithm will result in some
speed hit, although I'm not certain how much. I will probably try this
out at some point and see what the speed hit, but in case I drop the
ball I thought I'd throw this out there as something we should at least
look at. In most cases, I'll take accuracy over raw speed (within reason).
-tim
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;
More information about the Numpy-discussion
mailing list