３体問題

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 ：但し、抜粋。　
計算時間　１２時間２７分２２秒　(10進ＢＡＳＩＣによる１００桁区切り計算にて)
ＣＰＵ：Celelon　９００ＭＨz　

1.
7320508075 6887729352 7446341505 8723669428 0525381038 0628055806 9794519330 1690880003 7081146186 7572485756
7562614141 5406703029 9699450949 9895247881 1655512094 3736485280 9323190230 5582067974 8201010846 7492326501
5312343266 9033228866 5067225466 8921837971 2270471316 6036786158 8019049986 5373798593 8946765034 7506576050
7566183481 2960610094 7602187190 3250831458 2952395983 2997789824 5082887144 6383291734 7224163984 5878553976
6795806381 8353666110 8431737808 9437831610 2088305524 9016700235 2071114428 8695990956 3657970871 6849807289
9493296484 2830207864 0860398873 8697537582 3173178313 9599298300 7838702877 0539133695 6331210370 7264019249
1067682311 9928837564 1141422016 7427521023 7299427083 1059898459 4759876642 8889779614 7837958390 2288548529
0357603385 2808064381 9723446610 5968972287 2865264153 8226646984 2002119548 4155278441 1812865345 0703519165
0016689294 4154808460 7127714399 9762926834 6295774383 6189511012 7148638746 9765459824 5178855097 5379013880
6649619119 6222295711 0555242923 7231921977 3826256163 1468842032 8537166829 3864961191 7049738836 3954959381
1000 桁目

! √３の計算 百桁区切り     十進ＢＡＳＩＣ
! １／√３の漸化式
! Xn+1： = 1.5*Xn*(1 - Xn^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 OPT5(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

!最終計算 √3＝3×（１／√3）
FOR I=0 to BLOCK
LET XN0(I)=3*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： = 1.5*Xn*(1 - Xn^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

! ②1*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

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

! ④Xn*(1 - 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

! ⑤1.5*Xn*(1 - Xn^2)の計算
for I=0 TO BLOCK
for J=0 TO BLOCK-I
LET XN1(I+J)=AXN02m3X(I)*OPT5(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 OPT5(I)=0
next I

LET A(0)=1
LET OPT5(0)=1 ! 1.5
LET OPT5(1)=5*10^99 ! 1.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 OPT5(I)=0
next I

LET XN0(0)=0 ! 0.6が始まり 1／√3≒0.577
LET XN0(1)=6*10^99
LET A(0)=1
LET OPT5(0)=1 ! 1.5
LET OPT5(1)=5*10^99 ! 1.5
return

END