[ruby-list:1325] Re: Rational and Complex

> >いえ、結局私は assembler もできないので、Ｃで書いてました。(^^;
> >お題は「x < 2^32 の整数 x の sqrt を求めよ」というもので、

私はassembler会議入ってなかったのでこのゲームは知りませんでしたが
かって MS-DOS2.1 の上の N88-Basic V1.0 の機械語サブルーチンの中で
sqrt を計算する必要が出て、やはりループ無しで作った事があります。
初期値の精度を上げるのも大変だったので４回黙ってやった。

> >    a = 143257 - 202860544/(1582+(x>>22));
> 
>      ^ この数なんだろう??

Newton方で開平すると、一回毎に精度が倍になるのです。
ですから 
> >    a = (a + x/a)/2;
> >    b = (a + x/a)/2;

の２回で片付けるために、2^30 <= x < 2^32 の範囲の全てに
対して最悪 目的精度の 1/4 の有効桁精度を保証する近似値を
用意したのでしょうね。
  ( 1 〜 2^2 ) ＊ 2^30 
の近似値だから、 sqrt( 1 〜 2^2 ) * 2^15 って事で。

> すね. これって, Bignumにも適応できますかね?? つまり, 数の範囲の制限を
> なくせます??

上で述べた私の経験では 近似値の初期値として重要なのは 桁数 なのです。
sqrt(x) の初期値を x そのもので始めても答えはでるのですが、
初期値として桁が合っていないから有効桁数が一桁もない(ことが普通)
ので、何回近似を行うかが設計できないから、結局 誤差 を見ながら
終了判定 という事になります。
初期値の桁が合ってると計算できる。

私がやった時は 初代ナノピコ教室じゃないから判り安くして

 X * 2^2n  と  X * 2^(2n+1)  にわけ、( 1<= X < 2 )

前者( X*2^2n)    のとき  (1 + x/2) * 2^n  を 		X=1+x
後者( X*2(2n+1)) のとき  (1.5 + x/2) * 2^n  を

初期値にしてました。
	アセンブラですから 2^2n 2^(2n+1) の切り分けは 7bit目を見りゃいいし、
	2^2n から 2^n へは指数部を符号付で右シフト
	 (1 + x/2) ってのは 仮数部を 1bit 右シフトですむ
		(１ の部分はメモリ上にないから)
	(1.5 + x/2) は 仮数部を右シフトして 1bit たてる。
これで精度は最悪 ２進でほぼ２桁。４回やると 32bit 精度になる
で、仮数部 24bit だから十分。

(1.5 + x/2)、 (1 + x/2) に当たる部分をもっと精度良くするのが
assembler会議での競争だったのでしょう。
あの 初期値の出し方はよーわからん。
	おせーて ＞ shin 

Bignum って整数ですか？ 無限精度の実数？
Bignum のその時の桁数に応じた初期値の有効桁数の最悪値が
決定できるのならループの終了は誤差でなく、回数でできますね。
最悪値 と 最良値 の幅が大きいと無駄なループも出るけど

a = 143257 - 202860544/(1582+(x>>22)); そのままは 桁数が
ある範囲って事だから使えませんね。
＃ 桁に付いての合わせ込は簡単そうだな。
＃ ん？ 2^22 未満だとどうするのかな。