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

  integer, parameter :: iend=10
  ! 判定誤差設定 <-- これ，小数点以下の比較になってる
  integer, parameter :: check_digits = 8
  integer, parameter :: ignore = -8

  real(DP) :: check_eps

  check_eps = 1.0D0**(- check_digits) * EPSILON(1.0D0)

 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  =  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)
    real(DP) :: res
    integer, intent(in) :: i
    if ( i > 0 ) then
      res = i * ansi(i-1)
    else
      res = 1.0D0
    end if
 end function ansi

 recursive function ansr(i) result(res)
   real(DP) :: 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
