!--
!  Copyright(C) 1996 Takuya OOURA <ooura@mmm.t.u-tokyo.ac.jp>
!
!  Redistribution and use in source and binary forms, with or without
!  modification, are permitted provided that the following conditions
!  are met:
!  1. Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!  2. Redistributions in binary form must reproduce the above copyright
!     notice, this list of conditions and the following disclaimer in the
!     documentation and/or other materials provided with the distribution.
!!
!  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!  ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
!  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
!  A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE HOLDERS OR
!  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
!  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
!  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
!  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
!  LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
!  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
!  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!++
!
!= gamma
!
! Version:: $Id$
!
!== 概要
!
! このファイルは gamma 関数に関する計算用の関数を提供する
!
!== 関数
!
! factrl  :: 与えられた整数に対する階乗(ガンマ関数)の値を返す
! gammaln :: 与えられた正の実数に対する階乗の対数(ガンマ関数の対数)の値を返す
!
function gammaln(xx)
  !
  ! Returns the vaule ln[Γ(x)] for x >= 0
  !
  real(8) ::  gammaln, xx
  real(8) ::  dlgamma
  external dlgamma

  gammaln = dlgamma(xx)

  return
end function gammaln

function factrl(n)
  !
  ! Returns the value n! as a floating point number
  !
  integer, intent(IN) :: n
  real(8)             :: factrl
  real(8)             :: gammaln
  external gammaln

  factrl=exp(gammaln(n+1.0d0))

  return

end function factrl