DRM法によるラマヌジャン型公式π計算FORTRANプログラム

    2003/01/10 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------

1. 概要

 DRM法でのラマヌジャン型π計算ルーチンである。桁数nに対して、演算量はn・log(n)2〜3
に比例する。πを求める公式の2.2 ラマヌジャン型公式を参照。
 本計算には多数桁計算部品として 乗除算加減算及び π用I/Oの ルーチンが必要である。
 時間測定用ルーチンと XCLOCKして使用している。 このルーチンの仕様はHitachi
SR8000 FORTRANに合わせてある。この部分はシステムに合わせて変更が必要。

2. プログラム

C=================================================================C
C  DRM公式によるπ計算プログラム ( Ramanujan 版 , FORTRAN )
C             2003/01/18 Yasunori Ushiro ( 後 保範 )
C=================================================================C
C  Pi Value Computation Program by Ramanujan with DRM Method      C
C      Pi = K0*K1*SQRT(K2)/(K3*(Y+K4))                            C
C        Pi = 0                                                   C
C        do i=1,ML                                                C
C          Pi = (C + A*Y)/B  with A,B,C by DRM Method             C
C        end do                                                   C
C   Using FMT for High-Precision Multiplication                   C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,  2003/01/18                      C
C        (Hitachi Ltd. and Waseda University)                     C
C=================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(ND=270000)
      DIMENSION A(0:ND), B(0:ND), C(0:ND), X(0:ND), PI(0:ND)
      DIMENSION W(0:4*ND+1), WK(6*ND)
      DIMENSION KD(0:5,3), XX(2,3)
C  Pi Ramanujan Constant
C              K0   K1    K2  K3     K4    K5
      DATA KD/128,10005,10005, 3, 13591409, 5,   ! ID=1
     1          1, 9801,    2, 4,     1103, 3,   ! ID=2
     2          1, 3528,    1, 1,     1123, 3/   ! ID=3
      DATA XX/3.0,53360.0, 4.0,99.0, 2.0,882.0/
      COMMON /VBASE/BASE,NORID
C  Set Input Parameter
      CALL PIIN(3,ND,ID,N,NC,NN,L10)
C  Set MLP
      ALG = DLOG(10.0D0)/(XX(1,ID)*DLOG(XX(2,ID)))
      MLP = 8.0*L10*(NC+4)*ALG/N + 1
      WRITE(6,320) MLP
      WRITE(1,320) MLP
      CALL XCLOCK(T1,5)                  ! Start Time
      ERR = 0.0
C  Main Loop
      VV = N/8
      DO k=1,MLP
        VS = (MLP-k)*VV + 1.0            ! Start Value
        CALL XCLOCK(T2,5)
C   Computing A,B,C by DRM Method 
        CALL DRMRMA(A,B,C,ID,VS,N,ER1,W,WK)
C   Pi = (C+A*Pi)/B or Pi = C/B (k=1)
        IF(k.EQ.1) THEN 
          CALL FVSET(0,W,N)
          W(0) = -N                      ! W=Floating(0)
        ELSE
          CALL FMTCFM(A,PI,W,N,ER2,WK)   ! W=A*Pi
        END IF
        CALL FADD(C,W,X,N,1)             ! X=C+W
        CALL FDIV(X,B,PI,N,ER3,W,WK)     ! Pi=X/B
        ERR = MAX(ERR,ER1,ER2,ER3)
        CALL XCLOCK(T3,5)
        WRITE(6,310) k,T3-T2
      END DO
C  Pi=K0*K1*SQRT(K2)/(K3*(K4+K5*Pi))
      K2 = KD(2,ID)
      IF(K2.EQ.1) THEN
        CALL FVSET(1,X,N)                ! X=Floating(1)
      ELSE
        CALL FCSQRT(K2,X,N,1,ER1,W,WK)   ! X=SQRT(K2)
      END IF
      CALL FCMT(X,KD(0,ID),A,N)          ! A=K0*X
      CALL FCMT(A,KD(1,ID),B,N)          ! B=K1*A
      CALL FVSET(KD(4,ID),A,N)           ! A=K4
      CALL FCMT(PI,KD(5,ID),C,N)         ! C=K5*Pi
      CALL FADD(A,C,X,N,1)               ! X=A+C
      CALL FCMT(X,KD(3,ID),C,N)          ! C=K3*X
      CALL FDIV(B,C,PI,N,ER2,W,WK)       ! Pi=B/C
      ERR = MAX(ERR,ER1,ER2)
