円周率の計算(検算用)

 投稿者:しばっち  投稿日:2012年 7月11日(水)21時33分43秒
  OPTION ARITHMETIC NATIVE
PUBLIC NUMERIC BIAS, KETA, SIGN, EPS
LET  BIAS = 0
LET  KETA = 10000
LET  SIGN = -BIAS - 1
LET  EPS = 0
LET  KETA = KETA + EPS
DIM X(-BIAS - 1 TO KETA),S(-BIAS-1 TO KETA)
!'π/4=2805*ATN(1/5257)-398*ATN(1/9466)+1950*ATN(1/12943)+1850*ATN(1/34208)+2021*ATN(1/44179)+2097*ATN(1/85353)+1484*ATN(1/114669)+1389*ATN(1/330182)+808*ATN(1/485298)
PRINT (KETA-EPS)*4;"桁の計算開始 ";TIME$
LET T=TIME
PRINT "4*2805*ATN(1/5257)の計算 ";TIME$
CALL SATN(X,5257)
CALL SMUL(X,2805*4)
CALL LCOPY(S,X)
PRINT "4*398*ATN(1/9466)の計算 ";TIME$
CALL SATN(X,9466)
CALL SMUL(X,398*4)
CALL LSUB2(S,X)
PRINT "4*1950*ATN(1/12943)の計算 ";TIME$
CALL SATN(X,12943)
CALL SMUL(X,1950*4)
CALL LADD2(S,X)
PRINT "4*1850*ATN(1/34208)の計算 ";TIME$
CALL SATN(X,34208)
CALL SMUL(X,1850*4)
CALL LADD2(S,X)
PRINT "4*2021*ATN(1/44179)の計算 ";TIME$
CALL SATN(X,44179)
CALL SMUL(X,2021*4)
CALL LADD2(S,X)
PRINT "4*2097*ATN(1/85353)の計算 ";TIME$
CALL SATN(X,85353)
CALL SMUL(X,2097*4)
CALL LADD2(S,X)
PRINT "4*1484*ATN(1/114669)の計算 ";TIME$
CALL SATN(X,114669)
CALL SMUL(X,1484*4)
CALL LADD2(S,X)
PRINT "4*1389*ATN(1/330182)の計算 ";TIME$
CALL SATN(X,330182)
CALL SMUL(X,1389*4)
CALL LADD2(S,X)
PRINT "4*808*ATN(1/485298)の計算 ";TIME$
CALL SATN(X,485298)
CALL SMUL(X,808*4)
CALL LADD2(S,X)
PRINT "計算終了 ";TIME$;TIME-T;"秒"
CALL DISPLAY(S,"PI_検算.TXT")
END

EXTERNAL  SUB SATN(S(),XX) !'ATN(1/X)
!'ATN(1/X)=X/(X^2+1)*(1+(2)/(3)*1/(X^2+1)+(2*4)/(3*5)*(1/(X^2+1)^2+(2*4*6)/(3*5*7)*(1/(X^2+1))^3+...
OPTION ARITHMETIC NATIVE
DIM SS(-BIAS-1 TO KETA),Y(-BIAS-1 TO KETA)
LET Y(0)=1
LET Y(SIGN)=1
CALL SDIV(Y,XX*XX+1)
CALL SMUL(Y,XX)
CALL LCOPY(S,Y)
DO
   LET N=N+1
   CALL LCOPY(SS,S)
   CALL SMUL(Y,2*N)
   CALL SDIV(Y,2*N+1)
   CALL SDIV(Y,XX*XX+1)
   CALL LADD2(S,Y)
LOOP UNTIL EQUAL(SS,S)<>0
END SUB

以下省略
「n乗根の計算」から必要サブルーチンをコピペしてください
 

戻る