The "compdiv_robust" algorithm rewritten for C99
And a comparison with Smith's method (1962) using the scaling process used in it.



Nov. 27 2020
Takayuki HOSODA (aka Lyuka)

Japanese edition is here.

The compdiv_robust() function in "Improved Complex Division in Scilab" by Michael Baudin and Robert L Smith (2011) works well even at Mc'Larlen's difficult division test. Using the scaling process used in it, Smith's method for complex division (1962) also seems to work well. This means that popular gcc -- GNU C compiler -- uses Smith's method for complex division (1962) by default, so you can improve complex division simply by adding this scaling process. Also, it's found that if 80-bit extended-double can be used, Smith's method (1962) and compdiv_improved() can be used in 64-bit double range without problem.

cdiv.c (C99)

NAME cdiv -- returns complex division. LIBRARY Math library (libm, -lm) SYNOPSIS #include <math.h> #include <complex.h> #include <float.h> double complex cdiv(double complex x, double complex y); RETURN VALUE cdiv(x, y) returns complex division x / y.

/* cdiv() - the "compdiv_robust()" function rewritten for C99 by Takayuki HOSODA (aka Lyuka). Rev. 0.3.4 (Nov. 27 2020) The algolithm of the function is described in the paper "Improved Complex Division in Scilab", Michael Baudin, Robert L. Smith, Jun. 2011 http://forge.scilab.org/index.php/p/compdiv/downloads/get/improved_cdiv_v0.3.pdf */ #include <math.h> #include <complex.h> #include <float.h> double complex cdiv(double complex x, double complex y) { double complex z; double r, t; #ifdef USE_SCALING double AB, CD, S; const double OV2 = DBL_MAX / 2; const double UNBeps = DBL_MIN * 2 / DBL_EPSILON; const double Be = 2 / (DBL_EPSILON * DBL_EPSILON); AB = fmax(fabs(creal(x)), fabs(cimag(x))); CD = fmax(fabs(creal(y)), fabs(cimag(y))); S = 1; /* pre-scaling */ if (AB >= OV2) { x = x / 2; S = S * 2; } // Scale down x if (CD >= OV2) { y = y / 2; S = S / 2; } // Scale down y if (AB < UNBeps) { x = x * Be; S = S / Be;} // Scale up x if (CD < UNBeps) { y = y * Be; S = S * Be;} // Scale up y #endif #ifdef USE_COMPDIV_IMPROVED /* cdiv_improved (modified - Dividing by t instead of multiplying the reciprocal of t seems better from the point of view of accuracy) */ if (cabs(cimag(y)) <= cabs(creal(y))) { r = cimag(y) / creal(y); t = (creal(y) + cimag(y) * r); if (r == 0) { z = (creal(x) + cimag(y) * (cimag(x) / creal(y))) / t + (cimag(x) - cimag(y) * (creal(x) / creal(y))) / t * I; } else { z = (creal(x) + cimag(x) * r) / t + (cimag(x) - creal(x) * r) / t * I; } } else { r = creal(y) / cimag(y); t = (creal(y) * r + cimag(y)); if (r == 0) { z = (creal(y) * (creal(x) / cimag(y)) + cimag(x)) / t + (creal(y) * (cimag(x) / cimag(y)) - creal(x)) / t * I; } else { z = (creal(x) * r + cimag(x)) / t + (cimag(x) * r - creal(x)) / t * I; } } #else #ifdef USE_SMITH /* Smith's method (1962) */ if (cabs(creal(y)) < cabs(cimag(y))) { r = creal(y) / cimag(y); t = (r * creal(y)) + cimag(y); z = ((r * creal(x) + cimag(x)) / t + (r * cimag(x) - creal(x)) / t * I); } else { r = cimag(y) / creal(y); t = creal(y) + (r * cimag(y)); z = ((creal(x) + r * cimag(x)) / t + (cimag(x) - r * creal(x)) / t * I); } #else /* Use compiler's complex division. */ /* Note that gcc with -std=-gnu89 (default) works fine, but -std=c99 does not. */ z = x / y; #endif #endif #ifdef USE_SCALING return S * z; // Post-scaling #else return z; #endif } #ifdef DEBUG #include <stdio.h> double complex cdiv(double complex x, double complex y); int main() { double complex c, x, y, z; double g, k; #ifdef USE_COMPDIV_IMPROVED #ifdef USE_SCALING char * myname = "compdiv_robust"; #else char * myname = "compdiv_improved"; #endif #else #ifdef USE_SMITH #ifdef USE_SCALING char * myname = "cdiv_smith_robust"; #else char * myname = "cdiv_smith"; #endif #else #if __STDC_VERSION__ >= 199901L char * myname = "cdiv_C99"; #else char * myname = "cdiv_cc"; #endif #endif #endif #if __FLT_EVAL_METHOD__ > 0 printf("cdiv: EXCESS_PRECISION_STANDARD enabled.\n"); #endif printf("cdiv = %s((1 + I) * g, (1 + I * k) * g), g = DBL_MAX / 2\n", myname); g = DBL_MAX / 2; for (k = 0.0; k <= 2.0; k += 0.5) { x = 1 + I; y = 1 + I * k ; c = cdiv(x, y); #ifdef USE_C99 z = (x * g) / (y * g); #else z = cdiv(x * g, y * g); #endif printf("k=%2.1f, Ref.=(%.5g, %.5g), %s->(%.5g, %.5g), error=(%.5g, %.5g)\n", k, creal(c), cimag(c), myname, creal(z), cimag(z), creal(c) - creal(z), cimag(c) - cimag(z)); } return 0; } #endif

