[Numpy-discussion] complex division

[Tim H. <tim...@co...>]

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;