√２ の話：その２１：Babylonian methodとNewton's method

★バビロニアの平方根とニュートン法の関係。
http://yamada-kuebiko.cocolog-nifty.com/blog/2015/02/post-ed0c.html
2015年2月10日 (火)：√２ の話：その２０：バビロニアの平方根
↑前の記事で、バビロニアの平方根の話を書きました。

http://cpplover.blogspot.jp/2010/11/blog-post_20.html
バビロニア人の方法（Babylonian method）
↑このページからの引用文では

しかも、速度もそれほど遅くない。
・・・
たったの6回の計算で、それなりの精度の平方根を求められる。

とありました。

★実はですね、このバビロニアの平方根とニュートン法は数学的に同じものなのです。
そのことを明示的(explicit)に述べたサイトは比較的少ないように思います。「匂わせる」だけでね。
私、ＥＸＰＬＩＣＩＴなのが好きな性質(たち)ですので、明示してしまいましょう。

この図の上の部分は既に使ったものです。
今回は、それに下の部分を付け加えました。
上の図でお分かりと思いますが、f(x) = x^2 - S にニュートン法を適用すると
Xn+1 = (Xn + S/Xn)/2
こうなります。

上の引用サイトで紹介されたアルゴリズムを私が式にしたものを書きました。

このアルゴリズムを式の形に整理すると
Xn+1 = (Xn + S/Xn)/2
こうなります。

こうでした。
ほら、同じ式になるのですよ。この式に到達する経路は異なるのですが、到達点は同じなのですね。
結果として数学的には等価なのです。

★私が最初に書いたニュートン法のプログラムでは
メインプログラムのループ中には
　　x1 = x - (f(x, n)/g(x))
としてニュートン法の原理そのものを書き
必要な関数 f(x) とその導関数 g(x) は外に書きました。
これによって、関数を書き換えるだけで、いろいろな関数にニュートン法を適用できるようになります。
そういう汎用性を考慮した書き方だったのです。

プログラム①
! ニュートン法によって平方根を求める
DECLARE EXTERNAL FUNCTION f
DECLARE EXTERNAL FUNCTION g

INPUT PROMPT "n =" :n
LET x = 1
FOR i = 1 TO 15
   LET x1 = x - (f(x, n)/g(x))
   PRINT i;
   PRINT x1;
   PRINT x1*x1
   LET x = x1
NEXT i
END
!********************
EXTERNAL FUNCTION f(x, n)
LET f = x*x - n
END FUNCTION
!********************
EXTERNAL FUNCTION g(x)
LET g = 2*x
END FUNCTION
!********************

★今回、関数を書き込んでしまって、
Xn+1 = (Xn + S/Xn)/2
このように整理した形でプログラムを書き直してみました。

プログラム②
! バビロニアの平方根

INPUT PROMPT "n =" :n
LET x = 1
FOR i = 1 TO 15
   LET x1 = (x + n/x)/2
   PRINT i;
   PRINT x1;
   PRINT x1*x1
   LET x = x1
NEXT i
END

プログラム自体は簡潔になりましたが、他の関数に対応させようとすると、厄介なことになるでしょう。

さて、実行結果ですが、ループ回数６回までをお目にかけます
実行結果①
n =2
1  1.5  2.25
2  1.41666666666667  2.00694444444445
3  1.41421568627451  2.00000600730488
4  1.41421356237469  2.00000000000451
5  1.4142135623731  2.00000000000001
6  1.4142135623731  2.00000000000001
実行結果②
n =2
1  1.5  2.25
2  1.41666666666667  2.00694444444445
3  1.41421568627451  2.00000600730488
4  1.41421356237469  2.00000000000451
5  1.4142135623731  2.00000000000001
6  1.4142135623731  2.00000000000001

全く同じです。当たり前、数学的に全く同じことをしているのですから。

★

しかも、速度もそれほど遅くない。
・・・
たったの6回の計算で、それなりの精度の平方根を求められる。

惜しいなぁ。「たったの6回の計算で、それなりの精度の平方根を求められる。」どころではないのです。

私のトライアル↓
http://yamada-kuebiko.cocolog-nifty.com/blog/2015/02/post-fc81.html

★１１回のループで、どのくらい正しかったのか知りたくなりました。
・・・
最後のところです。
8215212822 9518488472
8215212822 951848847←この「７」が９９９桁目です。
ここまで正しかったのですね。

６回で終わらずに、更に５回ほど繰り返すとほぼ１０００桁正しくなっちゃうんですねぇ。
２次の収束、というものの恐ろしさは、このくらいやってみないとわかりません。

かつて、ＪＡＶＡの「BigDecimal」というのを使って√２をニュートン法で１０００桁以上求めてみたことがあるのです。
びっくりしましたねぇ。数学科の先生に報告したら、やっぱり２次の収束だものなぁ、と笑っておられましたね。

●ＪＡＶＡによる多倍長計算が載ってます↓興味のある方はどうぞ。
https://oku.edu.mie-u.ac.jp/~okumura/java2/keisan-big.html

http://mail2.nara-edu.ac.jp/~asait/java/big/big.htm