平方根 ニュートン法について

質問

こんにちは。趣味でプログラミングをしているものです。
さっそくですが、質問させていただきます。

Javaで文具店などで１０００円くらいで市販されている
電卓を再現しようとしているのですが、
質問させていただきたいのは、
平方根の計算についてです。
--------------------------------------------
まずは以下のメソッドのコードを見ていただきたいのですが・・・

このメソッドを用いて３の平方根を求めてみました。

private static BigDecimal sqrt( BigDecimal value )
{
BigDecimal two = new BigDecimal( "2" );
BigDecimal x = value.divide( two );
BigDecimal last_x = BigDecimal.ZERO;
BigDecimal gap = x.subtract( last_x );

BigDecimal range = BigDecimal.ONE.movePointLeft( 11 );//(a)
//int cnt = 0;//(b)
BigDecimal t;
while( gap.compareTo(range) > 0 )//(c)
//while( cnt < 20 )//(d)
{
last_x = new BigDecimal( x.toString() );
t = value.divide( x, 64, BigDecimal.ROUND_DOWN );
x = x.add( t ).divide( two );
gap = x.subtract( last_x );//(e)
//cnt ++; //(f)
}
return x;
}

(a)(c)(e) を用いて実行すると
1.732142857142857142857・・・
と小数部分の142857 が循環する値になりました。

(a)(c)(e) をコメントアウトして
(b)(d)(f) を用いて実行すると
1.732050807568877293527446341505872366942805253
810380628055806979425806427001953125
とwikipedia の「３の平方根」の記事の値と
記述されている範囲では同じ値が得られました。
--------------------------------------------
そこで質問なのですが、
(1)
(a) (c) (e) を用いた場合には、正しい値を得られないのでしょうか？
条件に設定する値などを変えてもだめなのでしょうか？
(2)
自分が作っている電卓では１２桁の表示を予定しているのですが
その場合、メソッド sqrt ないの cnt は while ループ内で
いくつまで、インクリメントすればよいのでしょうか？
(3)
インターネットで調べたこのアルゴリズムは
「ニュートン法」だそうですが、
y = x^2 - C の グラフ、接線などを書いてみて
ある程度理解できたのですが、

http://cpplover.blogspot.jp/2010/11/blog-post_20.html

上記サイトではこのアルゴリズムは
バビロニア人の方法（Babylonian method）
というものだそうですが、
バビロニア人は二次関数の微分はわかっていた、
ということでしょうか？
--------------------------------------------
ご存知のかた、教えていただけないでしょうか？
よろしくお願いします。

回答No.2

hirotn · ベストアンサー

示したtoString()による数字の表示例は、BigDecimalクラスのリファレンスより引用したものであります。
BigDecimalは、数字と小数点の有効桁数（スケール）を情報として持ちます。
#1の回答時は、一度toStringで渡しているため、スケールの情報も正しく渡るのか否か疑問がありました。クラスリファレンスの例のように帰ってしまうと２～３桁くらいで丸まってしまわないか、と。

混乱させてしまい申し訳ありません。

-----

再現して確認してみましたが、(a)(c)(e)の場合、

条件判定に使うgapの作り方
gap = x.subtract( last_x );//(e)

ループ終了判定
while( gap.compareTo(range) > 0 )//(c)

この条件がよくなさそうです。この条件によって、ループを2回しか回っていないので結果の精度が良くない。(b)(d)(f)の場合でも、ループを2回実行した場合の結果は、(a)(c)(e)と同様でした。

（１）条件判定に用いるgapの作り方
gap = x.subtract( last_x );//(e)

Xn (プログラム上変数x) と、 Xn-1 (同じくlast_x) の大小関係が保証できないので（一回目のxに依存する）、last_xとxの変化量は絶対値をとった方が良いでしょう。理由はrangeは常にプラスだからです。差分gapがマイナスの場合に不具合が出ます。これが今回、(c)においてループが2回で終わってしまった直接の原因です。

（２）ループ終了判定

(b)(d)(f)では、ループを必ず20回回るので、その分平方根の精度は良くなります。

確認方法は cnt ++; //(f) の下の行に以下のprintを追加して確認しました(toString()は省略できる）。
System.out.println(cnt + " Scale: lastx: " + last_x.scale() + "value: "+ last_x);
System.out.println(cnt + " Scale: t: " + t.scale() + "value: "+ t);
System.out.println(cnt + " Scale: x: " + x.scale() + "value: " + x);
System.out.println(cnt + " Scale: gap: " + gap.scale() + "value: " + gap);

ループ20回目のgapは、20 Scale: gap: 81 value: -2.5806427001953125E-65
すなわち小数点以下81桁まで有効で、-2.58*10^(-64)乗。rangeは10^(-11)乗なので、全然小さい。

（１）（２）より、このように変えてみました。

BigDecimal range = new BigDecimal(1.0e-64);//(a)

while (gap.abs().compareTo(range) > 0) //(c')

結果：
7Scale: lastx: 67value: 1.732050807568877293527446341505872366942805253810380628
0558069794125
7Scale: t: 64value: 1.732050807568877293527446341505872366942805253810380628
0558069794
7Scale: x: 68value: 1.732050807568877293527446341505872366942805253810380628
05580697940625
7Scale: gap: 68value: -6.25E-66

sqrt(3) = 1.73205080756887729352744634150587236694280525381038062805580697940625

ループ7回で抜けてくれました。

回答No.3

f272

(1）
正しい値が求められないのはアルゴリズムが間違っているせいです。
gap = x.subtract( last_x );//(e)
を
gap = x.subtract( last_x ).abs();//(e)
にしましょう。

(2)
それはvalueの値によります。
cntの値で制御するよりも、しっかりと必要な精度を満たしているかどうかを確認すべきでしょう。

(3)
微分という概念はわからなくても、このようなアルゴリズムは発見できるでしょう。
√Cを求めるということはC=(a+b)^2となるようなa+bを求めることですが
a^2がＣに近ければbは非常に小さな数になります。このとき
C=(a+b)^2=a^2+2ab+b^2=a^2+2ab　（b^2は非常に小さな数だから）
となって
b=(C-a^2)(2a)=C/(2a)-a/2
で
√C=a+b=a/2+C/(2a)=(a+C/a)/2
ですね。
微分を知らなくても大丈夫です。

お礼

回答していただき、ありがとうございます。
(1)
絶対値をとるのですか。
変数の内容を確認するべきでした。
わからなかったかもしれないですが・・・
(2)
(a)(c)(e) の組み合わせで正しいアルゴリズムを用いる
べきなのですね。

貴重な時間をさいての回答ありがとうございました。

回答No.1

hirotn

last_x = new BigDecimal( x.toString() );(**)
→last_x.clone(x);
(a)(c)(e)の場合、ここで精度が落ちていないか気になります。
BigDecimal.toString()が返す値は、例えば、以下の様な値だからです。

[123,0] "123"
[-123,0] "-123"
[123,-1] "1.23E+3"
[123,-3] "1.23E+5"
[123,1] "12.3"
[123,5] "0.00123"
[123,10] "1.23E-8"
[-123,12] "-1.23E-10"

逆に、(b)(d)(f)の方法の場合は、(**)で精度を壊さずに(e)でlast_x を使っているので問題ないのかと推測。

お礼

さっそく回答していただきありがとうございます。
解説していただけたところによると、
toString() で返される文字列は
IEEE754 などでの、仮数部、指数部の情報を
持っているものであり、
仮数部の表示桁数に制限があるための
丸め誤差が生じる
ということでよろしいでしょうか？
おそれいりますが、もう少し教えていただけないでしょうか？
よろしくお願いします。