Ubasicから十進BASICへのお願い

 投稿者：GAI  投稿日：2009年11月16日(月)15時45分28秒
  220   point -10:word -40:M=1850:dim C(M):C(0)=1:S=1/2+14#i
250   for K=1 to M:for J=K to 1 step -1:C(J)=(C(J)+C(J-1))/2:next:next
300   ' find zeta zero
330   repeat
350     Z=fnZeta(S):H=Z/S:W=fnZeta(S+H):S+=H/(1-W/Z)
360     print using(8,20),S
380   until abs(Z)<1/10^18
390   end
800   ' zeta function
810   fnZeta(X)
820   local J,U
830   for J=1 to M:U+=(-1)^(J-1)*C(J)/J^X:next:U/=1-2^(1-X)
880   return(U)



のUbasicプログラムを十進BASICへ書き直して頂けませんか？ 		 

Re: Ubasicから十進BASICへのお願い

 投稿者：SECOND  投稿日：2009年11月16日(月)17時45分59秒
  > No.725[元記事へ]

GAIさんへ

!直訳です。精度は、かなり落ちます。uBASIC の有効桁(約2600桁)

OPTION ARITHMETIC COMPLEX
OPTION BASE 0
!                                 !point -10 !変数の小数部 10*4.8桁  設定不可。
!                                 !word -40  !変数の長さ   40*4.8桁  設定不可。
LET M=1850                                   !※短くても内部計算は常に2600桁。
DIM C(M)
LET C(0)=1
LET S=COMPLEX(1/2,14)             !S=1/2+14#i
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
!' find zeta zero
DO                                !repeat
   LET Z=Zeta(S)                  !Z=fnZeta(S)
   LET H=Z/S
   LET W=Zeta(S+H)                !W=fnZeta(S+H)
   LET S=S+H/(1-W/Z)              !S+=H/(1-W/Z)
   PRINT S                        !print using(8,20),S
   ! PRINT USING "########.#################### ########.####################":re(S),im(S);
   ! PRINT "#i"
LOOP UNTIL ABS(Z)< 1e-12 ! 1e-14  !until abs(Z)< 1/10^18  !十進BASICで ^18は無理。
STOP                              !END

!' zeta function
FUNCTION Zeta(X)                  !fnZeta(X) ※'fn'は、ユーザー定義関数 接頭語
   local J,U
   FOR J=1 TO M
      LET U=U+(-1)^(J-1)*C(J)/J^X !U+=(-1)^(J-1)*C(J)/J^X
   NEXT J
   LET U=U/(1-2^(1-X))            !U/=1-2^(1-X)
   LET Zeta=U                     !return(U)
END FUNCTION

END

! uBASIC での実行結果。
!run
!       0.50758974993024427495     +14.13558444989598483799#i
!       0.49997720386661645733     +14.13471625740151952033#i
!       0.49999999984206884658     +14.13472514154016738957#i
!       0.50000000000000000001     +14.13472514173469379043#i
!       0.50000000000000000000     +14.13472514173469379046#i
!OK
 		 

Re: Ubasicから十進BASICへのお願い

 投稿者：荒田浩二  投稿日：2009年11月17日(火)10時33分24秒
  > No.725[元記事へ]

1000桁モードで実行できるようにしました。
複素数は実部と虚部を分けて計算しています。超越関数は50桁の精度で計算しています。

結果
        .50844135172596863716       14.13022666245748735689
        .49997302379385895921       14.13474225141619583149
        .49999999977346757559       14.13472514196906007421
        .50000000000000000000       14.13472514173469379049
        .50000000000000000000       14.13472514173469379046
745 秒

OPTION ARITHMETIC DECIMAL_HIGH
220 ! point -10:word -40
    DECLARE EXTERNAL FUNCTION LOG,EXP
    DECLARE FUNCTION SIN,COS
    PUBLIC NUMERIC prec
    print time$
    LET t0=time
    LET prec=50 ! 精度
    LET M=1850
    DIM C(0 TO M),L(M),inv_fac(0 TO prec+1),S(2),Z(2),H(2),W(2),SH(2)
    LET C(0)=1
    LET S(1)=1/2
    LET S(2)=14
250 for K=1 to M
       LET L(K)=LOG(K)
       for J=K to 1 step -1
          LET C(J)=(C(J)+C(J-1))/2 ! C(K)=2^(-K)?
       next J
    next K
    LET inv_fac(0)=1 ! 階乗の逆数
    FOR i=1 TO prec+1
       LET inv_fac(i)=inv_fac(i-1)/i
    NEXT i
300 !' find zeta zero
330 DO ! repeat
350    CALL fnZeta(S,Z)
       LET H(1)=(Z(1)*S(1)+Z(2)*S(2))/(S(1)^2+S(2)^2) ! H=Z/S
       LET H(2)=(Z(2)*S(1)-Z(1)*S(2))/(S(1)^2+S(2)^2)
       MAT SH=S+H
       CALL fnZeta(SH,W) ! W=fnZeta(S+H)
       LET SH(1)=1-(W(1)*Z(1)+W(2)*Z(2))/(Z(1)^2+Z(2)^2)
       LET SH(2)=-(W(2)*Z(1)-W(1)*Z(2))/(Z(1)^2+Z(2)^2)
       LET S(1)=S(1)+(H(1)*SH(1)+H(2)*SH(2))/(SH(1)^2+SH(2)^2) ! S+=H/(1-W/Z)
       LET S(2)=S(2)+(H(2)*SH(1)-H(1)*SH(2))/(SH(1)^2+SH(2)^2)
