円周率の計算　〜１万桁程度

　シャープの式　　　　　　　　収束は遅い。

　　π/6=atn(１/√3)

　　 atn(x)=ｘ-(ｘ^3)/3+(x^5)/5-(ｘ^7)/7+・・・
　　　　　　 =ｘ(1-ｘ^2/3+ｘ^4/5-ｘ^6/7+・・・　　　　　　：x=１/√3

　　
　　π=2*√3･(1-1/3･3+1/3^2･5-1/3^3･7・・・)

　　tan(π/a)=b/√c　　a,b,cは正の整数となるものは
　　以下、２つだけ。

　　tan(π/4)=1
　　tan(π/6)=1/√3

　　尚、tan(π/3)=√3は、収束しないので円周率計算には使えない。

!■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
! シャープの式による計算 １００桁処理版
!
!tan(π/6)=1/√3
!atn(x)=ｘ-(ｘ^3)/3+(x^5)/5-(ｘ^7)/7+・・・
! =ｘ(1-ｘ^2/3+ｘ^4/5-ｘ^6/7+・・・
!x=１/√3
!
! π=2*√3･(1-1/3･3+1/3^2･5-1/3^3･7・・・)
!
!■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■

OPTION ARITHMETIC RATIONAL
OPTION base 0

do
INPUT prompt "桁数=":KTA_2
IF KTA_2/1000<>INT(KTA_2/1000) THEN PRINT "1000桁単位で入力してください"
loop until KTA_2/1000=INT(KTA_2/1000)

! 初期設定――――――――――――――――――――――――――
LET BLOKS=100
LET JJJ=10^BLOKS
LET QQ=JJJ-1

LET KTA=KTA_2/BLOKS*1.4 !40％マージン

print "計算開始時刻=";date$&" "&time$
let T1$=time$
!print "桁数=";KTA_2
! ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★

! ――配列を準備する―商・余・Σ・基底――――――――――――――――――――――――
dim SYO(KTA*2),JOYO(KTA*2),Σ(KTA*2),BASE(KTA*2),RT3(KTA*2),Z(KTA*2)
dim XN0(KTA*2),XN1(KTA*2),A(KTA*2),XN02(KTA*2)
dim AXN02(KTA*2),AXN02m3(KTA*2),AXN02m3X(KTA*2)
dim OPT5(KTA*2)
! ――――――――――――――――――――――――――――――――――――――――

GOSUB 10000 !√３の計算

!級数部分の計算
! π=2*√3･(1-1/3･3+1/3^2･5-1/3^3･7・・・)
! ■■■(1-1/3･3+1/3^2･5-1/3^3･7・・・)■■■■
LET BASE(0)=1
LET Σ(0)=1

LET kou=0
DO
LET kou=kou+100 !(2kou-1)・3^kou≒KTA_2桁数 となるkou
! log(2*kou-1)+kou*log3>KTA_2
! kou>KTA_2/log3=KTA_2/0.477
loop until kou>KTA_2/0.477

PRINT "kou=";kou
FOR k=2 to kou
GOSUB 1000 !resetする SYO JOYO

for I =0 to KTA
LET JOYO(I)=BASE(I) !前のBASEをＬｏａｄする
! PRINT "B:";BASE(I)
NEXT I

! 3の割算を先に処理
LET DIV_X=3 !3で割る
GOSUB 600 !SYOに答えが帰ってくる

! '割ったあと ベースを保管
for I =0 to KTA
LET BASE(I)=SYO(I)
LET JOYO(I)=SYO(I)
NEXT I

LET DIV_X=2*k-1 !2*k-1で割る
GOSUB 600 !SYOに答えが帰ってくる
! ◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆

for I =0 to KTA
! PRINT "k=";k
! PRINT "SYO(I)=";SYO(I)

LET Σ(I)=Σ(I)+(-1)^(k-1)*SYO(I)
! PRINT "Σ(I)=";Σ(I)
NEXT I


next k

!桁処理
GOSUB 800 !桁借り
GOSUB 900 !桁上げ


500
FOR I=0 TO KTA
LET RT3(I)=XN0(I)
NEXT I

!√3部分の計算

for I =0 to KTA
FOR J=0 to KTA
if I+J=<KTA then LET Z(I+J)=Σ(I)*2*RT3(J)+Z(I+J)
NEXT J
NEXT I

for I=0 to KTA
LET Σ(I)=Z(I)
NEXT I

GOSUB 900 !桁上げ

!印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷
PRINT "π計算終了時刻=";date$;" ";time$
PRINT ""
print str$(Σ(0));"."
for I=1 TO KTA
LET p$=repeat$("0",100-len(str$(Σ(I))))&str$(Σ(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




600 ! *DIV_X
! ■■■■■■■■■■割算■■■■■■■■■■■■■■■■■■■■■■■■
! *******ソロバン・ルーチン開始*********************
for I= 1 to KTA
! 商を計算
LET SYO(I)=int((JJJ*JOYO(I-1) + JOYO(I)) / DIV_X)
! 余りを算出
LET JOYO(I)=(JJJ*JOYO(I-1) + JOYO(I))-SYO(I)*DIV_X
next I
! ******ソロバン・ルーチン終了！！******************
! ***商の結果を余りへ置換え*****
for I=1 to KTA
LET JOYO(I) =SYO(I)
next I
! ***置換え終了!****************
return


800 ! *KETASYORI
! ■■■■■■■■■負の数になるときの桁借りルーチン■■■■■■■■■■■■■

DO
LET f_ken=0
for K_D=1 TO KTA*2
if Σ(K_D)<0 then
LET Σ(K_D-1)=Σ(K_D-1)-1 !上から借りて
LET Σ(K_D)=Σ(K_D)+JJJ !その桁に足す
LET f_ken=1
end if
NEXT K_D
loop while f_ken=1
return
! ■■■■■■■■■桁借りルーチン２ 終了！■■■■■■■■■■■■■■■■■



! ■桁上げ処理■■■■■■■■■■■■■
900 ! *KEA_UP
FOR G_KTA=KTA to 1 step -1
LET Σ(G_KTA-1)=Σ(G_KTA-1)+int(Σ(G_KTA) / JJJ)
LET Σ(G_KTA)=MOD( Σ(G_KTA) , JJJ )
next G_KTA
return
! ■■■■■■■■■■■■■■


! ■■リセットルーチン■■■■■
1000 ! *RE_SET
for I =0 to KTA
LET SYO(I)=0
LET JOYO(I)=0
next I
return
! ■■■■■■■■■■■■■■■




10000 !■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
! √３の計算 百桁区切り √３の計算 百桁区切り √３の計算 百桁区切り √３の計算 百桁区切り √３の計算 百桁区切り √３の計算 百桁区切り
! １／√３の漸化式
! Xn+1： = 1.5*Xn*(1 - Xn^2)
!

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


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

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

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

! 桁処理 --桁上げ処理
Do
LET OBF=0
for I= KTA 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 "√３の計算終了時刻=";date$;" ";time$
!印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷印刷
return



10100 !
! Xn+1： = 1.5*Xn*(1 - Xn^2)

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

! �A1*Xn^2の計算
! for I=0 TO KTA
! for J=0 TO KTA-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 KTA
LET AXN02m3(I)=AXN02m3(I)-XN02(I) ! 係数１のため �A項の計算を飛ばして
next I

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

! LET XN1(I)=AXN02m3X(I)*1.5
next J
next I

11000 ! 桁処理１ --負記号の削除
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


12000 ! 桁処理２ --桁上げ処理
DO
LET OBF=0
for I= KTA 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


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

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

return


20000 !初期設定 初期設定 初期設定 初期設定 初期設定 初期設定 初期設定 初期設定 初期設定 初期設定
for I=0 To KTA
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(3)) ! 1／√3≒0.577
LET XN0(1)=int((1/sqr(3)-int(1/sqr(3)))*1000000)*10^(100-6)
LET A(0)=1
LET OPT5(0)=1 ! 1.5
LET OPT5(1)=5*10^99 ! 1.5
return

END