#:set WITH_QP = False #:set WITH_XDP = False #:include "common.fypp" #:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_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 REAL_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 REAL_KINDS_TYPES module procedure ingamma_low_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor #:for k1, t1 in REAL_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 REAL_KINDS_TYPES module procedure l_ingamma_low_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor #:for k1, t1 in REAL_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 REAL_KINDS_TYPES module procedure ingamma_up_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor #:for k1, t1 in REAL_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 REAL_KINDS_TYPES module procedure l_ingamma_up_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor #:for k1, t1 in REAL_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 REAL_KINDS_TYPES module procedure regamma_p_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor #:for k1, t1 in REAL_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 REAL_KINDS_TYPES module procedure regamma_q_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor #:for k1, t1 in REAL_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 REAL_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 REAL_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 REAL_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 CMPLX_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 :: zero = 0.0_${k2}$, 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 REAL_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 CMPLX_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 :: res ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$ real, parameter :: zero_k2 = 0.0 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_k2 case (one) res = zero_k2 case (two : ) res = l_gamma(n + 1, 1.0D0) end select end function l_factorial_${t1[0]}$${k1}$ #:endfor #:for k1, t1 in REAL_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 REAL_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, m ${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 REAL_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 REAL_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 REAL_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 REAL_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 REAL_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 REAL_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 REAL_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 REAL_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 REAL_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 REAL_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 REAL_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 REAL_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