# cdiv test results. # g = DBL_MAX / 2; z = cdiv((1 + I) * g, (1 + I * k) * g); # cdiv_smith -- Smith's method (1962) % gcc -O test_cdiv.c -DDEBUG -DUSE_SMITH -lm ; ./a.out cdiv = cdiv_smith((1 + I) * g, (1 + I * k) * g), g = DBL_MAX / 2 k=0.0, Ref.=(1, 1), cdiv_smith->(1, 1), error=(0, 0) k=0.5, Ref.=(1.2, 0.4), cdiv_smith->(1.2, 0.4), error=(0, 0) k=1.0, Ref.=(1, 0), cdiv_smith->(1, 0), error=(0, 0) k=1.5, Ref.=(0.76923, -0.15385), cdiv_smith->(0, -0), error=(0.76923, -0.15385) k=2.0, Ref.=(0.6, -0.2), cdiv_smith->(0, -0), error=(0.6, -0.2) # cdiv_smith -- Smith's method (1962) using 80-bit extended double % gcc -O -mfpmath=387 test_cdiv.c -DDEBUG -DUSE_SMITH -lm ; ./a.out cdiv: EXCESS_PRECISION_STANDARD enabled. cdiv = cdiv_smith((1 + I) * g, (1 + I * k) * g), g = DBL_MAX / 2 k=0.0, Ref.=(1, 1), cdiv_smith->(1, 1), error=(0, 0) k=0.5, Ref.=(1.2, 0.4), cdiv_smith->(1.2, 0.4), error=(0, 0) k=1.0, Ref.=(1, 0), cdiv_smith->(1, 0), error=(0, 0) k=1.5, Ref.=(0.76923, -0.15385), cdiv_smith->(0.76923, -0.15385), error=(0, -2.7756e-17) k=2.0, Ref.=(0.6, -0.2), cdiv_smith->(0.6, -0.2), error=(0, 0) # cdiv_smith_robust -- Smith's method (1962) with pre-post-scaling % gcc -O test_cdiv.c -DDEBUG -DUSE_SCALING -DUSE_SMITH -lm ; ./a.out cdiv = cdiv_smith_robust((1 + I) * g, (1 + I * k) * g), g = DBL_MAX / 2 k=0.0, Ref.=(1, 1), cdiv_smith_robust->(1, 1), error=(0, 0) k=0.5, Ref.=(1.2, 0.4), cdiv_smith_robust->(1.2, 0.4), error=(0, 0) k=1.0, Ref.=(1, 0), cdiv_smith_robust->(1, 0), error=(0, 0) k=1.5, Ref.=(0.76923, -0.15385), cdiv_smith_robust->(0.76923, -0.15385), error=(0, -5.5511e-17) k=2.0, Ref.=(0.6, -0.2), cdiv_smith_robust->(0.6, -0.2), error=(0, 0) # compdiv_improved (2011) % gcc -O test_cdiv.c -DDEBUG -DUSE_COMPDIV_IMPROVED -lm ; ./a.out cdiv = compdiv_improved((1 + I) * g, (1 + I * k) * g), g = DBL_MAX / 2 k=0.0, Ref.=(1, 1), compdiv_improved->(1, 1), error=(0, 0) k=0.5, Ref.=(1.2, 0.4), compdiv_improved->(1.2, 0.4), error=(0, 0) k=1.0, Ref.=(1, 0), compdiv_improved->(1, 0), error=(0, 0) k=1.5, Ref.=(0.76923, -0.15385), compdiv_improved->(0, -0), error=(0.76923, -0.15385) k=2.0, Ref.=(0.6, -0.2), compdiv_improved->(0, -0), error=(0.6, -0.2) # compdiv_improved (2011) using 80-bit extended double % gcc -O -mfpmath=387 test_cdiv.c -DDEBUG -DUSE_COMPDIV_IMPROVED -lm ; ./a.out cdiv: EXCESS_PRECISION_STANDARD enabled. cdiv = compdiv_improved((1 + I) * g, (1 + I * k) * g), g = DBL_MAX / 2 k=0.0, Ref.=(1, 1), compdiv_improved->(1, 1), error=(0, 0) k=0.5, Ref.=(1.2, 0.4), compdiv_improved->(1.2, 0.4), error=(0, 0) k=1.0, Ref.=(1, 0), compdiv_improved->(1, 0), error=(0, 0) k=1.5, Ref.=(0.76923, -0.15385), compdiv_improved->(0.76923, -0.15385), error=(0, -2.7756e-17) k=2.0, Ref.=(0.6, -0.2), compdiv_improved->(0.6, -0.2), error=(0, 0) # compdiv_robust (2011) -- compdiv_improved (2011) with pre-post-scaling % gcc -O test_cdiv.c -DDEBUG -DUSE_SCALING -DUSE_COMPDIV_IMPROVED -lm ; ./a.out cdiv = compdiv_robust((1 + I) * g, (1 + I * k) * g), g = DBL_MAX / 2 k=0.0, Ref.=(1, 1), compdiv_robust->(1, 1), error=(0, 0) k=0.5, Ref.=(1.2, 0.4), compdiv_robust->(1.2, 0.4), error=(0, 0) k=1.0, Ref.=(1, 0), compdiv_robust->(1, 0), error=(0, 0) k=1.5, Ref.=(0.76923, -0.15385), compdiv_robust->(0.76923, -0.15385), error=(0, -5.5511e-17) k=2.0, Ref.=(0.6, -0.2), compdiv_robust->(0.6, -0.2), error=(0, 0) # C99 complex division % gcc -O test_cdiv.c -DDEBUG -DUSE_C99 -std=c99 -lm ; ./a.out cdiv = cdiv_C99((1 + I) * g, (1 + I * k) * g), g = DBL_MAX / 2 k=0.0, Ref.=(1, 1), cdiv_C99->(1, 1), error=(0, 0) k=0.5, Ref.=(1.2, 0.4), cdiv_C99->(inf, 0.4), error=(-inf, 5.5511e-17) k=1.0, Ref.=(1, 0), cdiv_C99->(inf, 0), error=(-inf, 0) k=1.5, Ref.=(0.76923, -0.15385), cdiv_C99->(inf, -0.15385), error=(-inf, -2.7756e-17) k=2.0, Ref.=(0.6, -0.2), cdiv_C99->(inf, -0.2), error=(-inf, -2.7756e-17)