C  Output       
      CALL PIOUT(PI,NC,L10,T1,ERR)
C
      STOP
  310 FORMAT(1H ,'  Ramanujan Step=',I3,'  Time(s)=',F8.2)
  320 FORMAT(1H ,'  DRM Steps=',I3)
      END
C==================================================================C
      SUBROUTINE DRMRMA(A,B,C,ID,VS,NN,ERR,W,WK)
C==================================================================C
C   Outer Routine of DRM Method for Ramanujan                      C
C      Set Initial Value to A,B,C                                  C
C      Compute A,B,C by DRM and Floating                           C
C------------------------------------------------------------------C
C   ARGUMENTS :                                                    C
C      A(0:NN)   R*8, Work, Data A                                 C
C      B(0:NN)   R*8, Work, Data B                                 C
C      C(0:NN)   R*8, Work, Data X                                 C
C      ID        I*4, In,   Type of Ramanujan (1,2,3)              C
C      VS        R*8, I/O,  Start Value of Computation Terms       C
C      NN        I*4, In,   Number of Data Gloup                   C
C      ERR       R*8, Out,  Maximum Rounding Error to Convert Int. C
C      W(0:4*NN) R*8, Work, Work for DRM                           C
C      WK(6*NN) R*8, Work, Work for Multiplication                 C
C------------------------------------------------------------------C
C    Written by Yasunori Ushiro,  2003/01/18                       C
C        (Hitachi Ltd. and Waseda University)                      C
C==================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(0:NN), B(0:NN), C(0:NN)
      DIMENSION W(0:4*NN), WK(6*NN)
      DIMENSION IKA(6,3), IKB(12,3), CK(2,3)
C  A-Data
      DATA IKA/2,-1,6,-5,6,-1, 2,-1,4,-3,4,-1, 2,-1,4,-3,4,-1/
C  B-Data
      DATA IKB/3,0, 3,0, 1,0, 0,106720, 0,106720, 0,106720,
     1         2,0, 1,0, 1,0, 0,39204,  0,198,    0,198,
     2         2,0, 1,0, 1,0, 0,3528,   0,3528,   0,1/
C  C-Data
      DATA CK/545140134.D0,13591409.D0, 26390.D0, 1103.D0,
     1        21460.D0, 1123.D0/
C  Set Initial Value to A,B,C
      N  = 8                ! Number of Elements for A,B,C
      kk = 1                ! Start Address of A,B,C
      SG = 1                ! Sign of C
      DO k=1,NN/N
        IF(ID.NE.2) SG = -SG         ! Change Sign of C 
C   Set A-Data
        VS1 = VS + 1.0
        CALL VRMSET(A(kk),VS1,IKA(1,ID),N,3)   ! Set A
C   Set B-Data
        CALL VRMSET(B(kk),VS,IKB(1,ID),N,6)    ! Set B
C   Set C-Data 
        CALL SETV(0,C(kk),N)
        C(kk+N-1) = VS
        CALL NORM(C(kk),N)                 !  Normalize C
        CALL CMLT(C(kk),CK(1,ID),N)        !  C=C*IKC(1)
        C(kk+N-1) = C(kk+N-1) + CK(2,ID)   !  C=C+IKC(2)
        CALL NORM(C(kk),N)	  	  	     !  Normalize C
        CALL CMLT(C(kk),SG,N)              !  C=C*IS
        CALL NORM(C(kk),N)	  	  	     !  Normalize C
C   Count Up
        VS = VS + 1.0
        kk = kk + N
      END DO
