🚀 go-pugleaf

近岡＠大門高校（富山県）です。

hirayama@etl.go.jp (Fuminori Hirayama) wrote:
>平山＠電総研といいます．
>
>In article <52vvmp$psc@inetsv1.daimon-hs.daimon.toyama.jp>, <chikaoka@daimon-hs.daimon.toyama.jp> wrote:
>
>> 近岡＠大門高校（富山県）です。
>（略）
>> (3)のNewton法による平方根のプログラム例は、
>> 　　（C言語風、ただし、aのn乗をa^nであらわす。）
>（略）
>>   /* Newton法 */
>>   do {
>>     t=s; s=(t^2+x)/(2*t); /* # ここが問題 */ 
>>   } while (s<t);
>> 
（中略）
>>　　プログラムの#印の所は、
>>　　    t=s; s=(t+x/t)/2;
>>　　とすると、高速でオーバーフローの心配も少なくなります。
>
>これを方法(A)として，
>
>> このNewton法はループ１回ごとに有効桁数が２倍になります。
>> しかし、プログラム中の#印の所を
>>     t=s; u=s*s; t=t-2*t*(u-x)/(3*u+x);
>> とすれば、ループ１回ごとに 有効桁数が３倍になるような気がします。
>こっちを(B)とします．

[1] この式のままだと、ほぼ１００％（例外あり）、
有効桁数の最後の１ビット分だけ大きくなります。
（浮動小数点数でも、固定小数点数でもです。）
誤差を無くすには、
  s=t-2*t*(u-x)/(3*u+x);
を、
  s=t*(u+3*x)/(3*u+x);
に修正する必要があります。
このとき、分子・分母のオーバーフローが心配です。

>(B)の繰り返しの式を変形すると
>
>    (t(n+1)-sqrt(x))=(t(n)-sqrt(x))^3/(3*t(n)^2+x)
>
>となって，t(n)がxに十分近いときには１回のiterationで桁数は３倍になる
>と思います．だから同じ桁数（精度）まで計算するためには，(A)に比べて
>2/3の繰り返し回数ですむことになります．でも計算時間を考えると，
>1iterationあたりだいたい
>
>方法(A)    除算１回
>方法(B)    自乗１回＋乗算１回＋除算１回
>
>となるので，結局(A)の方が速くなるのではないかと思います．

[2] 方法(A)３回と、方法(B)２回とでほぼ同じ効果が得られるので、
　方法(A)３回　　除算３回
　方法(B)２回    自乗２回＋乗算２回＋除算２回
で、除算１回が（自乗２回＋乗算２回）より速いかどうかが問題になるのですね。

従って、方法(B)が効果を発揮するのは、
　[2] 多倍長計算などで乗法の４倍以下の速度で除法ができない場合。
　[1] 分子・分母の計算でオーバーフローしないこと。
の２条件を同時にみたす時に限られるのですね。

>それにしても(B)は面白い方法ですね．参考文献等はあるのでしょうか？

ご存じの方はぜひ教えてください。
私自身は、立方根のアルゴリズムでループ１回で有効桁数が３倍になるものが
　奥村『C言語による最新アルゴリズム』（技術評論社）
にのっていたので、平方根でもできないかと考えた結果、
この方法を思いつきました。

近岡宣吉　Chikaoka, NobuYoshi
E-mail:chikaoka@daimon-hs.daimon.toyama.jp