DRM法によるATAN公式π計算FORTRANプログラム

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

1. 概要

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

2. プログラム

C=================================================================C
C  DRM公式によるπ計算プログラム ( ATAN 版 , FORTRAN )
C             2003/01/18 Yasunori Ushiro ( 後 保範 )
C=================================================================C
C  Pi Value Computation Program by Atan with DRM Method           C
C    Pi = M1*Atan(1/X1) + M2*Atan(1/X2) + M3*Atan(1/X3)           C
C    Atan(1/X) = 1/X - 1/(3*X**3) + 1/(5*X**5) - 1/(7*X**7) +..   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), X(0:ND), Y(0:ND), PI(0:ND)
      DIMENSION E(0:ND), W(0:4*ND+1), WK(6*ND)
      DIMENSION INO(5), IDATA(2,4,5)
C  Pi ATAN Constant
      DATA INO/2,3,3,4,4/
      DATA IDATA/16, 5, -4,239,  0,0,     0,0,
     1           48,18, 32, 57, -20,239,  0,0,
     2           32,10, -4,239, -16,515,  0,0,
     3          176,57, 28,239, -48,682, 96,12943,
     4           48,49, 128,57, -20,239, 48,110443/
      COMMON /VBASE/BASE,NORID
C  Set Parameter (N;FFT Size, NC;Element Number; NN;Decimal Pi)
      CALL PIIN(2,ND,ID,N,NC,NN,L10)   ! 2;Atan with DRM
C  Set Computation Parameter
      CALL XCLOCK(T1,5)             ! Start Time
      ERR = 0.0
      CALL SETV(0,PI,N+1)           !  PI=float(0)
      DO k=1,INO(ID)
        CALL XCLOCK(T2,5)
        IX  = IDATA(2,k,ID)         !  Get IX ; ATAN(1/IX)
        IB  = IDATA(1,k,ID)         !  Get IB ; IB*ATAN(1/IX)
        XX  = IX
        ALG = DLOG(10.0D0)/DLOG(XX)       ! Terms by a Digits 
        MLP = 2.0*L10*(NC+4)*ALG/N + 1    ! Number of DRM steps
C    Y=ATAN(1/IX)
        CALL DRMATN(A,B,X,E,Y,IX,N,MLP,ER1,W,WK)
        ERR = MAX(ERR,ER1)
C    Pi=Pi+IB*Y
        CALL FCMT(Y,IB,X,N)         !  X=IB*Y
        CALL FADD(PI,X,Y,N,1)       !  Y=Pi+X
        CALL COPY(Y,PI,N+1)         !  Pi=Pi+IB*Y
        CALL XCLOCK(T3,5)
        WRITE(6,310) IX,MLP,T3-T2
      END DO
C  Output Pi
      CALL PIOUT(PI,NC,L10,T1,ERR)
C
      STOP
  310 FORMAT(1H ,'  ATAN(1/',I6,'),  DRM Step=',I3,',  Time(s)=',F8.2)
      END
C==================================================================C
      SUBROUTINE DRMATN(A,B,X,E,Y,IX,NN,MLP,ERR,W,WK)
C==================================================================C
C   Compute Y=ATAN(1/IX) by DRM Method                             C
C      Set Y,E,VS and Loop for MLP                                 C
C        Compute Y=Y+A/(B*E)                                       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      X(0:NN)   R*8, Work, Data X                                 C
C      E(0:NN)   R*8, Work, Data E                                 C
C      Y(0:NN)   R*8, Out,  Data Y=ATAN(1/IX)                      C
C      IX        I*4, In,   Value IX for ATAN(1/IX)                C
C      NN        I*4, In,   Number of Data Gloup                   C
C      MLP       I*4, In,   Number of Loops                        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), X(0:NN), Y(0:NN)
      DIMENSION W(0:4*NN), WK(6*NN)
C  Set Y,E
      CALL SETV(0,Y,NN+1)          !  Y=0
      CALL FVSET(1,E,NN)           !  E=Float(1)
      CALL FCDIV(E,IX,E,NN)        !  E=Float(IX)
      VS = 1.0
C  Main Loop
      ERR = 0.0
      DO k=1,MLP
C   Compute Y=Y+A/(B*E),   Compute A,B,X and E=E*X
        CALL DRMAT2(A,B,X,E,Y,IX,VS,NN,ER1,W,WK)
        ERR = MAX(ERR,ER1)
      END DO