#1: g = DBL_MAX / k; x = k + I; y = 1 + I; z = compdiv_robust(x * g, y * g) #1, k=2^0, Ref.=(1, 0), compdiv_robust->(1, 0) #1, k=2^1, Ref.=(1.5, -0.5), compdiv_robust->(1.5, -0.5) #1, k=2^26, Ref.=(33554432.5, -33554431.5), compdiv_robust->(33554432.5, -33554431.5) #1, k=2^27, Ref.=(67108864.5, -67108863.5), compdiv_robust->(67108864.5, -67108863.5) #1, k=2^52, Ref.=(2251799813685248, -2251799813685248), compdiv_robust->(2251799813685248, -2251799813685248) #1, k=2^53, Ref.=(4503599627370496, -4503599627370496), compdiv_robust->(4503599627370497, -4503599627370496) #1: g = DBL_MAX / k; x = k + I; y = 1 + I; z = cdiv_smith_robust(x * g, y * g) #1, k=2^0, Ref.=(1, 0), cdiv_smith_robust->(1, 0) #1, k=2^1, Ref.=(1.5, -0.5), cdiv_smith_robust->(1.5, -0.5) #1, k=2^26, Ref.=(33554432.5, -33554431.5), cdiv_smith_robust->(33554432.5, -33554431.5) #1, k=2^27, Ref.=(67108864.5, -67108863.5), cdiv_smith_robust->(67108864.5, -67108863.5) #1, k=2^52, Ref.=(2251799813685248, -2251799813685248), cdiv_smith_robust->(2251799813685248, -2251799813685248) #1, k=2^53, Ref.=(4503599627370496, -4503599627370496), cdiv_smith_robust->(4503599627370497, -4503599627370496) #2: g = DBL_MAX / k; x = 1 + k * I; y = 1 + I; z = compdiv_robust(x * g, y * g) #2, k=2^0, Ref.=(1, 0), compdiv_robust->(1, 0) #2, k=2^1, Ref.=(1.5, 0.5), compdiv_robust->(1.5, 0.5) #2, k=2^26, Ref.=(33554432.5, 33554431.5), compdiv_robust->(33554432.5, 33554431.5) #2, k=2^27, Ref.=(67108864.5, 67108863.5), compdiv_robust->(67108864.5, 67108863.5) #2, k=2^52, Ref.=(2251799813685248, 2251799813685248), compdiv_robust->(2251799813685248, 2251799813685248) #2, k=2^53, Ref.=(4503599627370496, 4503599627370496), compdiv_robust->(4503599627370497, 4503599627370496) #2: g = DBL_MAX / k; x = 1 + k * I; y = 1 + I; z = cdiv_smith_robust(x * g, y * g) #2, k=2^0, Ref.=(1, 0), cdiv_smith_robust->(1, 0) #2, k=2^1, Ref.=(1.5, 0.5), cdiv_smith_robust->(1.5, 0.5) #2, k=2^26, Ref.=(33554432.5, 33554431.5), cdiv_smith_robust->(33554432.5, 33554431.5) #2, k=2^27, Ref.=(67108864.5, 67108863.5), cdiv_smith_robust->(67108864.5, 67108863.5) #2, k=2^52, Ref.=(2251799813685248, 2251799813685248), cdiv_smith_robust->(2251799813685248, 2251799813685248) #2, k=2^53, Ref.=(4503599627370496, 4503599627370496), cdiv_smith_robust->(4503599627370497, 4503599627370496) #3: g = DBL_MAX / k; x = 1 + I; y = k + I; z = compdiv_robust(x * g, y * g) #3, k=2^0, Ref.=(1, 0), compdiv_robust->(1, 0) #3, k=2^1, Ref.=(0.6, 0.2), compdiv_robust->(0.6, 0.2) #3, k=2^26, Ref.=(1.490116141589226e-08, 1.490116097180305e-08), compdiv_robust->(1.490116141589226e-08, 1.490116097180305e-08) #3, k=2^27, Ref.=(7.450580652434979e-09, 7.450580541412677e-09), compdiv_robust->(7.450580652434979e-09, 7.450580541412677e-09) #3, k=2^52, Ref.=(2.220446049250314e-16, 2.220446049250313e-16), compdiv_robust->(2.220446049250314e-16, 2.220446049250313e-16) #3, k=2^53, Ref.=(1.110223024625157e-16, 1.110223024625156e-16), compdiv_robust->(1.110223024625157e-16, 1.110223024625156e-16) #3: g = DBL_MAX / k; x = 1 + I; y = k + I; z = cdiv_smith_robust(x * g, y * g) #3, k=2^0, Ref.=(1, 0), cdiv_smith_robust->(1, 0) #3, k=2^1, Ref.=(0.6, 0.2), cdiv_smith_robust->(0.6, 0.2) #3, k=2^26, Ref.=(1.490116141589226e-08, 1.490116097180305e-08), cdiv_smith_robust->(1.490116141589226e-08, 1.490116097180305e-08) #3, k=2^27, Ref.=(7.450580652434979e-09, 7.450580541412677e-09), cdiv_smith_robust->(7.450580652434979e-09, 7.450580541412677e-09) #3, k=2^52, Ref.=(2.220446049250314e-16, 2.220446049250313e-16), cdiv_smith_robust->(2.220446049250314e-16, 2.220446049250313e-16) #3, k=2^53, Ref.=(1.110223024625157e-16, 1.110223024625156e-16), cdiv_smith_robust->(1.110223024625157e-16, 1.110223024625156e-16) #4: g = DBL_MAX / k; x = 1 + I; y = 1 + k * I; z = compdiv_robust(x * g, y * g) #4, k=2^0, Ref.=(1, 0), compdiv_robust->(1, 0) #4, k=2^1, Ref.=(0.6, -0.2), compdiv_robust->(0.6, -0.2) #4, k=2^26, Ref.=(1.490116141589226e-08, -1.490116097180305e-08), compdiv_robust->(1.490116141589226e-08, -1.490116097180305e-08) #4, k=2^27, Ref.=(7.450580652434979e-09, -7.450580541412677e-09), compdiv_robust->(7.450580652434979e-09, -7.450580541412677e-09) #4, k=2^52, Ref.=(2.220446049250314e-16, -2.220446049250313e-16), compdiv_robust->(2.220446049250314e-16, -2.220446049250313e-16) #4, k=2^53, Ref.=(1.110223024625157e-16, -1.110223024625156e-16), compdiv_robust->(1.110223024625157e-16, -1.110223024625156e-16) #4: g = DBL_MAX / k; x = 1 + I; y = 1 + k * I; z = cdiv_smith_robust(x * g, y * g) #4, k=2^0, Ref.=(1, 0), cdiv_smith_robust->(1, 0) #4, k=2^1, Ref.=(0.6, -0.2), cdiv_smith_robust->(0.6, -0.2) #4, k=2^26, Ref.=(1.490116141589226e-08, -1.490116097180305e-08), cdiv_smith_robust->(1.490116141589226e-08, -1.490116097180305e-08) #4, k=2^27, Ref.=(7.450580652434979e-09, -7.450580541412677e-09), cdiv_smith_robust->(7.450580652434979e-09, -7.450580541412677e-09) #4, k=2^52, Ref.=(2.220446049250314e-16, -2.220446049250313e-16), cdiv_smith_robust->(2.220446049250314e-16, -2.220446049250313e-16) #4, k=2^53, Ref.=(1.110223024625157e-16, -1.110223024625156e-16), cdiv_smith_robust->(1.110223024625157e-16, -1.110223024625156e-16)

1. "Improved Complex Division in Scilab", Michael Baudin, Robert L.Smigh, Jun. 2011
2. "Complex division", Robert L. Smith, 1962

Lambert W function and exponent of Lambert W function for C99