C99用
C99用に書き直した "compdiv_robust" アルゴリズム
及び、そこで使用されているスケーリング処理を使用した Smith の方法(1962)との比較

\frac{p+\mathrm{i}q}{r+\mathrm{i}s}= \frac{pr+qs}{r^2+s^2} + \mathrm{i}\frac{qr - ps}{r^2 + s^2}

2020年11月27日
細田 隆之 (aka Lyuka)
English English edition is here.

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.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 の計算例

# 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)
「拡張マクラーレンの難しい除算テスト (Extended Mc'Larlen's difficult division test)」の cdiv() のテスト結果抜粋

#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

www.finetune.co.jp [Mail] © 2000 Takayuki HOSODA.