TSSRC

　　覚え書２

　　　√２　の計算　

　　　平方根 √A のニュートン法による漸化式

　　　X_n+1 ：= (X_n + A / X_n) / 2

　　　これを希望する桁になるまで繰り返し計算する。

　　　繰り返す回数は、１００万桁≒２の２０乗　２０回。
　　　１億桁なら、≒２の２６．５乗→２７回。

　　　A / X_nは、多倍長÷多倍長の形なので計算が複雑になる。
――――――――――――――――――――――――――――――――――――――――――――――――――
　　　1 / √A のニュートン法による漸化式

X_n+1： = X_n * (3 - A * X_n²) / 2 　　　　　<乗算回数　5回　係数含>

但し、Ａ＝２のときには X_n * (１．５ - X_n²)　　<乗算回数　２回>

　　　これだと回数は多いが　多倍長÷多倍長はなくなる。
　　　※２の割算があるが、０．５倍と考えれば、すべて乗算になる。
　　　

　　　 √A＝　A 　＊　（1 / √A）　　であるから、
　　　1 / √Aがわかれば、√Aが　乗算でもとめられる。
――――――――――――――――――――――――――――――――――――――――――――――――――

　　　　√２　１００万桁.txt (約100kＢ）：但し、抜粋。　
　　　　計算時間　１１時間４０分５秒　(10進ＢＡＳＩＣによる１００桁区切り計算にて)
　　　　ＣＰＵ：Celelon　９００ＭＨz　

――――――――――――――――――――――――――――――――――――――――――――――――――

! √２の計算百桁区切り
! １／√２： Xn+1 = Xn * (3 - A * Xn^2) / 2
!
OPTION ARITHMETIC DECIMAL_HIGH
DO
DO
input PROMPT "桁数は？ (1000単位で)":KTA
loop while KTA<2000
loop until KTA/1000=int(KTA/1000)

print "計算開始 ";date$;" ";time$
print "桁数=";KTA;"桁"
! 100桁区切り
LET BLOCK=int(KTA*1.1/100) ! マージン10％とする
print "BLOCK=";BLOCK
option base 0
dim XN0(BLOCK),XN1(BLOCK),A(BLOCK),XN02(BLOCK)
dim AXN02(BLOCK),AXN02m3(BLOCK),AXN02m3X(BLOCK)
dim ZPT5(BLOCK)

LET KS=0
gosub 20000 ! まず初期化

10 !
DO
LET KS=KS+1
loop while KTA > 2^(KS-1)

print "KS=";KS+1
for K=0 to KS+1
print "計算";K;"回目";TIME$
gosub 100
next K

!最終計算 √２＝２×（１／√２）
FOR I=0 to BLOCK
LET XN0(I)=2*XN0(I)
next I

20 ! 桁処理 --桁上げ処理
Do
LET OBF=0
for I=BLOCK TO 1 step -1
LET XN0(I-1)=XN0(I-1)+int(XN0(I)/10^100)
LET XN0(I)=XN0(I)-10^100*int(XN0(I)/10^100)
if XN0(I)>10^100 then LET OBF=1
NEXT I
loop while OBF=1

!印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷
print "終了！"
print "計算終了 ";date$;" ";time$

print str$(XN0(0));"."
for I=1 TO BLOCK
LET p$=repeat$("0",100-len(str$(XN0(I))))&str$(XN0(I))
print MID$(p$,1,10);" ";MID$(p$,11,10);" ";MID$(p$,21,10);" ";MID$(p$,31,10);" ";MID$(p$,41,10);" ";
print MID$(p$,51,10);" ";MID$(p$,61,10);" ";MID$(p$,71,10);" ";MID$(p$,81,10);" ";MID$(p$,91,10);" "
if int(I/10)=I/10 then print I*100;"桁目"

next I
!印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷
stop

100 !
! Xn+1 = Xn * (3 - A * Xn^2) / 2

! ①Xn^2の計算
for I=0 TO BLOCK
for J=0 TO BLOCK-I
LET XN02(I+J)=XN0(I)*XN0(J)+XN02(I+J)
next J
next I

! ②A * Xn^2の計算
for I=0 TO BLOCK
for J=0 TO BLOCK-I
LET AXN02(I+J)=A(I)*XN02(J)+AXN02(I+J)
next J
next I

! ③３－A * Xn^2の計算
LET AXN02m3(0)=3
for I=0 TO BLOCK
LET AXN02m3(I)=AXN02m3(I)-AXN02(I)
next I

! ④Xn *（３－A * Xn^2）の計算
for I=0 TO BLOCK
for J=0 TO BLOCK-I
LET AXN02m3X(I+J)=XN0(I)*AXN02m3(J)+AXN02m3X(I+J)
next J
next I

! ⑤Xn *（３－A * Xn^2）／２の計算
for I=0 TO BLOCK
for J=0 TO BLOCK-I
LET XN1(I+J)=AXN02m3X(I)*ZPT5(J)+XN1(I+J)
next J
next I

1000 ! 桁処理１ --負記号の削除
DO
LET MNSF=0
if XN1(1)<0 then
LET XN1(0)=XN1(0)+int(XN1(1)/10^100)-1
LET XN1(1)=XN1(1)-10^100*int(XN1(1)/10^100)+10^100
end if
for I=BLOCK TO 1 step -1
if XN1(I)<0 then
LET MNSF=1
LET XN1(I-1)=XN1(I-1)+int(XN1(I)/10^100)-1
LET XN1(I)=XN1(I)-10^100*int(XN1(I)/10^100)+10^100
end if
NEXT I
LOOP while MNSF=1

2000 ! 桁処理２ --桁上げ処理
DO
LET OBF=0
for I=BLOCK TO 1 step -1
LET XN1(I-1)=XN1(I-1)+int(XN1(I)/10^100)
LET XN1(I)=XN1(I)-10^100*int(XN1(I)/10^100)
if XN1(I)>10^100 then
LET OBF=1
end if
NEXT I
loop while OBF=1

3000 ! XN1 => XN0 へ入替
for I=0 to BLOCK
LET XN0(I)=XN1(I)
NEXT I

!２回目以降の初期化
for I=0 To BLOCK
LET A(I)=0
! XN0(I)=0 ２回目の初期化では『0』にしない。
LET XN1(I)=0
LET XN02(I)=0
LET AXN02(I)=0
LET AXN02m3(I)=0
LET AXN02m3X(I)=0
LET ZPT5(I)=0
next I

LET A(0)=2
LET ZPT5(1)=5*10^99 ! 0.5

return

20000 !初期設定
for I=0 To BLOCK
LET A(I)=0
LET XN0(I)=0
LET XN1(I)=0
LET XN02(I)=0
LET AXN02(I)=0
LET AXN02m3(I)=0
LET AXN02m3X(I)=0
LET ZPT5(I)=0
next I
LET XN0(0)=1
LET A(0)=2
LET ZPT5(1)=5*10^99 ! 0.5
return

END