!=====================    円周率πの計算    ==============================!
!       π ＝ 32*Arctan(1/10) - 4*Arctan(1/239) - 16*Arctan(1/515)        !
! ========================================================================!
   MODULE multimal
      INTEGER :: mal = 10**5, prec
!
      INTERFACE OPERATOR(.plus.)
         MODULE PROCEDURE add
      END INTERFACE
      INTERFACE OPERATOR(.minus.)
         MODULE PROCEDURE subtract
      END INTERFACE
      INTERFACE OPERATOR(.times.)
         MODULE PROCEDURE multiply
      END INTERFACE
      INTERFACE OPERATOR(.div.)
         MODULE PROCEDURE divide
      END INTERFACE
!
   CONTAINS
!     --------------------------------------------------------------------
      FUNCTION add(x, y) RESULT(z)                   !  z = x + y
         INTEGER, INTENT(IN) :: x(0:prec), y(0:prec)
         INTEGER :: z(0:prec), i, zi, r
            r = 0
         DO i = prec, 0, -1
            zi = x(i) + y(i) + r
            r = zi/mal                               ! 繰り上がり
            z(i) = zi - r*mal
         END DO
      END FUNCTION add
!     --------------------------------------------------------------------
      FUNCTION subtract(x, y) RESULT(z)              !  z = x - y
         INTEGER, INTENT(IN) :: x(0:prec), y(0:prec)
         INTEGER :: z(0:prec), i, zi, r
            r = 1
         DO i = prec, 0, -1
            zi = x(i) + (mal - 1 - y(i)) + r
            r = zi/mal                               ! 繰り上がり
            z(i) = zi - r*mal
         END DO
      END FUNCTION subtract
!     --------------------------------------------------------------------
      FUNCTION multiply(x, s) RESULT(z)              !  z = x * s
         INTEGER, INTENT(IN) :: x(0:prec), s
         INTEGER :: z(0:prec), i, r, zi
            r = 0
         DO i = prec, 0, -1
            zi   = x(i)*s + r
            r    = zi/mal                            ! 繰り上がり
            z(i) = zi - r*mal
         END DO
      END FUNCTION multiply
!     --------------------------------------------------------------------
      FUNCTION divide(x, s) RESULT(z)                !  z = x / s
         INTEGER, INTENT(IN) :: x(0:prec), s
         INTEGER :: z(0:prec), i, r, zi
            r = 0 
         DO i = 0, prec
            zi   = x(i) + r*mal
            z(i) = zi/s
            r    = zi - s*z(i)                       ! 剰り
         END DO
      END FUNCTION divide
   END MODULE multimal
! ========================================================================
   PROGRAM pi_main                                   !!!    Main Program
      USE multimal
      INTEGER :: digit
      INTEGER, ALLOCATABLE :: pi(:)
!
      PRINT *," How many digits? (0<digit<100,000): "
      READ*, digit
      prec = CEILING(digit/5.0) + 1
      ALLOCATE( pi(0:prec) )
!
      pi =        (arctan( 10) .times. 32)       &
          .minus. (arctan(239) .times.  4)       &
          .minus. (arctan(515) .times. 16)
!
      PRINT "(2X, '***** πの計算 ( ', I6, ' 桁 ) *****')", digit
      PRINT "(' π = 3.', I5.5, 9I6.5/(7X, 10I6.5))", pi(1:prec-1)
   CONTAINS
!     --------------------------------------------------------------------
      FUNCTION arctan(k) RESULT(a)                   !!!   a = Arctan(1/k) 
         INTEGER, INTENT(IN)  :: k
         INTEGER :: a(0:prec), unity(0:prec), n 
         unity = (/ 1, (0, n = 1, prec) /)           ! unity = 1.00000 00000...
            a = 0
         DO n = INT( 0.5*digit/LOG10(REAL(k)) )+1, 1, -1 
            a = (((unity .div.(2*n+1)) .minus. a) .div. k) .div. k
         END DO
            a = (unity .minus. a) .div. k            ! for n = 0 
      END FUNCTION arctan
   END PROGRAM pi_main