C  Compute A,B,X by DRM
      M = NN/(2*N) 
      CALL LOG2F(2*M,MEP)                  !  MEP=log2(M)
      ERR = 0.0
      DO k=1,MEP
        CALL DRMRM1(A(1),B(1),C(1),N,M,ER1,W,WK)   ! DRM Step
        ERR = MAX(ERR,ER1)
        N = N*2
        M = M/2
      END DO
C  Change A,B,C from Fixed to Floating-Type
      CALL FLOAT(A,NN)
      CALL FLOAT(B,NN)
      CALL FLOAT(C,NN)
C
      RETURN
      END
C==================================================================C
      SUBROUTINE VRMSET(A,VS,IB,N,NO)
C==================================================================C
C   Set Initial Value for A-Data and B-Data                        C
C   ARGUMENTS :                                                    C
C      A(N)      R*8, I/O,  Input and Output Data                  C
C      VS        R*8, I/O,  Start Value of Computation Terms       C
C      IB(2,LP)  I*4, In,   Constant Table                         C
C      N         I*4, In,   Number of Elements in A                C
C      NO        I*4, In,   Number of Constant Table               C
C------------------------------------------------------------------C
C    Written by Yasunori Ushiro,  2003/01/18                       C
C        (Hitachi Ltd. and Waseda University)                      C
C==================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(N), IB(2,NO)
C
      CALL SETV(0,A,N)                 ! Clear A
      A(N) = 1.0                       ! A=1 by Fixed Type 
      DO k=1,NO
        X = IB(1,k)*VS + IB(2,k)
        CALL CMLT(A,X,N)               ! A=A*X by Fixed Type
      END DO
C
      RETURN
      END
C==================================================================C
      SUBROUTINE DRMRM1(A,B,C,N,M,ERR,W,WK)
C==================================================================C
C   Inner Routine of DRM Method for Ramanujan                      C
C      A(1-2,k) = A(1,k)*A(2,k)                                    C
C      B(1-2,k) = B(1,k)*B(2,k)                                    C
C      C(1-2,k) = A(1,k)*C(2,k) + B(2,k)*C(1,k) , for k=1,M        C
C------------------------------------------------------------------C
C   ARGUMENTS :                                                    C
C      A(N,2,M)  R*8, I/O,  Data A                                 C
C      B(N,2,M)  R*8, I/O,  Data B                                 C
C      C(N,2,M)  R*8, I/O,  Data C                                 C
C      N         I*4, In,   Number of Computation Elements         C
C      M         I*4, In,   Number of Data Gloup                   C
C      ERR       R*8, Out,  Maximum Rounding Error to Convert Int. C
C      W(N,4)    R*8, Work, Work for DRM                           C
C      WK(6*N)   R*8, Work, Work for Multiplication                C
C------------------------------------------------------------------C
C    Written by Yasunori Ushiro,  2003/01/18                       C
C        (Hitachi Ltd. and Waseda University)                      C
C==================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(N,2,M), B(N,2,M), C(N,2,M)
      DIMENSION W(N,4), WK(6*N)
C  Main Loop for M
      ERR = 0.0
      DO k=1,M
C   C(1-2,k) = A(1,k)*C(2,k) + B(2,k)*C(1,k)
        CALL MTCFM(A(1,1,k),C(1,2,k),W(1,1),N,ER1,WK)  ! W1=A1*C2
        CALL MTCFM(B(1,2,k),C(1,1,k),W(1,3),N,ER2,WK)  ! W3=B2*C1
        CALL ADD(W(1,1),W(1,3),C(1,1,k),2*N,1)         ! C=W1+W3
C   A(1-2,k) = A(1,k)*A(2,k) 
        CALL MTCFM(A(1,1,k),A(1,2,k),W,N,ER3,WK)       ! W=A1*A2
        CALL COPY(W,A(1,1,k),2*N)                      ! A=W 
C   B(1-2,k) = B(1,k)*B(2,k)
        CALL MTCFM(B(1,1,k),B(1,2,k),W,N,ER4,WK)       ! W=B1*B2
        CALL COPY(W,B(1,1,k),2*N)                      ! B=W
        ERR = MAX(ERR,ER1,ER2,ER3,ER4)
      END DO
C
      RETURN
      END