Complex Division Patch (reprise)

[Stephen L Moshier]

Here is the complex divide method from C9X.

       [#8]  EXAMPLE 2 Division  of  two  double  _Complex operands
       could be implemented as follows.

       #include <math.h>
       #include <complex.h>

       /* Divide z / w ... */
       double complex _Cdivd(double complex z, double complex w)
       {
               #pragma STDC FP_CONTRACT OFF
               double a, b, c, d, logbw, denom, x, y;
               int ilogbw = 0;
               a = creal(z); b = cimag(z);
               c = creal(w); d = cimag(w);
               logbw = logb(fmax(fabs(c), fabs(d)));
               if (isfinite(logbw)) {
                       ilogbw = (int)logbw;
                       c = scalbn(c, -ilogbw);
                       d = scalbn(d, -ilogbw);
               }
               denom = c * c + d * d;
               x = scalbn((a * c + b * d) / denom, -ilogbw);
               y = scalbn((b * c - a * d) / denom, -ilogbw);
               /*
                * Recover infinities and zeros that computed
                * as NaN+iNaN; the only cases are non-zero/zero,
                * infinite/finite, and finite/infinite, ...
                */
               if (isnan(x) && isnan(y)) {
                       if ((denom == 0.0) &&
                               (!isnan(a) || !isnan(b))) {
                               x = copysign(INFINITY, c) * a;
                               y = copysign(INFINITY, c) * b;
                       }
                       else if ((isinf(a) || isinf(b)) &&
                               isfinite(c) && isfinite(d)) {
                               a = copysign(isinf(a) ? 1.0 : 0.0, a);
                               b = copysign(isinf(b) ? 1.0 : 0.0, b);
                               x = INFINITY * ( a * c + b * d );
                               y = INFINITY * ( b * c - a * d );
                       }
                       else if (isinf(logbw) &&
                               isfinite(a) && isfinite(b)) {
                               c = copysign(isinf(c) ? 1.0 : 0.0, c);
                               d = copysign(isinf(d) ? 1.0 : 0.0, d);
                               x = 0.0 * ( a * c + b * d );
                               y = 0.0 * ( b * c - a * d );
                       }
               }
               return x + I * y;
       }

       [#9] Scaling the denominator alleviates  the  main  overflow
       and  underflow  problem,  which  is  more  serious  than for
       multiplication.  In the spirit of the multiplication example
       above,  this  code  does  not  defend  against  overflow and
       underflow in the calculation of the numerator.  Scaling with
       the  scalbn  function,  instead  of  with division, provides
       better roundoff characteristics.