[Numpy-discussion] complex division
Charles R Harris
charlesr.harris at gmail.com
Sun Feb 19 16:13:02 EST 2006
Hmm...
The new algorithm does look better with respect to overflow and
underflow, but I wonder if it is not a bit of overkill. It seems to me
that the same underflow/overflow problems attend complex
multiplication, which is pretty much all that goes on in the original
algorithm.
One thing I do know is that division is expensive. I wonder if one
division and two multiplications might be cheaper than two divisions.
I'll have to check that out.
Chuck
On 2/19/06, Tim Hochberg <tim.hochberg at cox.net> wrote:
>
> 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;
>
>
>
>
> -------------------------------------------------------
> This SF.net email is sponsored by: Splunk Inc. Do you grep through log files
> for problems? Stop! Download the new AJAX search engine that makes
> searching your log files as easy as surfing the web. DOWNLOAD SPLUNK!
> http://sel.as-us.falkag.net/sel?cmd=lnk&kid=103432&bid=230486&dat=121642
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>
More information about the NumPy-Discussion
mailing list