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@[...].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@[...].net
>  https://lists.sourceforge.net/lists/listinfo/numpy-discussion
> 


-------------------------------------------------------
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_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@[...].net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion