#:include "common.fypp"
#:set R_KINDS_TYPES = [KT for KT in REAL_KINDS_TYPES if KT[0] in ["sp","dp"]]
#:set C_KINDS_TYPES = [KT for KT in CMPLX_KINDS_TYPES if KT[0] in ["sp","dp"]]
#:set CI_KINDS_TYPES = INT_KINDS_TYPES + C_KINDS_TYPES
module stdlib_specialfunctions_gamma
use iso_fortran_env, only : qp => real128
use stdlib_kinds, only : sp, dp, int8, int16, int32, int64
use stdlib_error, only : error_stop
implicit none
private
integer(int8), parameter :: max_fact_int8 = 6_int8
integer(int16), parameter :: max_fact_int16 = 8_int16
integer(int32), parameter :: max_fact_int32 = 13_int32
integer(int64), parameter :: max_fact_int64 = 21_int64
#:for k1, t1 in R_KINDS_TYPES
${t1}$, parameter :: tol_${k1}$ = epsilon(1.0_${k1}$)
#:endfor
real(qp), parameter :: tol_qp = epsilon(1.0_qp)
public :: gamma, log_gamma, log_factorial
public :: lower_incomplete_gamma, log_lower_incomplete_gamma
public :: upper_incomplete_gamma, log_upper_incomplete_gamma
public :: regularized_gamma_p, regularized_gamma_q
interface gamma
!! Gamma function for integer and complex numbers
!!
#:for k1, t1 in CI_KINDS_TYPES
module procedure gamma_${t1[0]}$${k1}$
#:endfor
end interface gamma
interface log_gamma
!! Logarithm of gamma function
!!
#:for k1, t1 in CI_KINDS_TYPES
module procedure l_gamma_${t1[0]}$${k1}$
#:endfor
end interface log_gamma
interface log_factorial
!! Logarithm of factorial n!, integer variable
!!
#:for k1, t1 in INT_KINDS_TYPES
module procedure l_factorial_${t1[0]}$${k1}$
#:endfor
end interface log_factorial
interface lower_incomplete_gamma
!! Lower incomplete gamma function
!!
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
module procedure ingamma_low_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
module procedure ingamma_low_${t1[0]}$${k1}$
#:endfor
end interface lower_incomplete_gamma
interface log_lower_incomplete_gamma
!! Logarithm of lower incomplete gamma function
!!
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
module procedure l_ingamma_low_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
module procedure l_ingamma_low_${t1[0]}$${k1}$
#:endfor
end interface log_lower_incomplete_gamma
interface upper_incomplete_gamma
!! Upper incomplete gamma function
!!
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
module procedure ingamma_up_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
module procedure ingamma_up_${t1[0]}$${k1}$
#:endfor
end interface upper_incomplete_gamma
interface log_upper_incomplete_gamma
!! Logarithm of upper incomplete gamma function
!!
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
module procedure l_ingamma_up_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
module procedure l_ingamma_up_${t1[0]}$${k1}$
#:endfor
end interface log_upper_incomplete_gamma
interface regularized_gamma_p
!! Regularized (normalized) lower incomplete gamma function, P
!!
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
module procedure regamma_p_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
module procedure regamma_p_${t1[0]}$${k1}$
#:endfor
end interface regularized_gamma_p
interface regularized_gamma_q
!! Regularized (normalized) upper incomplete gamma function, Q
!!
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
module procedure regamma_q_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
module procedure regamma_q_${t1[0]}$${k1}$
#:endfor
end interface regularized_gamma_q
interface gpx
! Incomplete gamma G function.
! Internal use only
!
#:for k1, t1 in R_KINDS_TYPES
module procedure gpx_${t1[0]}$${k1}$ !for real p and x
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
module procedure gpx_${t1[0]}$${k1}$${k2}$ !for integer p and real x
#:endfor
#:endfor
end interface gpx
interface l_gamma
! Logarithm of gamma with integer argument for designated output kind.
! Internal use only
!
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
module procedure l_gamma_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
end interface l_gamma
contains
#:for k1, t1 in INT_KINDS_TYPES
impure elemental function gamma_${t1[0]}$${k1}$(z) result(res)
${t1}$, intent(in) :: z
${t1}$ :: res, i
${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$
if(z <= zero) call error_stop("Error(gamma): Gamma function argument" &
//" must be positive integer.")
if(z > max_fact_${k1}$) call error_stop("Error(gamma): Gamma function" &
//" integer argument is greater than the upper limit from which an"&
//" integer overflow will be generated. Suggest switch to high " &
//" precision or convert to real data type")
res = one
do i = one, z - one
res = res * i
end do
end function gamma_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in C_KINDS_TYPES
#:if k1 == "sp"
#:set k2 = "dp"
#:elif k1 == "dp"
#:set k2 = "qp"
#:endif
#:set t2 = "real({})".format(k2)
impure elemental function gamma_${t1[0]}$${k1}$(z) result(res)
${t1}$, intent(in) :: z
${t1}$ :: res
integer :: i
real(${k1}$), parameter :: zero_k1 = 0.0_${k1}$
${t2}$, parameter :: half = 0.5_${k2}$, &
one = 1.0_${k2}$, pi = acos(- one), sqpi = sqrt(pi)
complex(${k2}$) :: y, x, sum
#:if k1 == "sp"
#! for single precision input, using double precision for calculation
integer, parameter :: n = 10
${t2}$, parameter :: r = 10.900511_${k2}$
${t2}$, parameter :: d(0 : n) = [2.48574089138753566e-5_${k2}$, &
1.05142378581721974_${k2}$, &
-3.45687097222016235_${k2}$, &
4.51227709466894824_${k2}$, &
-2.98285225323576656_${k2}$, &
1.05639711577126713_${k2}$, &
-1.95428773191645870e-1_${k2}$, &
1.70970543404441224e-2_${k2}$, &
-5.71926117404305781e-4_${k2}$, &
4.63399473359905637e-6_${k2}$, &
-2.71994908488607704e-9_${k2}$]
! parameters from above referenced source.
#:elif k1 == "dp"
#! for double precision input, using quadruple precision for calculation
integer, parameter :: n = 24
${t2}$, parameter :: r = 25.617904_${k2}$
${t2}$, parameter :: d(0 : n)= &
[1.0087261714899910504854136977047144166e-11_${k2}$, &
1.6339627701280724777912729825256860624_${k2}$, &
-1.4205787702221583745972794018472259342e+1_${k2}$, &
5.6689501646428786119793943350900908698e+1_${k2}$, &
-1.3766376824252176069406853670529834070e+2_${k2}$, &
2.2739972766608392140035874845640820558e+2_${k2}$, &
-2.7058382145757164380300118233258834430e+2_${k2}$, &
2.39614374587263042692333711131832094166e+2_${k2}$, &
-1.6090450559507517723393498276315290189e+2_${k2}$, &
8.27378183187161305711485619113605553100e+1_${k2}$, &
-3.2678977082742592701862249152153110206e+1_${k2}$, &
9.89018079175824824537131521501652931756_${k2}$, &
-2.2762136356329318377213053650799013041_${k2}$, &
3.93265017303573867227590563182750070164e-1_${k2}$, &
-5.0051054352146209116457193223422284239e-2_${k2}$, &
4.57142601898244576789629257292603538238e-3_${k2}$, &
-2.8922592124650765614787233510990416584e-4_${k2}$, &
1.20833375377219592849746118012697473202e-5_${k2}$, &
-3.1220812187551248389268359432609135033e-7_${k2}$, &
4.55117045361638520378367871355819524460e-9_${k2}$, &
-3.2757632817493581828033170342853173968e-11_${k2}$, &
9.49784279240135747819870224486376897253e-14_${k2}$, &
-7.9480594917454410117072562195702526836e-17_${k2}$, &
1.04692819439870077791406760109955648941e-20_${k2}$, &
-5.8990280044857540075384586350723191533e-26_${k2}$]
! parameters from above referenced source.
#:endif
if(abs(z % im) < tol_${k1}$) then
res = cmplx(gamma(z % re), kind = ${k1}$)
return
end if
if(z % re < zero_k1) then
x = cmplx(abs(z % re), - z % im, kind = ${k1}$)
y = x - one
else
y = z - one
end if
sum = cmplx(d(0), kind = ${k2}$)
do i = 1, n
sum = sum + d(i) / (y + i)
end do
y = exp((y + half) * log(y + half + r) - y) * sum
y = y * 2 / sqpi !Re(z) > 0 return
if(z % re < zero_k1 ) then
y = - pi / (sin(pi * x) * x * y) !Re(z) < 0 return
end if
res = y
end function gamma_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
impure elemental function l_gamma_${t1[0]}$${k1}$(z) result(res)
!
! Logarithm of gamma function for integer input
!
${t1}$, intent(in) :: z
real :: res
${t1}$ :: i
${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$
if(z <= zero) call error_stop("Error(log_gamma): Gamma function" &
//" argument must be positive integer.")
select case(z)
case (one)
res = 0.0
case (two :)
res = 0.0
do i = one, z - one
res = res + log(real(i))
end do
end select
end function l_gamma_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
impure elemental function l_gamma_${t1[0]}$${k1}$${k2}$(z, x) result(res)
!
! Logarithm of gamma function for integer input with defined precision output
!
${t1}$, intent(in) :: z
${t2}$, intent(in) :: x
${t2}$ :: res
${t1}$ :: i
${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$
${t2}$, parameter :: zero_k2 = 0.0_${k2}$
if(z <= zero) call error_stop("Error(log_gamma): Gamma function" &
//" argument must be positive integer.")
select case(z)
case (one)
res = zero_k2
case (two :)
res = zero_k2
do i = one, z - one
res = res + log(real(i, ${k2}$))
end do
end select
end function l_gamma_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in C_KINDS_TYPES
#:if k1 == "sp"
#:set k2 = "dp"
#:elif k1 == "dp"
#:set k2 = "qp"
#:endif
#:set t2 = "real({})".format(k2)
impure elemental function l_gamma_${t1[0]}$${k1}$(z) result (res)
!
! log_gamma function for any complex number, excluding negative whole number
! "Computation of special functions", Shanjie Zhang & Jianmin Jin, 1996, p.48
! "Computing the principal branch of log-gamma", D.E.G. Hare,
! J. of Algorithms, 25(2), 1997 p. 221–236
!
! Fortran 90 program by Jim-215-Fisher
!
${t1}$, intent(in) :: z
${t1}$ :: res, z1, z2
real(${k1}$) :: d
integer :: m, i
complex(${k2}$) :: zr, zr2, sum, s
real(${k1}$), parameter :: z_limit = 10_${k1}$, zero_k1 = 0.0_${k1}$
integer, parameter :: n = 20
${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$, &
pi = acos(-one), ln2pi = log(2 * pi)
${t2}$, parameter :: a(n) = [ &
.8333333333333333333333333333333333333333E-1_${k2}$,&
-.2777777777777777777777777777777777777778E-2_${k2}$,&
.7936507936507936507936507936507936507937E-3_${k2}$,&
-.5952380952380952380952380952380952380952E-3_${k2}$,&
.8417508417508417508417508417508417508418E-3_${k2}$,&
-.1917526917526917526917526917526917526918E-2_${k2}$,&
.6410256410256410256410256410256410256410E-2_${k2}$,&
-.2955065359477124183006535947712418300654E-1_${k2}$,&
.1796443723688305731649384900158893966944E+0_${k2}$,&
-.1392432216905901116427432216905901116427E+1_${k2}$,&
.1340286404416839199447895100069013112491E+2_${k2}$,&
-.1568482846260020173063651324520889738281E+3_${k2}$,&
.2193103333333333333333333333333333333333E+4_${k2}$,&
-.3610877125372498935717326521924223073648E+5_${k2}$,&
.6914722688513130671083952507756734675533E+6_${k2}$,&
-.1523822153940741619228336495888678051866E+8_${k2}$,&
.3829007513914141414141414141414141414141E+9_${k2}$,&
-.1088226603578439108901514916552510537473E+11_${k2}$,&
.3473202837650022522522522522522522522523E+12_${k2}$,&
-.1236960214226927445425171034927132488108E+14_${k2}$]
! parameters from above reference
z2 = z
if(z % re < zero_k1) then
z2 = cmplx(abs(z % re), - z % im, kind = ${k1}$) + 1
end if
d = hypot(z2 % re, z2 % im)
z1 = z2
m = 0
if(d <= z_limit) then !for small |z|
m = ceiling(z_limit - d)
z1 = z2 + m
end if
zr = one / z1
zr2 = zr * zr
sum = (((a(20) * zr2 + a(19)) * zr2 + a(18)) * zr2 + a(17)) * zr2
sum = (((sum + a(16)) * zr2 + a(15)) * zr2 + a(14)) * zr2
sum = (((sum + a(13)) * zr2 + a(12)) * zr2 + a(11)) * zr2
sum = (((sum + a(10)) * zr2 + a(9)) * zr2 + a(8)) * zr2
sum = (((sum + a(7)) * zr2 + a(6)) * zr2 + a(5)) * zr2
sum = (((sum + a(4)) * zr2 + a(3)) * zr2 + a(2)) * zr2
sum = (sum + a(1)) * zr + ln2pi / 2 - z1 + (z1 - 0.5_${k2}$) * log(z1)
if(m /= 0) then
s = cmplx(zero, zero, kind = ${k2}$)
do i = 1, m
s = s + log(cmplx(z1, kind = ${k2}$) - i)
end do
sum = sum - s
end if
if(z % re < zero_k1) then
sum = log(pi) - log(sin(pi * z)) - sum
m = ceiling((2 * z % re - 3) / 4)
sum % im = sum % im + 2 * pi * m * sign(1.0_${k1}$, z % im)
end if
res = cmplx(sum, kind = ${k1}$)
end function l_gamma_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
impure elemental function l_factorial_${t1[0]}$${k1}$(n) result(res)
!
! Log(n!)
!
${t1}$, intent(in) :: n
real(dp) :: res
${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$
real(dp), parameter :: zero_dp = 0.0_dp
if(n < zero) call error_stop("Error(l_factorial): Logarithm of" &
//" factorial function argument must be non-negative")
select case(n)
case (zero)
res = zero_dp
case (one)
res = zero_dp
case (two : )
res = l_gamma(n + 1, 1.0_dp)
end select
end function l_factorial_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in R_KINDS_TYPES
#:if k1 == "sp"
#:set k2 = "dp"
#:elif k1 == "dp"
#:set k2 = "qp"
#:endif
#:set t2 = "real({})".format(k2)
impure elemental function gpx_${t1[0]}$${k1}$(p, x) result(res)
!
! Approximation of incomplete gamma G function with real argument p.
!
! Based on Rémy Abergel and Lionel Moisan "Algorithm 1006, Fast and
! Accurate Evaluation of a Generalized Incomplete Gamma Function", ACM
! Transactions on Mathematical Software, March 2020.
!
! Fortran 90 program by Jim-215-Fisher
!
${t1}$, intent(in) :: p, x
integer :: n, m
${t2}$ :: res, p_lim, a, b, g, c, d, y, ss
${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$
${t2}$, parameter :: dm = tiny(1.0_${k2}$) * 10 ** 6
${t1}$, parameter :: zero_k1 = 0.0_${k1}$
if(p <= zero_k1) call error_stop("Error(gpx): Incomplete gamma" &
//" function must have a positive parameter p")
if(x < -9.0_${k1}$) then
p_lim = 5.0_${k1}$ * (sqrt(abs(x)) - 1.0_${k1}$)
elseif(x >= -9.0_${k1}$ .and. x <= zero_k1) then
p_lim = zero_k1
else
p_lim = x
endif
if(x < zero_k1 .and. p < p_lim .and. abs(anint(p) - p) > tol_${k1}$) &
call error_stop("Error(gpx): Incomplete gamma function with " &
//"negative x must come with a whole number p not too small")
if(p >= p_lim) then !use modified Lentz method of continued fraction
!for eq. (15) in the above reference.
a = one
b = p
g = a / b
c = a / dm
d = one / b
n = 2
do
if(mod(n, 2) == 0) then
a = (one - p - n / 2) * x
else
a = (n / 2) * x
end if
b = p - one + n
d = d * a + b
if(d == zero) d = dm
c = b + a / c
if(c == zero) c = dm
d = one / d
y = c * d
g = g * y
n = n + 1
if(abs(y - one) < tol_${k2}$) exit
end do
else if(x >= zero_k1) then !use modified Lentz method of continued
!fraction for eq. (16) in the reference.
a = one
b = x + one - p
g = a / b
c = a / dm
d = one / b
n = 2
do
a = (n - 1) * (1 + p - n)
b = b + 2
d = d * a + b
if(d == zero) d = dm
c = b + a / c
if(c == zero) c = dm
d = one / d
y = c * d
g = g * y
n = n + 1
if(abs(y - one) < tol_${k2}$) exit
end do
else !Algorithm 2 in the reference
m = nint(ss)
a = - x
c = one / a
d = p - one
b = c * (a - d)
n = 1
do
c = d * (d - one) / (a * a)
d = d - 2
y = c * (a - d)
b = b + y
n = n + 1
if(n > int((p - 2) / 2) .or. y < b * tol_${k2}$) exit
end do
if(y >= b * tol_${k2}$ .and. mod(m , 2) /= 0) b = b + d * c / a
g = ((-1) ** m * exp(-a + log_gamma(p) - (p - 1) * log(a)) + b) / a
end if
res = g
end function gpx_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
impure elemental function gpx_${t1[0]}$${k1}$${k2}$(p, x) result(res)
!
! Approximation of incomplete gamma G function with integer argument p.
!
! Based on Rémy Abergel and Lionel Moisan "Algorithm 1006, Fast and
! Accurate Evaluation of a Generalized Incomplete Gamma Function", ACM
! Transactions on Mathematical Software, March 2020.
!
${t1}$, intent(in) :: p
${t2}$, intent(in) :: x
${t2}$ :: res, p_lim, a, b, g, c, d, y
integer :: n
${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$
${t2}$, parameter :: dm = tiny(1.0_${k2}$) * 10 ** 6
${t1}$, parameter :: zero_k1 = 0_${k1}$, two = 2_${k1}$
if(p <= zero_k1) call error_stop("Error(gpx): Incomplete gamma " &
//"function must have a positive parameter p")
if(x < -9.0_${k2}$) then
p_lim = 5.0_${k2}$ * (sqrt(abs(x)) - 1.0_${k2}$)
else if(x >= -9.0_${k2}$ .and. x <= zero) then
p_lim = zero
else
p_lim = x
end if
if(real(p, ${k2}$) >= p_lim) then
a = one
b = p
g = a / b
c = a / dm
d = one / b
n = 2
do
if(mod(n, 2) == 0) then
a = (1 - p - n / 2) * x
else
a = (n / 2) * x
end if
b = p - 1 + n
d = d * a + b
if(d == zero) d = dm
c = b + a / c
if(c == zero) c = dm
d = one / d
y = c * d
g = g * y
n = n + 1
if(abs(y - one) < tol_${k2}$) exit
end do
else if(x >= zero) then
a = one
b = x + 1 - p
g = a / b
c = a / dm
d = one / b
n = 2
do
a = -(n - 1) * (n - 1 - p)
b = b + 2
d = d * a + b
if(d == zero) d = dm
c = b + a / c
if(c == zero) c = dm
d = one / d
y = c * d
g = g * y
n = n + 1
if(abs(y - one) < tol_${k2}$) exit
end do
else
a = -x
c = one / a
d = p - 1
b = c * (a - d)
n = 1
do
c = d * (d - one) / (a * a)
d = d - 2
y = c * ( a - d)
b = b + y
n = n + 1
if(int(n, ${k1}$) > (p - two) / two .or. y < b * tol_${k2}$) exit
end do
if(y >= b * tol_${k2}$ .and. mod(p, two) /= zero_k1) &
b = b + d * c / a
g = ((-1) ** p * exp(-a + l_gamma(p, one) - (p - 1) * log(a)) &
+ b ) / a
end if
res = g
end function gpx_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
impure elemental function ingamma_low_${t1[0]}$${k1}$(p, x) result(res)
!
! Approximation of lower incomplete gamma function with real p.
!
${t1}$, intent(in) :: p, x
${t1}$ :: res, s1, y
${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
if(x == zero) then
res = zero
else if(x > p) then
s1 = log_gamma(p)
y = one - exp(-x + p * log(x) - s1) * gpx(p, x)
res = exp(s1 + log(y))
else if(x <= p .and. x > zero) then
s1 = -x + p * log(x)
res = gpx(p, x) * exp(s1)
else
call error_stop("Error(Logarithm of upper incomplete gamma " &
//"function): negative x must be with integer p")
end if
end function ingamma_low_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
impure elemental function ingamma_low_${t1[0]}$${k1}$${k2}$(p, x) &
result(res)
!
! Approximation of lower incomplete gamma function with integer p.
!
${t1}$, intent(in) :: p
${t2}$, intent(in) :: x
${t2}$ :: res, s1, y
${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$
if(x == zero) then
res = zero
else if(x > real(p, ${k2}$)) then
s1 = l_gamma(p, one)
y = one - exp(-x + p * log(x) - s1) * gpx(p, x)
res = exp(s1 + log(y))
else if(x <= real(p, ${k2}$) .and. x > zero) then
s1 = -x + p * log(x)
res = gpx(p, x) * exp(s1)
else
s1 = -x + p * log(abs(x))
res = gpx(p, x) * exp(s1)
res = (-1) ** p * res
end if
end function ingamma_low_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
impure elemental function l_ingamma_low_${t1[0]}$${k1}$(p, x) result(res)
${t1}$, intent(in) :: p, x
${t1}$ :: res, s1, y
${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
if(x == zero) then
res = zero
else if(x > p) then
s1 = log_gamma(p)
y = one - exp(-x + p * log(x) - s1) * gpx(p, x)
res = s1 + log(y)
else if(x <= p .and. x > zero) then
s1 = -x + p * log(abs(x))
res = log(abs(gpx(p, x))) + s1
else
call error_stop("Error(Logarithm of upper incomplete gamma " &
//"function): negative x must be with integer p")
end if
end function l_ingamma_low_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
impure elemental function l_ingamma_low_${t1[0]}$${k1}$${k2}$(p, x) &
result(res)
${t1}$, intent(in) :: p
${t2}$, intent(in) :: x
${t2}$ :: res, s1, y
${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$
if(x == zero) then
res = zero
else if(x > real(p, ${k2}$)) then
s1 = l_gamma(p, one)
y = one - exp(-x + p * log(x) - s1) * gpx(p, x)
res = s1 + log(y)
else if(x <= real(p, ${k2}$)) then
s1 = -x + p * log(abs(x))
res = log(abs(gpx(p, x))) + s1
end if
end function l_ingamma_low_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
impure elemental function ingamma_up_${t1[0]}$${k1}$(p, x) result(res)
!
! Approximation of upper incomplete gamma function with real p.
!
${t1}$, intent(in) :: p, x
${t1}$ :: res, s1, y
${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
if(x == zero) then
res = gamma(p)
else if(x > p) then
s1 = -x + p * log(x)
res = gpx(p, x) * exp(s1)
else if(x <= p .and. x > zero) then
y = log_gamma(p)
s1 = -x + p * log(x) - y
res = (one - gpx(p, x) * exp(s1)) * exp(y)
else
call error_stop("Error(Logarithm of upper incomplete gamma " &
//"function): negative x must be with integer p")
end if
end function ingamma_up_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
impure elemental function ingamma_up_${t1[0]}$${k1}$${k2}$(p, x) &
result(res)
!
! Approximation of upper incomplete gamma function with integer p.
!
${t1}$, intent(in) :: p
${t2}$, intent(in) :: x
${t2}$ :: res, s1, y
${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$
if(x == zero) then
res = gamma(real(p, ${k2}$))
else if(x > real(p, ${k2}$)) then
s1 = -x + p * log(x)
res = gpx(p, x) * exp(s1)
else if(x <= real(p, ${k2}$) .and. x > zero) then
y = l_gamma(p, one)
s1 = -x + p * log(x) - y
res = gpx(p, x) * exp(s1)
res = (one - res) * exp(y)
else
y = l_gamma(p, one)
s1 = -x + p * log(abs(x)) - y
res = gpx(p, x) * exp(s1)
res = (one - (-1) ** p * res) * exp(y)
end if
end function ingamma_up_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
impure elemental function l_ingamma_up_${t1[0]}$${k1}$(p, x) result(res)
${t1}$, intent(in) :: p, x
${t1}$ :: res, s1, y
${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
if(x == zero) then
res = log_gamma(p)
else if(x > p) then
s1 = -x + p * log(x)
res = log(gpx(p, x)) + s1
else if(x <= p .and. x > zero) then
y= log_gamma(p)
s1 = -x + p * log(x) - y
res = gpx(p, x) * exp(s1)
res = log(one - res) + y
else
call error_stop("Error(Logarithm of upper incomplete gamma " &
//"function): negative x must be with integer p")
end if
end function l_ingamma_up_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
impure elemental function l_ingamma_up_${t1[0]}$${k1}$${k2}$(p, x) &
result(res)
${t1}$, intent(in) :: p
${t2}$, intent(in) :: x
${t2}$ :: res, s1, y
${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$
if(x == zero) then
res = l_gamma(p, one)
else if(x > real(p, ${k2}$)) then
s1 = -x + p * log(x)
res = log(gpx(p, x)) + s1
else if(x <= real(p, ${k2}$) .and. x > zero) then
y = l_gamma(p, one)
s1 = -x + p * log(x) - y
res = gpx(p, x) * exp(s1)
res = log(one - res) + y
else
y = l_gamma(p, one)
s1 = -x + p * log(abs(x)) + log(gpx(p, x))
res = (-1) ** p * exp(s1)
res = log(abs(exp(y) - res))
end if
end function l_ingamma_up_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
impure elemental function regamma_p_${t1[0]}$${k1}$(p, x) result(res)
!
! Approximation of regularized incomplete gamma function P(p,x) for real p
!
${t1}$, intent(in) :: p, x
${t1}$ :: res, s1
${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_p" &
//" function is not defined at x < 0")
if(x == zero) then
res = zero
else if(x > p) then
s1 = -x + p * log(x) - log_gamma(p)
res = one - exp(s1 + log(gpx(p,x)))
else if(x <= p) then
s1 = -x + p * log(abs(x)) - log_gamma(p)
res = exp(log(gpx(p, x)) + s1)
end if
end function regamma_p_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
impure elemental function regamma_p_${t1[0]}$${k1}$${k2}$(p, x) result(res)
!
! Approximation of regularized incomplete gamma function P(p,x) for integer p
!
${t1}$, intent(in) :: p
${t2}$, intent(in) :: x
${t2}$ :: res, s1
${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$
if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_p" &
//" function is not defined at x < 0")
if(x == zero) then
res = zero
else if(x > real(p, ${k2}$)) then
s1 = -x + p * log(x) - l_gamma(p, one)
res = one - exp(s1 + log(gpx(p,x)))
else if(x <= real(p, ${k2}$)) then
s1 = -x + p * log(abs(x)) - l_gamma(p, one)
res = exp(log(gpx(p, x)) + s1)
end if
end function regamma_p_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
#:for k1, t1 in R_KINDS_TYPES
impure elemental function regamma_q_${t1[0]}$${k1}$(p, x) result(res)
!
! Approximation of regularized incomplete gamma function Q(p,x) for real p
!
${t1}$, intent(in) :: p, x
${t1}$ :: res, s1
${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_q" &
//" function is not defined at x < 0")
if(x == zero) then
res = one
else if(x > p) then
s1 = -x + p * log(x) - log_gamma(p)
res = exp(s1 + log(gpx(p,x)))
else if(x <= p) then
s1 = -x + p * log(abs(x)) - log_gamma(p)
res = one - exp(log(gpx(p, x)) + s1)
end if
end function regamma_q_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:for k2, t2 in R_KINDS_TYPES
impure elemental function regamma_q_${t1[0]}$${k1}$${k2}$(p, x) result(res)
!
! Approximation of regularized incomplet gamma function Q(p,x) for integer p
!
${t1}$, intent(in) :: p
${t2}$, intent(in) :: x
${t2}$ :: res, s1
${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$
if(x < zero) call error_stop("Error(regamma_q): Regularized gamma_q" &
//" function is not defined at x < 0")
if(x == zero) then
res = one
else if(x > real(p, ${k2}$)) then
s1 = -x + p * log(x) - l_gamma(p, one)
res = exp(log(gpx(p,x)) + s1)
elseif(x <= real(p, ${k2}$)) then
s1 = -x + p * log(abs(x)) - l_gamma(p, one)
res = one - exp(s1 + log(gpx(p,x)))
end if
end function regamma_q_${t1[0]}$${k1}$${k2}$
#:endfor
#:endfor
end module stdlib_specialfunctions_gamma
