逆数と平方根を求める高次収束アルゴリズム

example of sqrt(2) by Newton-Raphson method

細田　隆之
2009年11月16日
2025年 5月 9日 改定

ニュートン法 (Newton-Raphson methode) は関数 f(x) を、近似値 xn の周りで１次近似したものであるが、 高次の近似により収束速度を高めることができる。

拡張精度や多倍長精度の計算では乗算の計算時間が支配的であるので、 目的とする桁数の計算結果を得るための乗算回数を少なくするために、 高次の収束を示す漸化式を利用するのが望ましい。 ここでは高次収束の方法を逆数や平方根の計算に適用した例を紹介する。


EnglishEnglish edition is here.

逆数

ニュートン法により、逆数 1/A を求める漸化式は次式で与えられる:

xn+1 = 2 xn - A xn2
有効桁は、繰り返しごとに２倍になる。

拡張精度や多倍長精度の計算では乗算の計算時間が支配的であるので、目的とする桁数の計算結果を得るための乗算回数を少なくするために、高次の収束を示す漸化式を利用するのが望ましい。

以下では 1/(1 - h) の級数展開を使用している。

1/A の 3 次収束の漸化式

hn は、xn の誤差とする。 繰り返しごとに有効桁は約３倍になる。

hn = 1 - A xn
xn+1 = xn (1 + hn + hn2)

1/A の 4 次収束の漸化式

hn = 1 - Axn
xn+1 = xn(1 + hn)(1 + hn2)

1/A の 5 次収束の漸化式

hn = 1 - A xn
xn+1 = xn (1 + (1 + hn2) (hn + hn2))

平方根

同様の手法が平方根の計算にも適用できる。平方根 √A のニュートン法による漸化式は次式で与えられる：

xn+1 = (xn + A / xn) / 2

しかし、この式には A/xn の除算を行わなければならないという不利な点がある。
そこで、1/√A のニュートン法を考えてみると、より良い漸化式が得られることに気がつく。

1/√A のニュートン法による漸化式

xn+1 = xn (3 - A xn2) / 2

つまり、3/2, 1/2 は定数であるので、上の式では乗算が必要とされるだけである。
√A を求めるには、B = 1/√A を求めてから、 √A = AB により √A を求める。

逆数の計算の場合と同様に高次収束の漸化式を利用するのが良い。
以下では、1/√(1 - h) の級数展開を使用している。

1/√A の 3 次収束の漸化式

hn は、xn の誤差と同じである。繰り返しごとに有効桁は３倍になる。

hn = 1 - Axn2
xn+1 = xn(1 + hn(4 + 3hn)/8)

1/√A の 4 次収束の漸化式

hn = 1 - A xn2
xn+1 = xn (1 + hn (8 + hn (6 + 5 hn)) / 16)

1/√A の 5 次収束の漸化式

hn = 1 - Axn2
xn+1 = xn(1 + (64hn + hn2(48 + 40hn + 35hn2))/128)

1/√A の 6 次収束の漸化式

hn = 1 - Axn2
xn+1 = xn(1 + hn(1/2 + hn(3/8 + hn(5/16 + hn(35/128 + hn(63/256))))))

その他の 6 次収束の漸化式

√A の 6 次収束の漸化式

1/√A の場合と同様に、√(1 - h) の級数展開を使用している。
hn の計算に除算が必要という重大なデメリットがあるのに要注意。

hn = 1 - A/xn2
xn+1 = xn(1 - hn(1/2 + hn(1/8 + hn(1/16 + hn(5/128 + hn(7/256))))))

1/∛A (= A-1/3) の 6 次収束の漸化式

hn = 1 - A xn3
xn+1 = xn + xn hn (1/3 + hn (2/9 + hn (14/81 + hn (35/243 + hn 91/729))))

1/∜A (= A-1/3) の 6 次収束の漸化式

hn = 1 - A xn4
xn+1 = xn + xn hn (1/4 + hn (5/32 + hn (15/128 + hn (195/2048 + hn 663/8192))))

プログラム例

Listing 1:
6 次収束の平方根の C での実装例
double _sqrt(double x) {
    double a, d, h;
    int    e;

    /* Lookup table for scaling factor based on exponent parity */
    static const double scf[3] = { M_SQRT2, 1.0, M_SQRT1_2 };

    if (x <  0) return NAN;     /* Return NaN for negative input */
    if (x == 0) return x;       /* sqrt(0) = 0 */

    /* Range reduction: x = d * 2^e, where d ∈ [0.5, 1) */
    d = frexp(x, &e);

    /* Initial approximation of 1 / sqrt(d) using 3rd-order polynomial */
    /* Approximation error |err| < 5e-4 over d ∈ [0.5, 1] */
    a = 2.605008 + d * (-3.6382977 + d * (2.9885788 - 0.9557557 * d));

    /* Scaling: convert 1/sqrt(d) to 1/sqrt(x) */
    a = ldexp(a * scf[1 + e % 2], -e / 2);

    /* Compute residual h = 1 - x * a^2, then apply correction */
    /* 6th-order convergence in final refinement */
    h = 1 - x * a * a;
    a *= (1 + h * (0.5 + h * (0.375 + h * (0.3125 + h * (0.2734375 + h * 0.24609375)))));

    return x * a; /* Return sqrt(x) = x * (1 / sqrt(x)) */
}

SEE ALSO

www.finetune.co.jp [Mail]