!    繧ｪ繧､繝ｩ繝ｼ繝ｻ繝槭け繝ｭ繝ｼ繝ｪ繝ｳ縺ｮ蜈ｬ蠑上ｒ逕ｨ縺�◆霑台ｼｼ
!    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

! 蛻�牡謨ｰ縺茎etmax縺ｫ縺ｪ繧九∪縺ｧ繝ｫ繝ｼ繝励ｒ蝗槭☆
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