!    オイラー・マクローリンの公式を用いた近似
!    f(x)=e^(x) の 0<=x<=1 における積分
!    分割数が1, 2, 4, 8, 16, 32の場合
!    2013/08/01 村橋究理基 新規作成

program EM
implicit none

integer :: k,i     ! ループ変数
integer :: kp      ! 階乗計算用変数
integer :: set     ! 台形の分割数
integer :: setmax  ! 分割数の上限
real(8) :: max     ! 積分区間の上端
real(8) :: min     ! 積分区間の下端
real(8) :: fmax    ! 被積分関数の引数maxの値
real(8) :: fmin    ! 被積分関数の引数minの値
real(8) :: bn      ! ベルヌーイ数用変数
real(8) :: div     ! 分割した一区間の幅
real(8) :: sum     ! 部分和を保存する変数
real(8) :: f       ! 関数の値を計算する変数
real(8) :: res     ! 近似計算の結果を入れる変数
real(8) :: error   ! 近似解と解析解の差から得られる誤差


! 積分区間の上端と下端を設定
min = 99.0d0
max = 100.0d0

! 積分区間の上端と下端における関数fの値を計算
fmin = exp(min)
fmax = exp(max)

write(*,*) '分割数, 数値解, 誤差'

set = 1
setmax = 600

! 分割数がsetmaxになるまでループを回す
do while(set < setmax)
   do k = 1, 15
      
      ! 和を0で初期化
      sum = 0.0d0
      
      ! 分割幅
      div = (max - min)/dble(set)
      
      ! 階乗[(2k)!]の計算
      kp = 1
      do i = 1, 2*k
         kp = kp * i
      enddo

      call bnum(2*k, bn)

      sum = sum + bn/dble(kp) * (div**(2*k)) * (fmax - fmin)

   enddo

   error = sum
   ! 数値解の抜き出し
   res = error + (exp(1.0d0) - 1.0d0)

   ! 結果の出力
   write(*,*) set, res, error
   !write(*,*) set, '&', k, '&', error, '\\ \hline'
   !write(*,*) k, '&', error, '\\ \hline'
   !write(*,*) error, '\\ \hline'
!enddo
   set = set * 2
enddo

endprogram

! ベルヌーイ数を求めるサブルーチン
! 定義式を書きなおして得られる漸化式から求める
subroutine bnum(num, res)
implicit none

integer,intent(in) :: num           ! 求めたいベルヌーイ数添字
integer :: i           ! ループ用ベルヌーイ数添字
integer :: j           ! ループ変数(部分和計算用)
integer :: kp, lp, mp  ! 階乗計算用変数
real(8) :: Cij         ! Combination計算用変数
real(8) :: sum         ! 部分和を保存する変数
real(8),allocatable :: be(:)  !ベルヌーイ数として要素が100ある配列
real(8),intent(out) :: res  ! 返すベルヌーイ数

allocate(be(0:num))

be(0) = 1.0d0
i = 0

do i = 1, num
   Cij = 1.0d0

   sum = Cij * be(0)

   do j = 1, i-1
      ! (i+1)C(j)の計算
      Cij = Cij * dble(i+2-j)/dble(j)

      sum = sum + Cij*be(j)
   enddo

   be(i) = -1.0d0/dble(i+1)*sum

enddo

res = be(num)

deallocate(be)

endsubroutine bnum
