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