!------------------------------------------------------------------------
! Copyright (c) 2011-2013 SPMODEL Development Group. All rights reserved.
!------------------------------------------------------------------------
!
!表題  dlgamma テストプログラム
!
!履歴  --
!      2011/02/18  佐々木洋平 dc_test を使用するように修正
!      2013/05/22  佐々木洋平 精度落ちを修正
!
program gamma_test
  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  implicit none
  real(8)  :: factrl
  real(8)  :: gammaln
  integer  :: i
  character(len=2) :: mess_int

  integer, parameter :: iend=10
  real(8), parameter :: pi = 3.141592653589793238d0
  !! 判定誤差設定
  integer, parameter :: check_digits = 12
  integer, parameter :: ignore = -13
  external factrl, gammaln

 call MessageNotify('M','gamma_test','test of gamma function')
 do i=0,iend
   write(mess_int, '(I2)') i
   call AssertEqual(                            &
     message = 'Gamma(  '//mess_int//')',       &
     answer  =  dble(ansi(i)),                  &
     check   =  factrl(i),                      &
     significant_digits = check_digits, ignore_digits = ignore   )
 end do

 do i=0,iend
   write(mess_int, '(I2)') i
   call AssertEqual(                      &
     message = 'Gamma('//mess_int//'.5)', &
     answer  =  ansr(i),                  &
     check   =  exp(gammaln(0.5D0+i)),    &
     significant_digits = check_digits, ignore_digits = ignore   )
 end do

 call MessageNotify('M','gamma_test','test of gamma function succeeded!')

contains
 recursive function ansi(i) result(res)
   integer(8) :: res
   integer, intent(in) :: i
   if ( i > 0 ) then
     res = i * ansi(i-1)
   else
     res = 1
   end if
 end function ansi

 recursive function ansr(i) result(res)
   real(8) :: res
   integer, intent(in) :: i
   if ( i > 0 ) then
     res = (i - 0.5d0) * ansr(i-1)
   else
     res = dsqrt(pi)
   end if
 end function ansr

end program gamma_test