|
投稿者：nu
投稿日：2015年11月 9日(月)13時27分24秒
|
|
|
|
! 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
REM カイ二乗分布の確率分布関数による計算。1000桁モードで実行のこと。数値積分ではないことに注意。ここからENDまで。
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 "chisquaredを小数で入力してください。(例)7.87943857662241735736"
INPUT chisquared
LET SUN=1
LET TUE=1
FOR i=1 TO 1000
IF TUE<10^(-1000) THEN EXIT FOR ! TUE<10^(-1000)の(かっこ内)の累乗の数値を変えて、計算精度を落として、早く終了させることができる。
LET TUE=TUE*chisquared/(n+2*i)
LET SUN=SUN+TUE
NEXT i
LET bb1000=INT(SQR(chisquared/2)^n/GAMMA(n)*EXP(-chisquared/2)*SUN*10^1000+0.5)/10^1000
PRINT bb1000
PRINT #1: bb1000
LET bbketa=INT(bb1000*10^keta+0.5)/10^keta
PRINT bbketa
PRINT #1: bbketa
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
gnuutera2012もといnuもといun
|
|
|