逆数と平方根を求める高次収束アルゴリズム
細田 隆之
2009年11月16日
example of sqrt(2) by Newton-Raphson method

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

English English edition is here.

逆数

ニュートン法により、逆数 1 / A を求める漸化式は
Xn+1 = 2 * Xn - A * Xn2
で与えられる。有効桁は、繰り返しごとに2倍になる。多桁の計算では乗算の計算時間が支配的であるので、目的とする桁数の計算結果を得るための乗算回数を少なくするために高次の収束をする漸化式を利用するのが望ましい。 ここでは 1 / (1 - h) の級数展開を使用している。

1 / A の3次収束の漸化式
hn = 1 - A * xn
xn+1 = xn * (1 + hn + hn2)

繰り返しごとに有効桁は3倍になる。hn は、xn の誤差と同じである。
1 / A の4次収束の漸化式
hn = 1 - A * xn
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 = A * B により √A を求める。
逆数の計算の場合と同様に高次収束の漸化式を利用するのが良い。ここでは、1 / √(1 - h) の級数展開を使用している。

1 / √A の 3次収束の漸化式
hn = 1 - A * xn2
xn+1 = xn * (1 + hn * (4 + 3 * hn) / 8)

hn は、xn の誤差と同じである。繰り返しごとに有効桁は3倍になる。
1 / √A の 4次収束の漸化式
hn = 1 - A * xn2
xn+1 = xn * (1 + hn * (8 + hn * (6 + 5 * hn)) / 16)
1 / √A の 5次収束の漸化式
hn = 1 - A * xn2
xn+1 = xn * (1 + (64 * hn + hn2 * (48 + 40 * hn + 35 * hn2)) / 128)
1 / √A の 6次収束の漸化式
hn = 1 - A * xn2
xn+1 = xn * (1 + hn * (1/2 + hn * (3/8 + hn * (5/16 + hn * (35/128 + hn * 63/256)))))

1 / √A の場合と同様に、√(1 - h) の級数展開を使用している。

√A の6次収束の漸化式
hn = 1 - A / xn2
xn+1 = xn * (1 - hn * (1/2 + hn * (1/8 + hn * (1/16 + hn * (5/128 + hn * 7/256)))))

1 / 3√A の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 / 4√A の6次収束の漸化式
hn = 1 - A * xn4
xn+1 = xn + xn * hn * (1/4 + hn * (5/32 + hn * (15/128 + hn * (195/2048 + hn * 663/8192))))

4次収束の √x のCプログラム例
#include <stdio.h>
#include <math.h>
#ifndef DBL_EPSILON
#define DBL_EPSILON     2.2204460492503131E-16
#endif

double
_sqrtinv(a)
    double a;
{
    double x, h, g;
    int    e;

    /* enhalf the exponent for a half digit initial accuracy */
    frexp(a, &e);
    x = ldexp(1.0, -e >> 1);

    /* 1/sqrt(a) 4th order convergence */
    g = 1.0;
    while(fabs( h = 1.0 - a * x * x) < fabs(g)) {
        x += x * (h * (8.0 + h * (6.0 + 5.0 * h)) / 16.0);
        g = h;
    }
    return(x);
}

double
sqrt(a)
    double(a);
{
    if(a < 0){
        fprintf(stderr, "sqrt: domain error\n");
        return(0);
    }
    return(a * _sqrtinv(a));
}
*Note: 今時のプロセッサではループの条件判断をしている間に中身の計算くらい出来るので、 ループを展開した方が良い。機械ε が10-16くらいならば3回分で十分である。
SEE ALSO: 算術幾何平均による円周率πの計算公式整数の平方根と√(x2 + y2) を求めるアルゴリズム, √(x2+y2)を求める高次収束アルゴリズム

www.finetune.co.jp [Mail]