C
      RETURN
      END
C==================================================================C
      SUBROUTINE DRMAT2(A,B,X,E,Y,IX,VS,NN,ERR,W,WK)
C==================================================================C
C   Outer Routine of DRM Method for ATAN                           C
C      Set Initial Value to A,B,X                                  C
C      Compute A,B,X by DRM and E=E*X                              C
C      Y=Y+A/(B*E) by 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      X(0:NN)   R*8, Work, Data X                                 C
C      E(0:NN)   R*8, I/O,  Data E                                 C
C      Y(0:NN)   R*8, I/O,  Data Y                                 C
C      IX        R*8, In,   Value IX for ATAN(1/IX)                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), X(0:NN), Y(0:NN), E(0:NN)
      DIMENSION W(0:4*NN), WK(6*NN)
C  Set Initial Value to A,B,X
      N = 4
      CALL SETV(0,A(1),NN)           !  Clear A
      CALL SETV(0,B(1),NN)           !  Clear B
      AV = 1.0
      DO k=1,NN/N
        B(k*N) = VS
        CALL NORM(B(k*(N-1)+1),N)    !  B(k)=VS+2*(k-1)
        A(k*N) = AV                  !  A(k)=(-1)**(k-1)
        VS     = VS + 2.0
        AV     = -AV
      END DO
      CALL SETV(0,X(1),N)            !  Clear X
      VX   = IX
      X(N) = VX*VX
      CALL NORM(X(1),N)              !  X=IX**2
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 DRMAT1(A(1),B(1),X(1),N,M,ER1,W,WK)
        ERR = MAX(ERR,ER1)
        N = N*2
        M = M/2
      END DO
C  Change Floating A,B,X from Fixed
      CALL FLOAT(A,NN)
      CALL FLOAT(B,NN)
      CALL FLOAT(X,NN)
C  E=E*X
      CALL FMTCFM(X,E,W,NN,ER1,WK)    ! W=E*X
      CALL COPY(W,E,NN+1)             ! E=E*X
C  Y=Y+A/(B*E)
      CALL FMTCFM(B,E,W,NN,ER2,WK)    ! W=B*E
      CALL COPY(W,X,NN+1)             ! X=B*E
      CALL FDIV(A,X,B,NN,ER3,W,WK)    ! B=A/X
      CALL FADD(Y,B,W,NN,1)           ! W=Y+A/X
      CALL COPY(W,Y,NN+1)             ! Y=Y+A(B*E)
      ERR = MAX(ERR,ER1,ER2,ER3)
C
      RETURN
  100 FORMAT(5F10.1)
      END
C==================================================================C
      SUBROUTINE DRMAT1(A,B,X,N,M,ERR,W,WK)
C==================================================================C
C   Inner Routine of DRM Method for ATAN                           C
C      W        = X(1)*B(2,k)                                      C
C      A(1-2,k) = A(1,k)*W + A(2,k)*B(1,k)                         C
C      B(1-2,k) = B(1,k)*B(2,k)       for k=1,M                    C
C        X(1-2)   = X(1)*X(1)                                      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      X(N,2)    R*8, I/O,   Data X                                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), X(N,2)
      DIMENSION W(N,4), WK(6*N)
C  Main Loop for M
      ERR = 0.0
      DO k=1,M
C   W=X*B(2,k)
        CALL MTCFM(X,B(1,2,k),W,N,ER1,WK)
        IF(W(N,1).NE.0.0D0) THEN
          WRITE(6,*) '** Error Stop **'
          STOP
        END IF
C   A(1-2,k)=A(1,k)*W+A(2,k)*B(1,k)
        CALL MTCFM(A(1,1,k),W(1,2),W(1,3),N,ER2,WK)  ! W(3)=A(1)*W(1)
        CALL MTCFM(A(1,2,k),B(1,1,k),W,N,ER3,WK)    ! W(1)=A(2)*B(1)
        CALL ADD(W,W(1,3),A(1,1,k),2*N,1)           ! A=W(1)+W(3)
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=B(1)*B(2)
        CALL COPY(W,B(1,1,k),2*N)                   ! B=W
        ERR = MAX(ERR,ER1,ER2,ER3,ER4)
      END DO
C  X(1-2)=X*X
      CALL COPY(X,W,N)                              ! W=X
      CALL MTCFM(W,W,X,N,ER1,WK)                    ! X=W*W
      ERR = MAX(ERR,ER1)
C
      RETURN
      END