| √Ｎの計算 |
平方根 √A のニュートン法による漸化式
Xn+1 ：= (Xn + A /
Xn) /
2
繰り返す回数は、１００万桁≒２の２０乗 ２０回。
１億桁なら、≒２の２６．５乗→２７回。
A /
Xnは、多倍長÷多倍長の形なので計算が複雑になる。
――――――――――――――――――――――――――――――――――――――――――――――――――
1
/ √A のニュートン法による漸化式
Xn+1： = Xn * (3 - A * Xn2) / 2 <乗算回数 5回 係数含>
これだと回数は多いが 多倍長÷多倍長はなくなる。
※２の割算があるが、０．５倍と考えれば、すべて乗算になる。
√A＝ A ＊ （1 / √A） であるから、
1 /
√Aがわかれば、√Aが 乗算でもとめられる。
――――――――――――――――――――――――――――――――――――――――――――――――――
! √Ｎの計算 百桁区切り 十進ＢＡＳＩＣ
! １／√Ｎの漸化式
! Xn+1： = Xn+1： = Xn * (3 - A * Xn2) / 2
!
OPTION ARITHMETIC
DECIMAL_HIGH
!――――――――入力条件―――――――――――――――
DO
DO
DO
input PROMPT
"N=":N
if int(N)<>N then print "正の整数しか計算できません"
loop until
int(N)=N
if int(sqr(N))=sqr(N) then print "平方数かも知れません"
loop while
int(sqr(N))=sqr(N)
if N>10^6 then print "1000000超えると計算できません"
loop while
N>1000000
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
!―――――――エンジン部分―――――――――――――
print "KS=";KS+1
for K=0
to KS+1
print "計算";K;"回目";TIME$
gosub 100 ! 計算ルーチン本体部分
next
K
!―――――――――――――――――――――――――――
!――――最終処理
√N＝N×（１／√N）―――――――――――
FOR I=0 to BLOCK
LET XN0(I)=N*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 * Xn2) / 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
! ②Ａ*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
! ③３－Ａ *
Xn^2の計算
LET AXN02m3(0)=3
for I=0 TO BLOCK
LET
AXN02m3(I)=AXN02m3(I)-AXN02(I)
next I
! ④Xn*(３ -Ａ* 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
!
⑤0.5*Xn*(３ -Ａ* 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)=N
LET OPT5(0)=0 ! 0.5
LET OPT5(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
OPT5(I)=0
next I
LET XN0(0)=int(1/sqr(N)) ! 1／√Ｎ
LET
XN0(1)=int((1/sqr(N)-int(1/sqr(N)))*1000)*10^(100-3)
LET A(0)=N
LET
OPT5(0)=0 ! 0.5
LET OPT5(1)=5*10^99 !
0.5
return
END