[Numpy-discussion] complex division

[Charles R Harris]

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
>