      PROGRAM integ
*     x = a から x = b までの区間の f(x) の積分を計算する
*     長方形近似を使う
      IMPLICIT NONE
      REAL*8 a, b, delta, dx, sum
      INTEGER k, n, kbhit, KBGETCH
      LOGICAL Ltmp1_
      CALL initial(a,b,sum,dx,delta,n)
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL sumf(a,b,sum,dx,delta)
         CALL output(sum,dx,n)
         CALL double(dx,delta,n)
         Ltmp1_ = kbhit() .EQ. 0
      END DO
      k = KBGETCH(0)
      END
C
      SUBROUTINE initial(a,b,sum,dx,delta,n)
      IMPLICIT NONE
      REAL*8 a, b, sum, dx, delta, f
      INTEGER n
      WRITE(*,'(A,$)') 'lower limit a = '
      READ(*,*) a
      WRITE(*,'(A,$)') 'upper limit b = '
      READ(*,*) b
      n = 2                                 ! 区間数の初期値
      dx = (b - a)/n
      delta = dx                            ! n = 2 の場合だけ delta = dx
      sum = f(a)                            ! 台形近似を使うときには (f(a) + f(b))/2 にする
      WRITE(*,*)
      WRITE(*,*) '       n rectangular approximation'
      WRITE(*,*)
      END
C
      SUBROUTINE sumf(a,b,sum,dx,delta)
      IMPLICIT NONE
      REAL*8 x, a, b, sum, dx, delta, f
      x = a + dx                            ! dx は積分の区間幅
      DO WHILE(x .LT. b)
         sum = sum + f(x)
         x = x + delta
      END DO
      END
C
      SUBROUTINE double(dx,delta,n)
*     区間数を2倍にする
      IMPLICIT NONE
      REAL*8 dx, delta
      INTEGER n
      delta = dx
      n = 2*n
      dx = 0.5D0*dx                         ! n > 2 では dx は delta に等しくない
      END
C
      SUBROUTINE output(sum,dx,n)
      IMPLICIT NONE
      REAL*8 sum, dx
      INTEGER n
      WRITE(*, '(I9,4X,F12.7)') n, sum*dx
      END
C
      FUNCTION f(x)
      IMPLICIT NONE
      REAL*8 f, x
      f = cos(x)
      END
