hypotenuse

√(x² + y²) を求める高次収束アルゴリズム

hypotenuse - 直角三角形の斜辺の長さを求める計算と同様の計算が、ベクトルの大きさや複素数の乗除算などで必要になることがしばしばある。多桁の計算を行う場合には、Newton-Raphson や Householder の方法などが使われることが多いが、高次の収束アルゴリズムを用いる場合には、その漸化式に与える初期値が問題になる。

例えば４次の収束アルゴリズムを用いたとしても、半桁程の精度の数を２桁の精度にするために漸化式の計算を１度行うのは非能率的に思える。また、最初に与える値が収束範囲外であれば計算が発散したり、正しい答えが得られないこともある。従って誤差の少ない初期値の与え方が問題と言える。一般にある区間での要求精度が決まっているときには、テーラー級数展開(Taylor series expansion)による近似よりもチェビシェフ近似(Chebyshev approximation)の方が少ない次数で近似できる。そのため、ここでは初期値を２次のチェビシェフ近似で計算しループ内の漸化式に与えている。係数は除数が２のべきになるように調整してある。この近似計算だけで３桁の精度が得られる。

下に示すプログラムは６次の収束を示し、ワーストケースでの誤差の収束の様子は無限精度の計算機で次のようになる。

1.71499e-3 :初期値見積もり時点の誤差
1.04585e-18 :繰り返し１回目の誤差
5.36747e-110 :繰り返し２回目の誤差
9.80773e-658 :繰り返し３回目の誤差

即ち、倍精度演算程度の18桁精度ならループ内の処理を１回行うだけで得られるということである。^*Note1
初期値の見積もり誤差を、１次／１次有理チェビシェフ近似、３次チェビシェフ近似、及び２次／２次のパデ近似(Pade approximation)を行った場合と併せて Fig.1 に示す。もし、荒い見積もりで良ければ、(129 + 53 * x)/128 の近似式でも誤差 1% 以下の初期値が得られる。

Fig.1 sqrt(1+x) の近似誤差

参考プログラムを下に示す。
このプログラムではループ内で係数が簡単で済むため √x の漸近計算を行っているが、多桁の数同士の除算が１箇所必要である。数千桁を越える数値計算の場合には除算の計算速度が問題になるため、1/√x の漸近計算にしたほうが除算を使用しないで済むので望ましい。もしループ内での除算が問題とならないのであれば、ループ内の計算をパデ近似にすることもできる。
余談だが、Moler-Morrison の方法は、√(1 + x) - 1 の１次／１次のパデ近似 2 h / (h + 4) を元に工夫したものと言えよう。

very_long_float
hypotenuse(p, q)
    very_long_float p;
    very_long_float q;
{
    very_long_float r, x, h;

    p = fabs(p);
    q = fabs(q);
    if (p == 0) return(q);
    if (q == 0) return(p);
    if(p < q) {
        r = p; p = q; q = r;
    }
    r = q / p;
    r *= r;
    x = 1 + r * (500 - r * 77) / 1024;  /* 2nd order chebyshev approximation (|err| < 0.09 %) */
    /* Other approximations.
     * x = (129 + r * 53) / 128;                          chebyshev[1] (|err| < 1   %)
     * x = (x * 147 / 214 + 1) / (x * 52 / 269 + 1);    chebyshev[1,1] (|err| < 0.04%)
     * x = ((r + 4) * r * 5 + 16) / ((r + 12) * r + 16);     pade[2,2] (|err| < 0.03%)
     * x = 1 + r * (255 - r * (56 - r * 13)) / 512;       chebyshev[3] (|err| < 0.02%)
     */
    r += 1;
    while(fabs(h = 1 - r / (x * x)) > VLF_EPSILON){
        x *= 1 - h * (128 + h * (32 + h * (16 + h * (10 + h * 7))))/256; /* Taylor series*/
        /* Altanative: Pade[2,2] approximation. (6th order convergence)
         * x *= 1 - h * ((h - 12) * h + 16) / ((6 * h - 32) * h + 32); 
         */
    }
    return(p * x);
}

ループ内で除算を使わない方法(1/√(x+1) の漸化式)。初期値も 1/√(x+1) のチェビシェフ近似で求めている。

Fig.2 1/sqrt(1+x) の近似誤差

収束の様子と誤差の例
-7.5034e-3 :初期値見積もり時点の誤差
8.02596e-14 :繰り返し１回目の誤差
1.20593e-79 :繰り返し２回目の誤差
1.38764e-474 :繰り返し３回目の誤差

very_long_float
hypotenuse(p, q)
    very_long_float p;
    very_long_float q;
{
    very_long_float r, x, h;

    p = fabs(p);
    q = fabs(q);
    if (p == 0) return(q);
    if (q == 0) return(p);
    if(p < q) {
        r = p; p = q; q = r;
    }
    r = q / p;
    r *= r;
    x = (1021 + r * (r * 149 - 444)) / 1024; /* chebyshev[2] (err < 0.38%) */
    r += 1;
    while(fabs(h = 1 - r * x * x) > VLF_EPILON){
        x += x * (h * (128 + h * (96 + h * (80 + h * (70 + h * 63))))/256);
    }
    return(p * x * r);
}

*Note1 :これはものの例えである。このアルゴリズムをそのまま倍精度演算に適用しても計算誤差のために16桁程度の精度しか得られない。ハードウェアで sqrt() がサポートされている場合には、下記アルゴリズムが誤差が少なく(1ULP以下)適している。

#ifndef M_SQRT2
#define M_SQRT2  1.4142135623730950488
#endif

double
_hypot(x, y)
    double x;
    double y;
{
    register double a;
    register double b;

    x = fabs(x);
    y = fabs(y);
    if(x == 0) return(y);
    if(y == 0) return(x);
    if(y > x) {
        a = x; x = y; y = a;
    }

    a = x - y;
    if(a > y){
        a = x / y;
        a += sqrt(1.0 + a * a);
    }else{
        a /= y;
        b = a * (a + 2.0);
        a += b / (M_SQRT2 + sqrt(2.0 + b));
        a += M_SQRT2 + 1.0;
    }
    return(x + y / a);
}

hypotenuse の語源はギリシャ語の hupo (下) + teinousa (伸) だそうです。

Nov. 21 2020 : ソースコード中 fabs()処理位置の誤植修正

REFERENCE:

C. B. Moler and D. Morrison, IBM Journal of Research and Development 27(6), 577-581, 1983