360    print using "-"&REPEAT$("#",7)&"."&REPEAT$("#",20)&" -"&REPEAT$("#",7)&"."&REPEAT$("#",20):S(1),S(2)
380 LOOP until SQR(Z(1)^2+Z(2)^2)<1/10^18
    print int(time-t0);"秒"
390 !　end
800 !' zeta function
810 SUB fnZeta(X(),U())
820    local J
       MAT U=ZER
830    for J=1 to M
          LET p1=EXP(X(1)*L(J))*COS(X(2)*L(J))
          LET p2=EXP(X(1)*L(J))*SIN(X(2)*L(J))
          LET pp=1/(p1^2+p2^2)
          LET U(1)=U(1)+(-1)^(J-1)*C(J)*p1*pp ! U+=(-1)^(J-1)*C(J)/J^X
          LET U(2)=U(2)-(-1)^(J-1)*C(J)*p2*pp
       next J
       LET p1=1-EXP((1-X(1))*L(2))*COS(-X(2)*L(2))
       LET p2=-EXP((1-X(1))*L(2))*COS(-X(2)*L(2))
       LET pp=1/(p1^2+p2^2)
       LET u1=(U(1)*p1+U(2)*p2)*pp
       LET U(2)=(U(2)*p1-U(1)*p2)*pp ! U/=1-2^(1-X)
       LET U(1)=u1
880 END SUB ! return(U)
    !
    FUNCTION SIN(x)
       LET x=MOD(x,2*PI)
       LET sum=0
       FOR i=0 TO prec/2
          LET sum=sum+(-1)^i*x^(2*i+1)*inv_fac(2*i+1)
       NEXT i
       LET SIN=sum
    END FUNCTION
    FUNCTION COS(x)
       LET x=MOD(x,2*PI)
       LET sum=0
       FOR i=0 TO prec/2
          LET sum=sum+(-1)^i*x^(2*i)*inv_fac(2*i)
       NEXT i
       LET COS=sum
    END FUNCTION
END
!
! 1000桁モードで利用する対数関数(十進BASIC添付"\BASICw32\Library\log.LIB"参照)
EXTERNAL FUNCTION LOG(x)
    OPTION ARITHMETIC DECIMAL_HIGH
    IF x<=0 THEN
       CAUSE EXCEPTION 3004
    ELSEIF x<1 THEN
       LET log=-log(1/x)
    ELSEIF x>3 THEN
       LET log=2*log(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<=eps(0)*10^(1000-prec) ! 精度prec
       LET log=2*t
    END IF
END FUNCTION
!
! 1000桁モードで利用する指数関数(十進BASIC添付"\BASICw32\Library\exp.LIB"参照)
EXTERNAL FUNCTION EXP(x)
    OPTION ARITHMETIC DECIMAL_HIGH
    FUNCTION s(y,n)
       LET t=y*x/n
       IF ABS(t)<=EPS(0)*10^(1000-prec) THEN ! 精度prec
          LET s=y+t
       ELSE
          LET s=y+s(t,n+1)
       END IF
    END FUNCTION
    LET EXP=s(1,1)
END FUNCTION 		 

Re: Ubasicから十進BASICへのお願い

 投稿者：山中和義  投稿日：2009年11月17日(火)11時36分42秒
  > No.725[元記事へ]

みなさん、リーマン予想にはまっていますね？！
!リーマンのゼータ関数
!ζ(s)=1/(1-2^(1-s))*∑[m=0,∞]{2^(-(m+1))*∑[j=0,m]{(-1)^j*comb(m,j)*(j+1)^(-s)}}

LET t0=TIME


LET M=50 !1850
DIM C(0 TO M) !∑2^(-(M+1))*∑comb(M,J)
LET C(0)=1
FOR K=1 TO M !パスカルの三角形のM段より、二項係数comb(m,j)を求める
   FOR J=K TO 1 STEP -1
      LET C(J)=(C(J)+C(J-1))/2 !※2^(-(m+1))も加味する
   NEXT J
NEXT K
FUNCTION fnZeta(S) !リーマンのゼータ関数
   local J,U
   LET U=0
   FOR J=0 TO M
      LET U=U+(-1)^J*C(J)/(J+1)^S
   NEXT J
   LET U=U/(1-2^(1-S))
   LET fnZeta=U
END FUNCTION


PRINT C(M) !debug

PRINT fnZeta(2) - PI*PI/6


PRINT "計算時間=";TIME-t0

END

私も、実数の範囲で計算させてみて、多桁の複素数計算の実装を考えようとしたのですが、
荒田浩二さんがすでに着手されたようで、、、

1000桁モードでは二項係数を求めるのに、苦戦しているようで時間がかかっています。

10進モードや2進モードや複素数モードの１７桁程度では、M=50で十分でしょう。
５０桁精度の場合は、M=100で十分でしょう。 		 

戻る