|
|
!ゼータ関数 50桁
OPTION ARITHMETIC DECIMAL_HIGH !1000桁モード
LET format$="##." & REPEAT$("#",49)
FOR i=-1 TO 0.9 STEP .1
LET z=fnZeta(i)
PRINT "zeta(";STR$(i)&")";
PRINT USING format$:z
NEXT i
END
EXTERNAL FUNCTION fnZeta(X)
OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0 !ゼロベース
LET M=200
DIM C(M)
LET C(0)=1
for K=1 to M
FOR J=K TO 1 STEP -1
LET C(J)=(C(J)+C(J-1))/2
NEXT j
NEXT k
for J=1 to M
LET U=U+((-1)^INT(J-1))*C(J)/Lpower(J,X)
! PRINT "\:";INT(J-1)
IF C(J)<1E-53 THEN EXIT FOR
NEXT j
LET fnZeta=U/(1-Lpower(2,(1-X)))
END FUNCTION
EXTERNAL FUNCTION Lpower(n,x)
OPTION ARITHMETIC DECIMAL_HIGH
LET Lpower=LEXP(x*LLOG(N))
END FUNCTION
EXTERNAL FUNCTION LLOG(X)
OPTION ARITHMETIC DECIMAL_HIGH
IF X<=0 THEN
CAUSE EXCEPTION 3004
ELSEIF X<1 THEN
LET LLOG=-LLOG(1/X)
ELSEIF X>3 THEN
LET LLOG=2*LLOG(SQR(X))
ELSE ! 1<=x<=3
LET H=(X-1)/(X+1) ! 0<=h<=0.5
LET T=0
LET N=1
LET K=H
LET H2=H^2
DO
LET T=T+K/N
LET N=N+2
LET K=K*H2
LOOP UNTIL K<=1E-120
LET LLOG=2*T
END IF
END FUNCTION
!1000桁モードで利用する指数関数
EXTERNAL FUNCTION LEXP(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 LEXP=s(1,1)
END FUNCTION
---------------------------
1000桁モードで零点の計算はできないものかと思い、小数点有りのゼータ関数を作成しました。
50桁でもかなり遅いです。少しでも速く成ったら良いなと思います。
http://blogs.yahoo.co.jp/donald_stinger
|
|