円周率.jp > Newton 法

Newton 法

多倍長での計算で除算や平方根の計算を行う際に用いられる Newton 法, ならびに Newton-Raphson 法を紹介する.

これらは本来,方程式 f (x)=0 の解を求める方法であるが,適当な f (x) を設定することで逆数や平方根といった値を求めることができる.

Newton 法

まずは Newton 法について説明する. なお,単に「ニュートン法」と言う場合に Newton-Raphson 法を指す場合もある,というより Newton-Raphson 法の方を指す方が多いので注意されたい.

Newton 法の流れ
図 1: Newton 法の流れ

Newton 法では図 1 にあるように yf(x) 上の適当な点 (a0f (a0)) からスタートし, (anf (an)) を通る傾き c の直線と x 軸の交点座標を求める事を繰り返して解を求める. 図 1 から分かる通り,定数 c の選び方によって収束の速さ,収束半径などが違ってくる.

数式を用いてアルゴリズムを記述すると

0 = c(an+1an) + f (an)
an+1an\frac{1}{c} f (an)
によって an+1 を求めていき,十分な精度で収束していったところで値を確定させる.

この中で,定数 c,初期値 a0 は任意に設定できるが, 多倍長計算を用いて逆数や平方根への適用を考える場合は, 初期値は 1/xsqrt(x) で得られる値をそのまま使うと便利である.

Newton-Raphson 法

Newton-Raphson 法の流れ
図 2: Newton-Raphson 法の流れ

Newton-Raphson 法は,Newton 法における定数 cf (x) の導関数 f' (an) で置き換えたもの, つまり図 1 の青線が常に曲線 yf (x) の接線になるようにしたもの(図 2)である.

別の見方をすれば,f (x) を xan で 1 次まで Taylor 展開近似させて解を得る方法でもある. とりあえず数式として表すと

an+1an\frac{f(a_n)}{f'(a_n)}
となる.

この方法は2次近似として知られる. つまり f (x)=0 の根を r とすると

|ran+1| ≦ C|ran|2

となる. 勿論 C |ran| < 1 でなければ収束の保証はない.

適用例

このNewton-Raphson 法を用いて, 冒頭に挙げた除算と平方根の近似値を求める漸化式を求める.

除算

除算は,正確には逆数の近似と乗算を組み合わせる事によって求める. なので Newton-Raphson 法を適用させて求めるのは逆数となる.

多倍長数 AB について B÷A を求めようとするとき, x=1/A が根になっている f (x) として

f (x) = \frac{1}{x}A

を定義する.これにNewton-Raphson 法を適用させることで漸化式

an+1an (2 − Aan)

を得ることができる. 倍精度浮動小数点数などの CPU 命令で容易に得られる逆数の近似値を初期値として用い, 必要な精度で上記の漸化式を必要な回数繰り返すことで目的の精度の逆数が得られる.

繰り返す際の誤差を計算すると

an+1\frac{1}{A} = −A \Bigl( an\frac{1}{A} \Bigr) 2

となり,2次近似であることが確認できる.

平方根

平方根も除算と同様に, 逆平方根の近似と乗算を組み合わせる事によって求める.

多倍長数 A について \sqrt{A} を求めようとするとき,

f (x) = \frac{1}{x%5e2}A
an+1\frac{1}{2} an (3−Aan2) = an\frac{1}{2} an (1−Aan2)

により \frac{1}{\sqrt{A}} を得, A × \frac{1}{\sqrt{A}}\sqrt{A} として平方根を得る.

f (x)=x2A とすることで直接 \sqrt{A} を求めることもできるが, その場合には演算コストの高い多倍長除算が必要となるため,平方根の逆数を求め,元の数との乗算により平方根を求める. 上記の計算において \frac{1}{2} との乗算は O(n) で行うことができるので,多倍長除算は必要ないことが分かる.

計算例

逆数の計算例として,7-1 (=0.14285714…) を計算する. 結果を知っている分,どの程度正確なのかが分かりやすいだろう.

本来コンピュータで計算をする場合は初期値も有効数字 16 桁程度の値を使うが,ここでは 0.1 を初期値とする. また,漸化式も含めた全体のアルゴリズムをプログラム風に書くと

A     := 7                 # 逆数を求める値
x     := 0.1               # 初期値設定
delta := 1 - x * A         # 誤差計算
while( |delta| > eps ) {   # 精度チェック
    x     := x + x * delta
    delta := 1 - x * A     # 誤差計算
}

と,これだけのことである.この時の x の値を計算すると 表 1 のようになる.

表 1.Newton 法の計算例
ループ x delta
0 回目 0.1 0.3
1 回目 0.13 0.09
2 回目 0.1417 0.0081
3 回目 0.14284077 0.00011461
4 回目 0.1428571409806497 0.0000000131354521
5 回目 0.14285714285714283249427116122937 0.00000000000000017254010187139441

また,平方根の計算例として 1/\sqrt{2} (=0.707106…)を計算する.

初期値は有効数字 1 桁分の計算として 0.7 とする. 全体のアルゴリズムは

A     := 2                 # 平方根を求める値
x     := 0.7               # 初期値設定
delta := 1 - x2 * A        # 誤差計算
while( |delta| > eps ) {   # 精度チェック
    x     := x + x * delta / 2
    delta := 1 - x2 * A    # 誤差計算
}

となる.この時の各変数の値は表 2 のようになる.

表 2.Newton 法での平方根の計算例
ループ x delta
0 回目 0.7 0.02
1 回目 0.707 0.000302
2 回目 0.707106757 0.000000068409885902
3 回目 0.707106781186546283451619907 0.000000000000003509934446881563618293673456195202622702

ここの例で各値は計算している桁全てを表示しているが, 実際にはそれほど必要ない. 具体的には 3 回目のループで得られた結果と, そのうち正しい部分だけに切り詰めた数とのそれぞれの 2 乗を計算してみると

0.7071067811865462834516199072 0.499999999999998245032776559218190853163271902398688649
0.707106781186542 0.4999999999999894

と,計算する桁数は倍以上違うのに 2 乗した結果の誤差には 1 桁しか違いがない, この様なことは Newton-Raphson 法の演算中でも言えることなので, 必要無い桁は省くことで無駄な計算を削減することができる. 各計算中に必要な桁数の見積もりと, それに基づいた計算量の削減については [JT04] が詳しい.