従来ATAN方式のπ計算FORTRANプログラム

C=================================================================C
C  Atan公式によるπ計算プログラム(O(N**2)版、FORTRANプログラム)   C
C　　　　　　　　　　　　　　　　　　　　　　2002/9/26  後　保範  C
C=================================================================C
C  Pi Value Computation Program by Atan Normal 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-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,  2002/09/08                      C
C        (Hitachi Ltd. and Waseda University)                     C
C=================================================================C
      PARAMETER(ND=1000010, L10V=5)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /DEC10/P10,D10,N2,L10
      DIMENSION IV(ND), IXD(ND), INO(5), IDATA(2,4,5)
      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/
      CHARACTER*12 FILE1
      DATA FILE1/'AtanOrg?.txt'/ 	 
C   Constant Set
      L10 = L10V                ! L10 = 5
      P10 = 10.0D0**L10         ! P10 = 10**5
      D10 = 1.0D0/P10           ! D10 = 1/10**5
C   Parameter Read
      WRITE(6,*) ' Type In ID, (ID=1,2,3,4,5)'
      WRITE(6,*) ' 1;Machin, 2;Gauss, 3;Klingen., 4;Stoemer, 5;Takano'
      READ(5,*) ID              ! 計算公式選択
      IF(ID.LE.0 .or. ID.GE.6) ID = 1	 
      WRITE(6,*) ' Type In N (N=Number of Pi Decimal-Digits)'
      READ(5,*) NN              ! 計算桁数入力
      N = NN/L10
      IF(N.GE.ND-3) N = ND - 3
      NN = N*L10
      N2 = N + 3                ! 計算要素数を求める
      WRITE(6,300) NN
C   Output File Open
      FILE1(8:8) = CHAR(ID+ICHAR('0'))
      OPEN(UNIT=1,FILE=FILE1)
      CALL XCLOCK(T1,5)
C  Clear
      DO i=1,N2
        IV(i) = 0
      END DO
C  Pi Computation Start
      DO k=1,INO(ID)             ! Atanの個数だけ反復計算
        CALL XCLOCK(T2,5)
        M = IDATA(1,k,ID)        ! 当該AtanのM
        X = IDATA(2,k,ID)        ! 当該Atan(1/X)のX
        CALL ATANV(M,X,IV,IXD)   ! π :　IV = IV + M*Atan(1/X)  
        CALL XCLOCK(T3,5)
        IX = X
        WRITE(6,310) IX,T3-T2
      END DO
C  Output 
      CALL NORMAL(IV)             ! 各要素のを10進5桁にそろえる
      CALL XCLOCK(T4,5)
      WRITE(6,320) T4-T1	        ! πの計算時間
      WRITE(1,*) '  Pi = 3.'
      DO k=2,N2-2,16              ! 結果をファイルに出力
        KEND = MIN(k+15,N2-2)       ! 1行に10進80桁で
        KETA = (k-2)*L10 + 1
        WRITE(1,500) KETA,(IV(i),i=k,KEND) 
      END DO
C
      STOP
  300 FORMAT(1H ,' ** Pi Computation Start by Atan Original Method ** '/
     1     '  Number of Pi Decimal-Digits=',I8)
  310 FORMAT(1H ,'  Atan(1/',I6,') Computation,  Time(s)=',F7.2)
  320 FORMAT(1H ,'  Pi Computation,  Time(s)=',F7.2) 
  500 FORMAT(1H ,I8,8(1X,I5.5,I5.5))
      END
C=================================================================C
      SUBROUTINE ATANV(M,X,IV,IXD)
C=================================================================C
C   M*Atan(1/X) Computation by Normal Method                      C
C     M       I*4, In,  Value of M for M*Atan(1/X)                C
C     X       R*8, In,  Value of X for M*Atan(1/X)                C
C     IV(N2)  I*4, Out, Computaion Value of Atan(1/X)             C
C     IXD(N2) I*4, Wk,  Work for Computaion 1/X**K                C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,  2002/09/08                      C
C        (Hitachi Ltd. and Waseda University)                     C
C=================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /DEC10/P10,D10,N2,L10
      DIMENSION IV(N2), IXD(N2)
C  Initial Constant
      X2   = X**2
      DX2  = 1.0D0/X2
      KOSU = N2*L10/DLOG10(X2)     ! 多数桁に必要なAtanの項数
C  Initial Value Set
      DO i=2,N2
        IXD(i) = 0
      END DO
      IXD(1) = M
C  Initial Computation, IV=IXD=M/X
      W2 = 0.0
      DO i=1,N2
        W  = IXD(i) + W2*P10
        W1 = DINT(W/X)
        W2 = W - W1*X
        IXD(i) = W1
        IV(i)  = IV(i) + IXD(i)     ! IV = 前のAtanの計算値 + M/X
      END DO	 
C  Main Computation
      LST = 1
      ISG = -1
      DO k=1,KOSU                   ! 計算項数だけ反復計算 
        W2  = 0.0
        WK2 = 0.0
        VK  = 2*k + 1
        DK  = 1.0D0/VK
        DO i=LST,N2                  ! 多数桁は1要素(I*4)に
          W      = IXD(i) + W2*P10     ! 10進5桁ごとに分けて記憶
          W1     = DINT(W*DX2)       ! DINT: 切り捨てて整数値にする
          W2     = W - W1*X2  
          IXD(i) = W1                ! W1 = 多数桁のM/X**(2k+1)の値
          WK     = W1 + WK2*P10
          WK1    = DINT(WK*DK)
          WK2    = WK - WK1*VK       ! 多数桁のW1/(2k+1)
          IV(i)  = IV(i) + ISG*WK1
        END DO
        IF(IXD(LST).EQ.0) THEN
          LST = LST + 1               !  無駄な計算を省く処理
          IF(IXD(LST).EQ.0) LST = LST + 1
        END IF
        ISG = -ISG
      END DO
C
      RETURN
      END
C=================================================================C
      SUBROUTINE NORMAL(IV)         ! 多数桁の各要素を10進5桁に収める
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /DEC10/P10,D10,N2,L10
      DIMENSION IV(N2)
C 
      DO i=N2,2,-1
        W       = IV(i)            ! 10進5桁を超えたものは、一つ前の
        W1      = DINT(W*D10)        ! 要素に桁上げ処理をする
        IV(i-1) = IV(I-1) + W1
        IV(i)   = W - W1*P10
      END DO	 
C
      DO i=N2,2,-1                 ! 再度正確な桁揃え
        IF(IV(i).GE.P10) THEN        ! 10**5以上の要素の処理
          IV(i)   = IV(i) - P10
          IV(i-1) = IV(i-1) + 1
        END IF
        IF(IV(i).LT.0) THEN
          IV(i)   = IV(i) + P10       ! 負要素の処理
          IV(i-1) = IV(i-1) - 1
        END IF
      END DO
C
      RETURN
      END