|
投稿者：nu
投稿日：2015年11月 9日(月)13時33分52秒
|
|
|
|
! 1000桁モードで利用する指数関数
DECLARE EXTERNAL FUNCTION EXP
REM カイ二乗分布の確率分布関数による計算。数値積分ではないことに注意。
!マチン（Machin）の公式 π/4=4*ArcTan(1/5)-ArcTan(1/239)
OPTION ARITHMETIC DECIMAL_HIGH !1000桁モード
LET p=0
!テイラー展開より、ArcTan(x)=x-x^3/3+x^5/5-x^7/7+ …
!第n項 {16/5^(2n-1)-4/239^(2n-1)}/(2n-1) 符号はnが奇数なら正
LET k=1
LET t=16/5
DO
LET last=p
LET p=p+t/k
LET t=t/(-5*5)
LET k=k+2
LOOP WHILE p<>last
PRINT (k-1)/2 !繰り返し回数 debug
LET k=1
LET t=4/239
DO
LET last=p
LET p=p-t/k
LET t=t/(-239*239)
LET k=k+2
LOOP WHILE p<>last
PRINT (k-1)/2 !繰り返し回数 debug
OPEN #1:NAME "D:\zzzzzzzzzzz.txt"
ERASE #1
LET a1=0
LET a2=0
LET b1=0
LET b2=400
PRINT "nを半角英数で入力してください。(例)2"
INPUT n
FUNCTION GAMMA(w)
IF MOD(w,2)=1 THEN GOTO 100
IF MOD(w,2)=0 THEN GOTO 200
100 LET V=1
FOR I=w TO 1 STEP -2
LET V=I/2*V
NEXT I
LET GAMMA=V*SQR(p)
GOTO 300
200 LET V=1
FOR I=w/2 TO 1 STEP -1
LET V=I*V
NEXT I
LET GAMMA=V
GOTO 300
300 END FUNCTION
PRINT GAMMA(n)
PRINT "小数点以下の桁数を半角英数で入力してください。(例)0.123456789→9"
INPUT keta
PRINT "alphaを小数で入力してください。(例)50%→0.5,75%→0.75,100%→1"
INPUT alpha
REM カイ二乗分布の確率分布関数による計算。1000桁モードで実行のこと。数値積分ではないことに注意。ここからENDまで。
DO
LET bm1=(b1+b2)/2
LET SUN=1
LET TUE=1
FOR J=1 TO 1000
IF TUE<10^(-1000) THEN EXIT FOR ! TUE<10^(-1000)の(かっこ内)の累乗の数値を変えて、計算精度を落として、早く終了させることができる。 ! TUE<10^(-12)の(かっこ内)の累乗の数値を変えて、計算精度を落として、早く終了させることができる。
LET TUE=TUE*bm1/(n+2*J)
LET SUN=SUN+TUE
NEXT J
LET bm1sm11000=INT(SQR(bm1/2)^n/GAMMA(n)*EXP(-bm1/2)*SUN*10^1000+0.5)/10^1000
IF bm1=a1 OR bm1=a2 OR bm1=am1 OR bm1=b1 OR bm1=b2 THEN EXIT DO
IF bm1sm11000<alpha THEN LET b1=bm1
IF alpha=<bm1sm11000 THEN LET b2=bm1
PRINT b1 ! 画面に表示する演算経過PRINT b1を消した方が数倍早くなる。
PRINT #1: b1 !残しておくと演算経過をD:\に残せて良い場合もあるが、場合によってはPRINT #1: b1を消した方が早くなる。
LET b11000=INT(b1*10^keta+0.5)/10^keta
LET b21000=INT(b2*10^keta+0.5)/10^keta
IF b11000=b21000 THEN LET bb=b11000
IF b11000=b21000 THEN EXIT DO
LOOP
PRINT bb
PRINT #1: bb
CLOSE #1
END
EXTERNAL FUNCTION EXP(x)
OPTION ARITHMETIC DECIMAL_HIGH
FUNCTION s(y,n)
LET t=y*x/n
IF ABS(t)<=EPS(0) THEN
LET s=y+t
ELSE
LET s=y+s(t,n+1)
END IF
END FUNCTION
LET EXP=s(1,1)
END FUNCTION
連続投稿すみません。もしかしたら読みにくい所があるかも知れないので、http://yutorinonatuyasumi.blog.fc2.com/blog-entry-123.html
に173.txt, 174.txtを置いておりますので、ダウンロードしてください。
gnuutera2012もといnuもといun
|
|
|