Michael Baudin と Robert L. Smith (2011) による "Improved Complex Division in Scilab" 中の compdiv_robust() 関数は、「マクラーレンのの難しい除算のテスト (Mc'Larlen's difficult division test)」でもうまくいきます。 その中で使用されているスケーリングプロセスを使用すると、複素除算のスミス方法(Smith's method for complex division, 1962) もうまく機能しているようです。 つまり、一般的な gcc(GNU Cコンパイラ)の複素除算ではデフォルトでスミスの方法を使用しているので、 このスケーリングプロセスを追加するだけで複素除算を改善することができるということになります。 また、80ビットの拡張精度浮動小数を使用できる場合は、スミスの方法と compdiv_improved() を 64ビットの倍精度浮動小数の範囲で問題なく使用できることがわかりました。
cdiv の計算例
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
「拡張マクラーレンの難しい除算テスト (Extended Mc'Larlen's difficult division test)」の cdiv() のテスト結果抜粋
# 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)
Lambert W function and exponent of Lambert W function for C99