√3の計算 |
平方根 √A のニュートン法による漸化式
Xn+1 := (Xn + A /
Xn) /
2
繰り返す回数は、100万桁≒2の20乗 20回。
1億桁なら、≒2の26.5乗→27回。
A /
Xnは、多倍長÷多倍長の形なので計算が複雑になる。
――――――――――――――――――――――――――――――――――――――――――――――――――
1
/ √A のニュートン法による漸化式
Xn+1: = Xn * (3 - A * Xn2) / 2 <乗算回数 5回 係数含>
但し、A=3のときには 1.5*Xn * (1 - Xn2) <乗算回数 2回>
これだと回数は多いが 多倍長÷多倍長はなくなる。
※2の割算があるが、0.5倍と考えれば、すべて乗算になる。
√A= A * (1 / √A) であるから、
1 /
√Aがわかれば、√Aが 乗算でもとめられる。
――――――――――――――――――――――――――――――――――――――――――――――――――
√3 100万桁.txt
(約100kB):但し、抜粋。
計算時間 12時間27分22秒 (10進BASICによる100桁区切り計算にて)
CPU:Celelon 900MHz
――――――――――――――――――――――――――――――――――――――――――――――――――
! √3の計算 百桁区切り
! 1/√3の漸化式
! 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×(1/√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
!
A1*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
! B1−1 *
Xn^2の計算
LET AXN02m3(0)=1
for I=0 TO BLOCK
LET
AXN02m3(I)=AXN02m3(I)-AXN02(I)
next I
! CXn*(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
!
D1.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 ! 桁処理1
--負記号の削除
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 ! 桁処理2 --桁上げ処理
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
!2回目以降の初期化
2回目以降の初期化2回目以降の初期化2回目以降の初期化2回目以降の初期化2回目以降の初期化
for I=0 To BLOCK
LET
A(I)=0
! XN0(I)=0 2回目の初期化では『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