From f4caf5d0d1af65c0b725cc44a30517686d74d9c7 Mon Sep 17 00:00:00 2001 From: jim-215-fisher Date: Mon, 17 Jan 2022 21:27:27 -0500 Subject: [PATCH 01/15] inital draft --- src/stdlib_specialfunctions_gamma.fypp | 1456 ++++++++++++++++++++++++ 1 file changed, 1456 insertions(+) create mode 100644 src/stdlib_specialfunctions_gamma.fypp diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp new file mode 100644 index 000000000..b2a95221e --- /dev/null +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -0,0 +1,1456 @@ +#:include "common.fypp" +#:set ALL_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +module specialfunctions_gamma + use stdlib_kinds, only : sp, dp, xdp, qp, int8, int16, int32, int64 + use stdlib_error, only : error_stop + + implicit none + private + + integer(int8), parameter :: max_fact_int8 = 5_int8 + integer(int16), parameter :: max_fact_int16 = 7_int16 + integer(int32), parameter :: max_fact_int32 = 12_int32 + integer(int64), parameter :: max_fact_int64 = 20_int64 + + #:for k1, t1 in REAL_KINDS_TYPES + ${t1}$, parameter :: tol_${k1}$ = epsilon(1.0_${k1}$) + #:endfor + + + + 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 + public :: inverse_regularized_gamma_p + + + + interface gamma + !! Gamma function for integer, real, and complex numbers + !! + #:for k1, t1 in ALL_KINDS_TYPES + module procedure gamma_${t1[0]}$${k1}$ + #:endfor + end interface gamma + + + + interface log_gamma + !! Logarithm of gamma function + !! + #:for k1, t1 in ALL_KINDS_TYPES + module procedure l_gamma_${t1[0]}$${k1}$ !1 dummy + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure l_gamma_${t1[0]}$${k1}$${k2}$ !2 dummies + #:endfor + #: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}$ !1 dummy + #:endfor + + #: for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure l_factorial_${t1[0]}$${k1}$${k2}$ !2 dummies + #:endfor + #:endfor + end interface log_factorial + + + + interface gpx + ! Evaluation of incomplete gamma function + ! + #:for k1, t1 in REAL_KINDS_TYPES + module procedure gpx_${t1[0]}$${k1}$ !for real s 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 s and real x + #:endfor + #:endfor + end interface gpx + + + + 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 inverse_regularized_gamma_p + !! Inverse regularized lower incomplete gamma function + !! + #:for k1, t1 in REAL_KINDS_TYPES + module procedure inv_regamma_p_${t1[0]}$${k1}$ + #:endfor + end interface inverse_regularized_gamma_p + + + + + + +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 REAL_KINDS_TYPES + #:if k1 == "sp" + #:set k2 = "dp" + #:else + #:set k2 = "qp" + #:endif + #:set t2 = "real({})".format(k2) + + impure elemental function gamma_${t1[0]}$${k1}$(z) result(res) + ! + ! Gamma function for any real number, excluding negative whole number. + ! Algorithm is based on Glendon Pugh, "An Analysis of The Lanczos Gamma + ! Approximation", The University of British Columbia, 2004 + ! + ! Fortran 90 program by Jim-215-Fisher + ! + ${t1}$, intent(in) :: z + ${t1}$ :: res + integer :: i + + ${t1}$, 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) + ${t2}$ :: 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. + + #:else + #! 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(z <= zero_k1 .and. abs(anint(z) - z) < tol_${k1}$) & + call error_stop("Error(gamma): Gamma function real argument must " & + //"NOT be negative whole number") + + y = abs(z) - 1 + x = y + half + sum = d(0) + + do i = 1, n + + sum = sum + d(i) / (y + i) + + end do + + x = exp(x * log(x + r) - y) * sum + x = x * 2 / sqpi !positive z return + + if(z < zero_k1) then + + x = - pi / (sin(pi * (y + 1)) * (y + 1) * x) !negative z return + + end if + + res = x + end function gamma_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + #:if k1 == "sp" + #:set k2 = "dp" + #:else + #: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. + + #:else + #! 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 + + y = z - one + + else + + x = cmplx(abs(z % re), - z % im, kind = ${k1}$) + y = x - 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 REAL_KINDS_TYPES + impure elemental function l_gamma_${t1[0]}$${k1}$(z) result (res) + ! + ! log_gamma function for absolute value of any real number excluding + ! negative whole number. + ! + ${t1}$, intent(in) :: z + ${t1}$ :: res, y + ${t1}$, parameter :: zero = 0.0_${k1}$ + + + if(z <= zero .and. abs(anint(z) - z) < tol_${k1}$) call error_stop( & + "Error(log_gamma): Gamma function real argument must NOT be" & + //" negative whole number") + + res = log(abs(gamma(z))) + end function l_gamma_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function l_gamma_${t1[0]}$${k1}$(z) result (res) + ! + ! log_gamma function for any complex number, excluding negative whole number + ! + ${t1}$, intent(in) :: z + ${t1}$ :: res + + res = log(gamma(z)) + 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!) with single precision result, n is integer + ! + ${t1}$, intent(in) :: n + real :: res + ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$ + + if(n < 0_${k1}$) call error_stop("Error(l_factorial): Logarithm of " & + //"factorial function argument must be non-negative") + + select case(n) + + case (zero) + + res = 0.0 + + case (one) + + res = 0.0 + + case (two : ) + + res = log_gamma(real(n + 1)) + + end select + end function l_factorial_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + impure elemental function l_factorial_${t1[0]}$${k1}$${k2}$(n,x) result(res) + ! + ! Log(n!) with defined prescision output, n is integer, x is a real + ! for specified kind + ! + ${t1}$, intent(in) :: n + ${t2}$, intent(in) :: x + ${t2}$ :: res + ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$ + ${t2}$, parameter :: zero_k2 = 0.0_${k2}$ + + 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 = log_gamma(real(n + 1, ${k2}$)) + + end select + end function l_factorial_${t1[0]}$${k1}$${k2}$ + + #:endfor + #: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 with negative x") + + 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") + + if(p >= p_lim) then !use modified Lentz method for 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}$) + + elseif(x >= -9.0_${k2}$ .and. x <= zero) then + + p_lim = zero + + else + + p_lim = x + + endif + + 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 + log_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 < zero) then + + s1 = -x + p * log(abs(x)) + res = gpx(p, x) * exp(s1) + res = (-1) ** nint(p) * res + + else if(x <= p) then + + s1 = -x + p * log(x) + res = gpx(p, x) * exp(s1) + + else + + s1 = log_gamma(p) + y = one - exp(-x + p * log(x) - s1) * gpx(p, x) + res = exp(s1 + log(y)) + + endif + 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 < zero) then + + s1 = -x + p * log(abs(x)) + res = gpx(p, x) * exp(s1) + res = (-1) ** p * res + + else if(x <= real(p, ${k2}$)) then + + s1 = -x + p * log(x) + res = gpx(p, x) * exp(s1) + + else + + s1 = log_gamma(p, one) + y = one - exp(-x + p * log(x) - s1) * gpx(p, x) + res = exp(s1 + log(y)) + + endif + 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 = -x + p * log(abs(x)) + res = log(abs(gpx(p, x))) + s1 + + else if(x > p) then + + s1 = log_gamma(p) + y = one - exp(-x + p * log(x) - s1) * gpx(p, x) + res = s1 + log(y) + + endif + 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 = -x + p * log(abs(x)) + res = log(abs(gpx(p, x))) + s1 + + else if(x > real(p, ${k2}$)) then + + s1 = log_gamma(p, one) + y = one - exp(-x + p * log(x) - s1) * gpx(p, x) + res = s1 + log(y) + + endif + 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 < zero) then + + y = log_gamma(p) + s1 = -x + p * log(abs(x)) - y + res = gpx(p, x) * exp(s1) + res = (one - (-1) ** nint(p) * res) * exp(y) + + else if(x <= p) then + + y = log_gamma(p) + s1 = -x + p * log(x) - y + res = (one - gpx(p, x) * exp(s1)) * exp(y) + + else + + s1 = -x + p * log(x) + res = gpx(p, x) * exp(s1) + + endif + 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 < zero) then + + y = log_gamma(p, one) + s1 = -x + p * log(abs(x)) - y + res = gpx(p, x) * exp(s1) + res = (one - (-1) ** p * res) * exp(y) + + else if(x <= real(p, ${k2}$)) then + + y = log_gamma(p, one) + s1 = -x + p * log(x) - y + res = gpx(p, x) * exp(s1) + res = (one - res) * exp(y) + + else + + s1 = -x + p * log(x) + res = gpx(p, x) * exp(s1) + + endif + 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 < zero) then + + y = log_gamma(p) + s1 = -x + p * log(abs(x)) - y + res = gpx(p, x) * exp(s1) + res = log(one - (-1) ** nint(p) * res) + y + + else if(x <= p) then + + y= log_gamma(p) + s1 = -x + p * log(x) - y + res = gpx(p, x) * exp(s1) + res = log(one - res) + y + + else + + s1 = -x + p * log(x) + res = log(gpx(p, x)) + s1 + + endif + 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 = log_gamma(p, one) + + else if(x < zero) then + + y = log_gamma(p, one) + s1 = -x + p * log(abs(x)) - y + res = gpx(p, x) * exp(s1) + res = log(one - (-1) ** p * res) + y + + else if(x <= real(p, ${k2}$)) then + + y = log_gamma(p, one) + s1 = -x + p * log(x) - y + res = gpx(p, x) * exp(s1) + res = log(one - res) + y + + else + + s1 = -x + p * log(x) + res = log(gpx(p, x)) + s1 + + endif + 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") + + s1 = -x + p * log(x) - log_gamma(p) + + if(x == zero) then + + res = zero + + else if(x > zero .and. x <= p) then + + res = exp(log(gpx(p, x)) + s1) + + else if(x > p) then + + res = one - exp(s1 + log(gpx(p,x))) + + endif + 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") + + s1 = -x + p * log(x) - log_gamma(p, one) + + if(x == zero) then + + res = zero + + else if(x > zero .and. x <= real(p, ${k2}$)) then + + res = exp(log(gpx(p, x)) + s1) + + else if(x > real(p, ${k2}$)) then + + res = one - exp(s1 + log(gpx(p,x))) + + endif + 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") + + s1 = -x + p * log(x) - log_gamma(p) + + if(x == zero) then + + res = one + + else if(x > zero .and. x <= p) then + + res = one - exp(log(gpx(p, x)) + s1) + + else if(x > p) then + + res = exp(s1 + log(gpx(p,x))) + + endif + 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") + + s1 = -x + p * log(x) - log_gamma(p, one) + + if(x == zero) then + + res = one + + elseif(x > zero .and. x <= real(p, ${k2}$)) then + + res = one - exp(s1 + log(gpx(p,x))) + + elseif(x > real(p, ${k2}$)) then + + res = exp(log(gpx(p,x)) + s1) + + endif + end function regamma_q_${t1[0]}$${k1}$${k2}$ + + #:endfor + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function inv_regamma_p_${t1[0]}$${k1}$(s, p) result(res) + ! + ! Return x such that p = P(s, x) for p between 0 and 1 + ! + ! Numerical Recipes, William H. Press, etc., 3rd ed., 2007, p.263 + ! + ! "Modified Halley's Method for Solving Nonlinear Functions with Convergence + ! of Order Six and Efficicency Index 1.8171", Waqas Nazeer, etc., + ! International Journal of Pure and Applied Mathematics, 111(1), 2016, p. 55 + ! + ${t1}$, intent(in) :: s + real, intent(in) :: p + ${t1}$ :: res + integer :: i + real :: pp + ${t1}$ :: x0, x0s, x1, x2, t, err, f1, f2 + ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ + + if(s <= zero) call error_stop("Error(inv_regamma_p): Inverse " & + //"regularized gamma_p function input s value must be positive") + + if(p >= 1.0) then + + res = max(100._${k1}$, s + 100._${k1}$ * sqrt(s)) + return + + elseif(p <= 0.) then + + res = zero + return + + end if + + if(s > one) then + + pp = p + + if(p > 0.5) pp = 1.0 - p + + t = sqrt(-2 * log(pp)) + x0 = (2.30753_${k1}$ + t * 0.27061_${k1}$) / (one + t * & + (0.99229_${k1}$ + t * 0.04481_${k1}$)) - t + + if(p < 0.5) x0 = - x0 + + x0 = max( 1.0e-3_${k1}$, s * (one - one / (s * 9._${k1}$) - x0 / & + (3._${k1}$ * sqrt(s))) ** 3 ) + + else + + t = one - s * (0.253_${k1}$ + s * 0.12_${k1}$) + + if(real(p, ${k1}$) < t) then + + x0 = (p / t) ** ( one / s ) + + else + + x0 = one - log(one - (p - t) / (one - t)) + + end if + + end if + + err = regularized_gamma_p(s, x0) - p + x0s = x0 + x2 = (x0 + x0s) / 2 + f1 = exp(-x2 + (s - 1) * log(x2) - log_gamma(x2)) + f2 = err / (f1 - err * ((s - 1) / x2 - 1) / 2) + x1 = x0 - f2 + + do + + if(abs(f2) <= tol_${k1}$) exit + + x0 = x1 + err = regularized_gamma_p(s, x0) - p + x0s = x0 - err / (f1 - err * ((s - 1) / x2 - 1) / 2) + x2 = (x0 + x0s) / 2 + f1 = exp(-x2 + (s - 1) * log(x2) - log_gamma(x2)) + f2 = err / (f1 - err * ((s - 1) / x2 - 1) / 2) + x1 = x0 - f2 + + end do + + res = x1 + end function inv_regamma_p_${t1[0]}$${k1}$ + + #:endfor + +end module specialfunctions_gamma From 50baa93948f25af9de73f9ba22ebc43b3bae3df0 Mon Sep 17 00:00:00 2001 From: jim-215-fisher Date: Thu, 27 Jan 2022 10:49:50 -0500 Subject: [PATCH 02/15] Add draft specialfunction_gamma --- doc/specs/stdlib_specialfunctions_gamma.md | 306 ++++++++++ src/CMakeLists.txt | 1 + src/Makefile.manual | 4 + src/stdlib_specialfunctions_gamma.fypp | 616 ++++++++------------- src/tests/specialfunctions/CMakeLists.txt | 10 + src/tests/specialfunctions/Makefile.manual | 9 + src/tests/stats/CMakeLists.txt | 28 +- 7 files changed, 555 insertions(+), 419 deletions(-) create mode 100644 doc/specs/stdlib_specialfunctions_gamma.md create mode 100644 src/tests/specialfunctions/CMakeLists.txt create mode 100644 src/tests/specialfunctions/Makefile.manual diff --git a/doc/specs/stdlib_specialfunctions_gamma.md b/doc/specs/stdlib_specialfunctions_gamma.md new file mode 100644 index 000000000..9233d60a0 --- /dev/null +++ b/doc/specs/stdlib_specialfunctions_gamma.md @@ -0,0 +1,306 @@ +--- +title: specialfunctions_gamma +--- + +# Special functions gamma + +[TOC] + +## `gamma` - Calculate gamma function with any number + +### Status + +Experimental + +### Description + +Intrinsic gamma function provides values for real type argument with single and double precision. Here the gamma function is extended to both integer and complex numbers. + +### Syntax + +`result = [[stdlib_specialfunctions_gamma(module):gamma(interface)]] (x)` + +### Class + +Elemental function + +### Arguments + +`x`: should be an integer or a complex type number + +### Return value + +The function returns a value with the same type and kind as input argument. + + +## `log_gamma` - calculate logarithm gamma function with any number + +### Status + +Experimental + +### Description + +Intrinsic log_gamma function provides values for real type arguments with single and double precisions. Here the log_gamma function is extended to both integer and complex numbers. + +### Syntax + +`result = [[stdlib_specialfunctions_gamma(module):log_gamma(interface)]] (x)` + +### Class + +Elemental function + +### Arguments + +`x`: Shall be an integer or a complex type number. + +### Return value + +The function returns a value with the same type and kind as input argument. For integer argument, the result is single precision real type. + +### Example + +```fortran +program demo_log_gamma + use stdlib_specialfunctions_gamma, only : lg => log_gamma + implicit none + + integer :: i + real :: x + real(real64) :: y + complex :: z + + i = 10 + x = 8.76 + y = x + z = (5.345, -3.467) + print *, lg(i) !default single precision output +!12.8018274 + + print *, lg(x) !same kind as input + +!10.0942659 + + print *, lg(y) !same kind as input + +!10.094265528673880 + + print *, lg(z) !same kind as input + +!(2.56165719, 0.549360633) + +end program demo_log_gamma +``` + +## `log_factorial` - calculate logarithm of a factorial + +### Status + +Experimental + +### Description + +Compute the logarithm of factorial, log(n!) + +### Syntax + +`result = [[stdlib_specialfunctions_gamma(module):log_factorial(interface)]] (x)` + +### Class + +Elemental function + +### Arguments + +`x`: Shall be an integer type number. + +### Return value + +The function returns a value with single precision real type. + + +## `lower_incomplete_gamma` - calculate lower incomplete gamma integral + +### Status + +Experimental + +### Description + +The lower incomplete gamma function is defined as: + +\gamma (p, x) = \int_{0}^{x}t^{p-1}e^{-t}dt, \; \; p >0,\; x,p\in \mathbb{R} + +### Syntax + +`result = [[stdlib_specialfunctions_gamma(module):lower_incomplete_gamma(interface)]] (p, x)` + +### Class + +Elemental function + +### Arguments + +`p`: is a positive integer or real type argument. + +`x`: is a real type argument. + +### Return value + +The function returns a real type value with the same kind as argument x. + + +## `upper_incomplete_gamma` - calculate upper incomplete gamma integral + +### Status + +Experimental + +### Description + +The upper incomplete gamma function is defined as: + +\\Gamma (p, x) = \int_{x}^{\infty }t^{p-1}e^{-t}dt, \; \; p >0,\; x,p\in \mathbb{R} + +### Syntax + +`result = [[stdlib_specialfunctions_gamma(module):upper_incomplete_gamma(interface)]] (p, x)` + +### Class + +Elemental function + +### Arguments + +`p`: is a positive integer or real type argument. + +`x`: is a real type argument. + +### Return value + +The function returns a real type value with the same kind as argument x. + + +## `log_lower_incomplete_gamma` - calculate logarithm of the lower incomplete gamma integral + +### Status + +Experimental + +### Description + +Compute the logarithm of the absolute value of the lower incomplete gamma function. + +### Syntax + +`result = [[stdlib_specialfunctions_gamma(module):log_lower_incomplete_gamma(interface)]] (p, x)` + +### Class + +Elemental function + +### Arguments + +`p`: is a positive integer or real type argument. + +`x`: is a real type argument. + +### Return value + +The function returns a real type value with the same kind as argument x. + + +## `log_upper_incomplete_gamma` - calculate logarithm of the upper incomplete gamma integral + +### Status + +Experimental + +### Description + +Compute the logarithm of the absolute value of the upper incomplete gamma function. + +### Syntax + +`result = [[stdlib_specialfunctions_gamma(module):log_upper_incomplete_gamma(interface)]] (p, x)` + +### Class + +Elemental function + +### Arguments + +`p`: is a positive integer or real type argument. + +`x`: is a real type argument. + +### Return value + +The function returns a real type value with the same kind as argument x. + + +## `regularized_gamma_p` - calculate the gamma quotient P + +### Status + +Experimental + +### Description + +The regularized gamma quotient p, also known as normalized incomplete gamma function, is defined as: + +P(p,x)=\gamma(p,x)/\Gamma(p) + +The values of regularized gamma P is in the range of [0, 1] + +### Syntax + +`result = [[stdlib_specialfunctions_gamma(module):regularized_gamma_p(interface)]] (p, x)` + +### Class + +Elemental function + +### Arguments + +`p`: is a positive integer or real type argument. + +`x`: is a real type argument. + +### Return value + +The function returns a real type value with the same kind as argument x. + + +## `regularized_gamma_q` - calculate the gamma quotient Q + +### Status + +Experimental + +### Description + +The regularized gamma quotient Q is defined as: + +Q(p,x)=\Gamma(p,x)/\Gamma(p)=1-P(p,x) + +The values of regularized gamma Q is in the range of [0, 1] + +### Syntax + +`result = [[stdlib_specialfunctions_gamma(module):regularized_gamma_q(interface)]] (p, x)` + +### Class + +Elemental function + +### Arguments + +`p`: is a positive integer or real type argument. + +`x`: is a real type argument. + +### Return value + +The function returns a real type value with the same kind as argument x. diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f1f9051dd..4e8de17dc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -28,6 +28,7 @@ set(fppFiles stdlib_sorting_ord_sort.fypp stdlib_sorting_sort.fypp stdlib_sorting_sort_index.fypp + stdlib_specialfunctions_gamma.fypp stdlib_stats.fypp stdlib_stats_corr.fypp stdlib_stats_cov.fypp diff --git a/src/Makefile.manual b/src/Makefile.manual index f61a6fef6..7b9ac2d80 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -30,6 +30,7 @@ SRCFYPP = \ stdlib_sorting_ord_sort.fypp \ stdlib_sorting_sort.fypp \ stdlib_sorting_sort_index.fypp \ + stdlib_specialfunctions_gamma.fypp \ stdlib_stats.fypp \ stdlib_stats_corr.fypp \ stdlib_stats_cov.fypp \ @@ -165,6 +166,9 @@ stdlib_sorting_sort.o: \ stdlib_sorting.o stdlib_sorting_sort_index.o: \ stdlib_sorting.o +stdlib_specialfunctions_gamma.o: \ + stdlib_kinds.o \ + stdlib_error.o stdlib_stats.o: \ stdlib_kinds.o stdlib_stats_corr.o: \ diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp index b2a95221e..5047aea50 100644 --- a/src/stdlib_specialfunctions_gamma.fypp +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -1,5 +1,7 @@ +#:set WITH_QP = False +#:set WITH_XDP = False #:include "common.fypp" -#:set ALL_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +#:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_KINDS_TYPES module specialfunctions_gamma use stdlib_kinds, only : sp, dp, xdp, qp, int8, int16, int32, int64 use stdlib_error, only : error_stop @@ -7,14 +9,15 @@ module specialfunctions_gamma implicit none private - integer(int8), parameter :: max_fact_int8 = 5_int8 - integer(int16), parameter :: max_fact_int16 = 7_int16 - integer(int32), parameter :: max_fact_int32 = 12_int32 - integer(int64), parameter :: max_fact_int64 = 20_int64 + 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) @@ -22,14 +25,13 @@ module specialfunctions_gamma public :: lower_incomplete_gamma, log_lower_incomplete_gamma public :: upper_incomplete_gamma, log_upper_incomplete_gamma public :: regularized_gamma_p, regularized_gamma_q - public :: inverse_regularized_gamma_p interface gamma - !! Gamma function for integer, real, and complex numbers + !! Gamma function for integer and complex numbers !! - #:for k1, t1 in ALL_KINDS_TYPES + #:for k1, t1 in CI_KINDS_TYPES module procedure gamma_${t1[0]}$${k1}$ #:endfor end interface gamma @@ -39,14 +41,8 @@ module specialfunctions_gamma interface log_gamma !! Logarithm of gamma function !! - #:for k1, t1 in ALL_KINDS_TYPES - module procedure l_gamma_${t1[0]}$${k1}$ !1 dummy - #:endfor - - #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in REAL_KINDS_TYPES - module procedure l_gamma_${t1[0]}$${k1}$${k2}$ !2 dummies - #:endfor + #:for k1, t1 in CI_KINDS_TYPES + module procedure l_gamma_${t1[0]}$${k1}$ #:endfor end interface log_gamma @@ -56,34 +52,12 @@ module specialfunctions_gamma !! Logarithm of factorial n!, integer variable !! #:for k1, t1 in INT_KINDS_TYPES - module procedure l_factorial_${t1[0]}$${k1}$ !1 dummy - #:endfor - - #: for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in REAL_KINDS_TYPES - module procedure l_factorial_${t1[0]}$${k1}$${k2}$ !2 dummies - #:endfor + module procedure l_factorial_${t1[0]}$${k1}$ #:endfor end interface log_factorial - interface gpx - ! Evaluation of incomplete gamma function - ! - #:for k1, t1 in REAL_KINDS_TYPES - module procedure gpx_${t1[0]}$${k1}$ !for real s 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 s and real x - #:endfor - #:endfor - end interface gpx - - - interface lower_incomplete_gamma !! Lower incomplete gamma function !! @@ -180,16 +154,35 @@ module specialfunctions_gamma - interface inverse_regularized_gamma_p - !! Inverse regularized lower incomplete gamma function - !! + interface gpx + ! Incomplete gamma G function. + ! Internal use only + ! #:for k1, t1 in REAL_KINDS_TYPES - module procedure inv_regamma_p_${t1[0]}$${k1}$ + 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 inverse_regularized_gamma_p + 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 + + @@ -223,119 +216,11 @@ contains - #:for k1, t1 in REAL_KINDS_TYPES - #:if k1 == "sp" - #:set k2 = "dp" - #:else - #:set k2 = "qp" - #:endif - #:set t2 = "real({})".format(k2) - - impure elemental function gamma_${t1[0]}$${k1}$(z) result(res) - ! - ! Gamma function for any real number, excluding negative whole number. - ! Algorithm is based on Glendon Pugh, "An Analysis of The Lanczos Gamma - ! Approximation", The University of British Columbia, 2004 - ! - ! Fortran 90 program by Jim-215-Fisher - ! - ${t1}$, intent(in) :: z - ${t1}$ :: res - integer :: i - - ${t1}$, 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) - ${t2}$ :: 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. - - #:else - #! 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(z <= zero_k1 .and. abs(anint(z) - z) < tol_${k1}$) & - call error_stop("Error(gamma): Gamma function real argument must " & - //"NOT be negative whole number") - - y = abs(z) - 1 - x = y + half - sum = d(0) - - do i = 1, n - - sum = sum + d(i) / (y + i) - - end do - - x = exp(x * log(x + r) - y) * sum - x = x * 2 / sqpi !positive z return - - if(z < zero_k1) then - - x = - pi / (sin(pi * (y + 1)) * (y + 1) * x) !negative z return - - end if - - res = x - end function gamma_${t1[0]}$${k1}$ - - #:endfor - - #:for k1, t1 in CMPLX_KINDS_TYPES #:if k1 == "sp" #:set k2 = "dp" - #:else + #:elif k1 == "dp" #:set k2 = "qp" #:endif #:set t2 = "real({})".format(k2) @@ -368,7 +253,7 @@ contains -2.71994908488607704e-9_${k2}$] ! parameters from above referenced source. - #:else + #:elif k1 == "dp" #! for double precision input, using quadruple precision for calculation integer, parameter :: n = 24 @@ -527,88 +412,123 @@ contains - #:for k1, t1 in REAL_KINDS_TYPES + + #: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 absolute value of any real number excluding - ! negative whole number. + ! log_gamma function for any complex number, excluding negative whole number + ! "Computation of special functions", Shanjie Zhang & Jianmin Jin, 1996, p.48 + ! + ! Fortran 90 program by Jim-215-Fisher ! ${t1}$, intent(in) :: z - ${t1}$ :: res, y - ${t1}$, parameter :: zero = 0.0_${k1}$ + ${t1}$ :: res, z1, z2 + real(${k1}$) :: d + integer :: m, i + complex(${k2}$) :: zr, zr2, sum, s + real(${k1}$), parameter :: z_limit = 7_${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 - if(z <= zero .and. abs(anint(z) - z) < tol_${k1}$) call error_stop( & - "Error(log_gamma): Gamma function real argument must NOT be" & - //" negative whole number") + if(z % re > zero_k1) then - res = log(abs(gamma(z))) - end function l_gamma_${t1[0]}$${k1}$ + z2 = z - #:endfor + else + z2 = cmplx(abs(z % re), - z % im, kind = ${k1}$) !for Re(z) < 0 + end if - #:for k1, t1 in CMPLX_KINDS_TYPES - impure elemental function l_gamma_${t1[0]}$${k1}$(z) result (res) - ! - ! log_gamma function for any complex number, excluding negative whole number - ! - ${t1}$, intent(in) :: z - ${t1}$ :: res + d = hypot(z2 % re, z2 % im) + z1 = z2 + m = 0 - res = log(gamma(z)) - end function l_gamma_${t1[0]}$${k1}$ + if(d <= z_limit) then !for small |z| - #:endfor + m = ceiling(z_limit - d) + z1 = z2 + m + end if + zr = one / z1 + zr2 = zr * zr - #:for k1, t1 in INT_KINDS_TYPES - impure elemental function l_factorial_${t1[0]}$${k1}$(n) result(res) - ! - ! Log(n!) with single precision result, n is integer - ! - ${t1}$, intent(in) :: n - real :: res - ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$ + 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(n < 0_${k1}$) call error_stop("Error(l_factorial): Logarithm of " & - //"factorial function argument must be non-negative") + if(m /= 0) then - select case(n) + s = cmplx(zero, zero, kind = ${k2}$) - case (zero) + do i = 1, m - res = 0.0 + s = s + log(cmplx(z1, kind = ${k2}$) - i) - case (one) + end do - res = 0.0 + sum = sum - s - case (two : ) + end if - res = log_gamma(real(n + 1)) + if(z % re < zero_k1) then - end select - end function l_factorial_${t1[0]}$${k1}$ + sum = -pi / (sin(pi * z2) * z2 * sum) !Re(Z) < 0 + + end if + + res = sum + 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_factorial_${t1[0]}$${k1}$${k2}$(n,x) result(res) + impure elemental function l_factorial_${t1[0]}$${k1}$(n) result(res) ! - ! Log(n!) with defined prescision output, n is integer, x is a real - ! for specified kind + ! Log(n!) ! ${t1}$, intent(in) :: n - ${t2}$, intent(in) :: x - ${t2}$ :: res + real :: res ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$ - ${t2}$, parameter :: zero_k2 = 0.0_${k2}$ + real, parameter :: zero_k2 = 0.0 if(n < zero) call error_stop("Error(l_factorial): Logarithm of" & //" factorial function argument must be non-negative") @@ -625,12 +545,11 @@ contains case (two : ) - res = log_gamma(real(n + 1, ${k2}$)) + res = l_gamma(n + 1, 1.0D0) end select - end function l_factorial_${t1[0]}$${k1}$${k2}$ + end function l_factorial_${t1[0]}$${k1}$ - #:endfor #:endfor @@ -662,7 +581,7 @@ contains ${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 with negative x") + //" function must have a positive parameter p") if(x < -9.0_${k1}$) then @@ -680,9 +599,9 @@ contains 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") + //"negative x must come with a whole number p not too small") - if(p >= p_lim) then !use modified Lentz method for continued fraction + if(p >= p_lim) then !use modified Lentz method of continued fraction !for eq. (15) in the above reference. a = one b = p @@ -805,7 +724,7 @@ contains p_lim = 5.0_${k2}$ * (sqrt(abs(x)) - 1.0_${k2}$) - elseif(x >= -9.0_${k2}$ .and. x <= zero) then + else if(x >= -9.0_${k2}$ .and. x <= zero) then p_lim = zero @@ -813,7 +732,7 @@ contains p_lim = x - endif + end if if(real(p, ${k2}$) >= p_lim) then @@ -907,7 +826,7 @@ contains if(y >= b * tol_${k2}$ .and. mod(p, two) /= zero_k1) & b = b + d * c / a - g = ((-1) ** p * exp(-a + log_gamma(p, one) - (p - 1) * log(a)) & + g = ((-1) ** p * exp(-a + l_gamma(p, one) - (p - 1) * log(a)) & + b ) / a end if @@ -933,24 +852,24 @@ contains res = zero - else if(x < zero) then + else if(x > p) then - s1 = -x + p * log(abs(x)) - res = gpx(p, x) * exp(s1) - res = (-1) ** nint(p) * res + s1 = log_gamma(p) + y = one - exp(-x + p * log(x) - s1) * gpx(p, x) + res = exp(s1 + log(y)) - else if(x <= p) then + else if(x <= p .and. x > zero) then s1 = -x + p * log(x) res = gpx(p, x) * exp(s1) else - s1 = log_gamma(p) - y = one - exp(-x + p * log(x) - s1) * gpx(p, x) - res = exp(s1 + log(y)) + s1 = -x + p * log(abs(x)) + res = gpx(p, x) * exp(s1) + res = (-1) ** nint(p) * res - endif + end if end function ingamma_low_${t1[0]}$${k1}$ #:endfor @@ -973,24 +892,24 @@ contains res = zero - else if(x < zero) then + else if(x > real(p, ${k2}$)) then - s1 = -x + p * log(abs(x)) - res = gpx(p, x) * exp(s1) - res = (-1) ** p * res + 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}$)) then + else if(x <= real(p, ${k2}$) .and. x > zero) then s1 = -x + p * log(x) res = gpx(p, x) * exp(s1) else - s1 = log_gamma(p, one) - y = one - exp(-x + p * log(x) - s1) * gpx(p, x) - res = exp(s1 + log(y)) + s1 = -x + p * log(abs(x)) + res = gpx(p, x) * exp(s1) + res = (-1) ** p * res - endif + end if end function ingamma_low_${t1[0]}$${k1}$${k2}$ #:endfor @@ -1010,18 +929,18 @@ contains res = zero - else if(x <= p) then - - s1 = -x + p * log(abs(x)) - res = log(abs(gpx(p, x))) + s1 - else if(x > p) then s1 = log_gamma(p) y = one - exp(-x + p * log(x) - s1) * gpx(p, x) res = s1 + log(y) - endif + else if(x <= p) then + + s1 = -x + p * log(abs(x)) + res = log(abs(gpx(p, x))) + s1 + + end if end function l_ingamma_low_${t1[0]}$${k1}$ #:endfor @@ -1043,18 +962,18 @@ contains res = zero - else if(x <= real(p, ${k2}$)) then - - s1 = -x + p * log(abs(x)) - res = log(abs(gpx(p, x))) + s1 - else if(x > real(p, ${k2}$)) then - s1 = log_gamma(p, one) + s1 = l_gamma(p, one) y = one - exp(-x + p * log(x) - s1) * gpx(p, x) res = s1 + log(y) - endif + 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 @@ -1075,14 +994,12 @@ contains res = gamma(p) - else if(x < zero) then + else if(x > p) then - y = log_gamma(p) - s1 = -x + p * log(abs(x)) - y + s1 = -x + p * log(x) res = gpx(p, x) * exp(s1) - res = (one - (-1) ** nint(p) * res) * exp(y) - else if(x <= p) then + else if(x <= p .and. x > zero) then y = log_gamma(p) s1 = -x + p * log(x) - y @@ -1090,10 +1007,12 @@ contains else - s1 = -x + p * log(x) + y = log_gamma(p) + s1 = -x + p * log(abs(x)) - y res = gpx(p, x) * exp(s1) + res = (one - (-1) ** nint(p) * res) * exp(y) - endif + end if end function ingamma_up_${t1[0]}$${k1}$ #:endfor @@ -1116,26 +1035,26 @@ contains res = gamma(real(p, ${k2}$)) - else if(x < zero) then + else if(x > real(p, ${k2}$)) then - y = log_gamma(p, one) - s1 = -x + p * log(abs(x)) - y + s1 = -x + p * log(x) res = gpx(p, x) * exp(s1) - res = (one - (-1) ** p * res) * exp(y) - else if(x <= real(p, ${k2}$)) then + else if(x <= real(p, ${k2}$) .and. x > zero) then - y = log_gamma(p, one) + y = l_gamma(p, one) s1 = -x + p * log(x) - y res = gpx(p, x) * exp(s1) res = (one - res) * exp(y) else - s1 = -x + p * log(x) + 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) - endif + end if end function ingamma_up_${t1[0]}$${k1}$${k2}$ #:endfor @@ -1155,14 +1074,12 @@ contains res = log_gamma(p) - else if(x < zero) then + else if(x > p) then - y = log_gamma(p) - s1 = -x + p * log(abs(x)) - y - res = gpx(p, x) * exp(s1) - res = log(one - (-1) ** nint(p) * res) + y + s1 = -x + p * log(x) + res = log(gpx(p, x)) + s1 - else if(x <= p) then + else if(x <= p .and. x > zero) then y= log_gamma(p) s1 = -x + p * log(x) - y @@ -1171,10 +1088,12 @@ contains else - s1 = -x + p * log(x) - res = log(gpx(p, x)) + s1 + y = log_gamma(p) + s1 = -x + p * log(abs(x)) - y + res = gpx(p, x) * exp(s1) + res = log(one - (-1) ** nint(p) * res) + y - endif + end if end function l_ingamma_up_${t1[0]}$${k1}$ #:endfor @@ -1194,28 +1113,28 @@ contains if(x == zero) then - res = log_gamma(p, one) + res = l_gamma(p, one) - else if(x < zero) then + else if(x > real(p, ${k2}$)) then - y = log_gamma(p, one) - s1 = -x + p * log(abs(x)) - y - res = gpx(p, x) * exp(s1) - res = log(one - (-1) ** p * res) + y + s1 = -x + p * log(x) + res = log(gpx(p, x)) + s1 - else if(x <= real(p, ${k2}$)) then + else if(x <= real(p, ${k2}$) .and. x > zero) then - y = log_gamma(p, one) + y = l_gamma(p, one) s1 = -x + p * log(x) - y res = gpx(p, x) * exp(s1) res = log(one - res) + y else - s1 = -x + p * log(x) - res = log(gpx(p, x)) + s1 + y = l_gamma(p, one) + s1 = -x + p * log(abs(x)) - y + res = gpx(p, x) * exp(s1) + res = log(one - (-1) ** p * res) + y - endif + end if end function l_ingamma_up_${t1[0]}$${k1}$${k2}$ #:endfor @@ -1236,21 +1155,22 @@ contains if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_p" & //" function is not defined at x < 0") - s1 = -x + p * log(x) - log_gamma(p) if(x == zero) then res = zero - else if(x > zero .and. x <= p) then - - res = exp(log(gpx(p, x)) + s1) - else if(x > p) then + s1 = -x + p * log(x) - log_gamma(p) res = one - exp(s1 + log(gpx(p,x))) - endif + 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 @@ -1271,21 +1191,22 @@ contains if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_p" & //" function is not defined at x < 0") - s1 = -x + p * log(x) - log_gamma(p, one) if(x == zero) then res = zero - else if(x > zero .and. x <= real(p, ${k2}$)) then - - res = exp(log(gpx(p, x)) + s1) - else if(x > real(p, ${k2}$)) then + s1 = -x + p * log(x) - l_gamma(p, one) res = one - exp(s1 + log(gpx(p,x))) - endif + 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 @@ -1305,21 +1226,22 @@ contains if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_q" & //" function is not defined at x < 0") - s1 = -x + p * log(x) - log_gamma(p) if(x == zero) then res = one - else if(x > zero .and. x <= p) then - - res = one - exp(log(gpx(p, x)) + s1) - else if(x > p) then + s1 = -x + p * log(x) - log_gamma(p) res = exp(s1 + log(gpx(p,x))) - endif + 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 @@ -1340,117 +1262,25 @@ contains if(x < zero) call error_stop("Error(regamma_q): Regularized gamma_q" & //" function is not defined at x < 0") - s1 = -x + p * log(x) - log_gamma(p, one) if(x == zero) then res = one - elseif(x > zero .and. x <= real(p, ${k2}$)) then - - res = one - exp(s1 + log(gpx(p,x))) - - elseif(x > real(p, ${k2}$)) then + else if(x > real(p, ${k2}$)) then + s1 = -x + p * log(x) - l_gamma(p, one) res = exp(log(gpx(p,x)) + s1) - endif - end function regamma_q_${t1[0]}$${k1}$${k2}$ + elseif(x <= real(p, ${k2}$)) then - #:endfor - #:endfor - - - - #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function inv_regamma_p_${t1[0]}$${k1}$(s, p) result(res) - ! - ! Return x such that p = P(s, x) for p between 0 and 1 - ! - ! Numerical Recipes, William H. Press, etc., 3rd ed., 2007, p.263 - ! - ! "Modified Halley's Method for Solving Nonlinear Functions with Convergence - ! of Order Six and Efficicency Index 1.8171", Waqas Nazeer, etc., - ! International Journal of Pure and Applied Mathematics, 111(1), 2016, p. 55 - ! - ${t1}$, intent(in) :: s - real, intent(in) :: p - ${t1}$ :: res - integer :: i - real :: pp - ${t1}$ :: x0, x0s, x1, x2, t, err, f1, f2 - ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ - - if(s <= zero) call error_stop("Error(inv_regamma_p): Inverse " & - //"regularized gamma_p function input s value must be positive") - - if(p >= 1.0) then - - res = max(100._${k1}$, s + 100._${k1}$ * sqrt(s)) - return - - elseif(p <= 0.) then - - res = zero - return - - end if - - if(s > one) then - - pp = p - - if(p > 0.5) pp = 1.0 - p - - t = sqrt(-2 * log(pp)) - x0 = (2.30753_${k1}$ + t * 0.27061_${k1}$) / (one + t * & - (0.99229_${k1}$ + t * 0.04481_${k1}$)) - t - - if(p < 0.5) x0 = - x0 - - x0 = max( 1.0e-3_${k1}$, s * (one - one / (s * 9._${k1}$) - x0 / & - (3._${k1}$ * sqrt(s))) ** 3 ) - - else - - t = one - s * (0.253_${k1}$ + s * 0.12_${k1}$) - - if(real(p, ${k1}$) < t) then - - x0 = (p / t) ** ( one / s ) - - else - - x0 = one - log(one - (p - t) / (one - t)) - - end if + 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}$ - err = regularized_gamma_p(s, x0) - p - x0s = x0 - x2 = (x0 + x0s) / 2 - f1 = exp(-x2 + (s - 1) * log(x2) - log_gamma(x2)) - f2 = err / (f1 - err * ((s - 1) / x2 - 1) / 2) - x1 = x0 - f2 - - do - - if(abs(f2) <= tol_${k1}$) exit - - x0 = x1 - err = regularized_gamma_p(s, x0) - p - x0s = x0 - err / (f1 - err * ((s - 1) / x2 - 1) / 2) - x2 = (x0 + x0s) / 2 - f1 = exp(-x2 + (s - 1) * log(x2) - log_gamma(x2)) - f2 = err / (f1 - err * ((s - 1) / x2 - 1) / 2) - x1 = x0 - f2 - - end do - - res = x1 - end function inv_regamma_p_${t1[0]}$${k1}$ - + #:endfor #:endfor end module specialfunctions_gamma diff --git a/src/tests/specialfunctions/CMakeLists.txt b/src/tests/specialfunctions/CMakeLists.txt new file mode 100644 index 000000000..caa3a96b5 --- /dev/null +++ b/src/tests/specialfunctions/CMakeLists.txt @@ -0,0 +1,10 @@ +### Pre-process: .fpp -> .f90 via Fypp + +# Create a list of the files to be preprocessed +set(fppFiles + test_specialfunctions_gamma.fypp +) + +fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) + +ADDTEST(specialfunctions_gamma) diff --git a/src/tests/specialfunctions/Makefile.manual b/src/tests/specialfunctions/Makefile.manual new file mode 100644 index 000000000..790744555 --- /dev/null +++ b/src/tests/specialfunctions/Makefile.manual @@ -0,0 +1,9 @@ +SRCFYPP = \ + test_specialfunctions_gamma.fypp + +SRCGEN = $(SRCFYPP:.fypp=.f90) + +$(SRCGEN): %.f90: %.fypp ../../common.fypp + fypp -I../.. $(FYPPFLAGS) $< $@ + +include ../Makefile.manual.test.mk diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 651286708..caa3a96b5 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -2,33 +2,9 @@ # Create a list of the files to be preprocessed set(fppFiles - test_mean.fypp - test_mean_f03.fypp - test_median.fypp - test_distribution_uniform.fypp - test_distribution_normal.fypp - test_distribution_exponential.fypp + test_specialfunctions_gamma.fypp ) fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) -ADDTEST(corr) -ADDTEST(cov) -ADDTEST(mean) -ADDTEST(median) -ADDTEST(moment) -ADDTEST(rawmoment) -ADDTEST(var) -ADDTEST(varn) -ADDTEST(random) -ADDTEST(distribution_uniform) -ADDTEST(distribution_normal) -ADDTEST(distribution_exponential) - -if(DEFINED CMAKE_MAXIMUM_RANK) - if(${CMAKE_MAXIMUM_RANK} GREATER 7) - ADDTEST(mean_f03) - endif() -elseif(f03rank) - ADDTEST(mean_f03) -endif() +ADDTEST(specialfunctions_gamma) From c2507a4f87842171a1800e86d080e04d068d0b7e Mon Sep 17 00:00:00 2001 From: jim-215-fisher Date: Thu, 27 Jan 2022 10:56:49 -0500 Subject: [PATCH 03/15] CMakeList file correction --- src/tests/stats/CMakeLists.txt | 30 +++++++++++++++++++++++++++--- 1 file changed, 27 insertions(+), 3 deletions(-) diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index caa3a96b5..ff9d45063 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -1,10 +1,34 @@ -### Pre-process: .fpp -> .f90 via Fypp +#### Pre-process: .fpp -> .f90 via Fypp # Create a list of the files to be preprocessed set(fppFiles - test_specialfunctions_gamma.fypp + test_mean.fypp + test_mean_f03.fypp + test_median.fypp + test_distribution_uniform.fypp + test_distribution_normal.fypp + test_distribution_exponential.fypp ) fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) -ADDTEST(specialfunctions_gamma) +ADDTEST(corr) +ADDTEST(cov) +ADDTEST(mean) +ADDTEST(median) +ADDTEST(moment) +ADDTEST(rawmoment) +ADDTEST(var) +ADDTEST(varn) +ADDTEST(random) +ADDTEST(distribution_uniform) +ADDTEST(distribution_normal) +ADDTEST(distribution_exponential) + +if(DEFINED CMAKE_MAXIMUM_RANK) + if(${CMAKE_MAXIMUM_RANK} GREATER 7) + ADDTEST(mean_f03) + endif() +elseif(f03rank) + ADDTEST(mean_f03) +endif() From bae36c42840d7dfcedd556924f285d84b0246cbc Mon Sep 17 00:00:00 2001 From: jim-215-fisher Date: Wed, 2 Feb 2022 14:01:07 -0500 Subject: [PATCH 04/15] update gamma function --- doc/specs/stdlib_specialfunctions_gamma.md | 173 ++++++++++++++++++--- src/stdlib_specialfunctions_gamma.fypp | 56 +++---- src/tests/CMakeLists.txt | 1 + src/tests/Makefile.manual | 1 + 4 files changed, 183 insertions(+), 48 deletions(-) diff --git a/doc/specs/stdlib_specialfunctions_gamma.md b/doc/specs/stdlib_specialfunctions_gamma.md index 9233d60a0..da1e17436 100644 --- a/doc/specs/stdlib_specialfunctions_gamma.md +++ b/doc/specs/stdlib_specialfunctions_gamma.md @@ -6,7 +6,7 @@ title: specialfunctions_gamma [TOC] -## `gamma` - Calculate gamma function with any number +## `gamma` - Calculate gamma function with any type of argument ### Status @@ -14,7 +14,11 @@ Experimental ### Description -Intrinsic gamma function provides values for real type argument with single and double precision. Here the gamma function is extended to both integer and complex numbers. +The gamma function is defined as the analytic continuation of a convergent improper integral function on the whole complex plane except zero and negative integers: + +\Gamma(z)=\int_{0}^{\infty}x^{z-1}e^{-x}dx, \;\; z\in \mathbb{C} \setminus 0, -1, -2, \cdots + +Fortran 2018 standard implements intrinsic gamma function of real type argument in single and double precisions. Here the gamma function is extended to both integer and complex arguments. The values of gamma function with integer arguments are exact. The values of gamma function with complex arguments are approximated in single and double precisions by using Lanczos approxiamtion. ### Syntax @@ -26,14 +30,54 @@ Elemental function ### Arguments -`x`: should be an integer or a complex type number +`x`: should be a positive integer or a complex type number ### Return value The function returns a value with the same type and kind as input argument. +### Example +```fortran +program demo_gamma + use iso_fortran_env, only : real64, int64 + use stdlib_specialfunctions_gamma, only : gamma + implicit none + + integer :: i + integer(int64) :: n + real :: x + real(real64) :: y + complex :: z + complex(real64) :: z1 + + i = 10 + n = 15_int64 + x = 2.5 + y = 4.3_real64 + z = (2.3, 0.6) + z1 = (-4.2_real64, 3.1_real64) + + print *, gamma(i) !integer gives exact result +! 362880 + + print *, gamma(n) +! 87178291200 + + print *, gamma(x) ! intrinsic function call +! 1.32934034 -## `log_gamma` - calculate logarithm gamma function with any number + print *, gamma(y) ! intrinsic function call +! 8.8553433604540341 + + print *, gamma(z) +! (0.988054395, 0.383354813) + + print *, gamma(z1) +! (-2.78916032990983999E-005, 9.83164600163221218E-006) +end program demo_gamma +``` + +## `log_gamma` - calculate natural logarithm of gamma function with any type of argument ### Status @@ -41,7 +85,9 @@ Experimental ### Description -Intrinsic log_gamma function provides values for real type arguments with single and double precisions. Here the log_gamma function is extended to both integer and complex numbers. +Due to the different branch cut structures and a different principal branch, natural logarithm of gamma function log_gamma(z) with complex argument is different from the ln(Gamma(z)). The two have the same real part but different imaginary part. + +Fortran 2018 standard implements intrinsic log_gamma function of absolute value of real type argument in single and double precision. Here the log_gamma function is extended to both integer and complex arguments. The values of log_gamma function with complex arguments are approximated in single and double precisions by using Stirling's approximation. ### Syntax @@ -53,43 +99,49 @@ Elemental function ### Arguments -`x`: Shall be an integer or a complex type number. +`x`: Shall be a positive integer or a complex type number. ### Return value -The function returns a value with the same type and kind as input argument. For integer argument, the result is single precision real type. +The function returns real single precision values for integer input arguments, while it returns complex values with the same kind as complex input arguments. ### Example ```fortran program demo_log_gamma - use stdlib_specialfunctions_gamma, only : lg => log_gamma + use iso_fortran_env, only : real64 + use stdlib_specialfunctions_gamma, only : log_gamma implicit none integer :: i real :: x real(real64) :: y complex :: z + complex(real64) :: z1 i = 10 x = 8.76 y = x z = (5.345, -3.467) - print *, lg(i) !default single precision output + z1 = z + print *, log_gamma(i) !default single precision output !12.8018274 - print *, lg(x) !same kind as input + print *, log_gamma(x) !intrinsic function call !10.0942659 - print *, lg(y) !same kind as input + print *, log_gamma(y) !intrinsic function call !10.094265528673880 - print *, lg(z) !same kind as input + print *, log_gamma(z) !same kind as input -!(2.56165719, 0.549360633) +!(2.56165648, -5.73382425) + print *, log_gamma(z1) + +!(2.5616575105114614, -5.7338247782852498) end program demo_log_gamma ``` @@ -101,7 +153,7 @@ Experimental ### Description -Compute the logarithm of factorial, log(n!) +Compute the natural logarithm of factorial, log(n!) ### Syntax @@ -113,12 +165,30 @@ Elemental function ### Arguments -`x`: Shall be an integer type number. +`x`: Shall be a positive integer type number. ### Return value -The function returns a value with single precision real type. +The function returns real type values with single precision. +### Example +```fortran +program demo_log_factorial + use iso_fortran_env, only : int64 + use stdlib_specialfunctions_gamma, only : lf => log_factorial + implicit none + integer :: n + + n = 10 + print *, lf(n) + +! 15.1044130 + + print *, lf(35_int64) + +! 92.1361771 +end program demo_log_factorial +``` ## `lower_incomplete_gamma` - calculate lower incomplete gamma integral @@ -130,7 +200,9 @@ Experimental The lower incomplete gamma function is defined as: -\gamma (p, x) = \int_{0}^{x}t^{p-1}e^{-t}dt, \; \; p >0,\; x,p\in \mathbb{R} +\gamma(p,x)=\int_{0}^{x}t^{p-1}e^{-t}dt, \;\; p > 0, x\in \mathbb{R} + +When x < 0, p must be positive integer. ### Syntax @@ -150,6 +222,25 @@ Elemental function The function returns a real type value with the same kind as argument x. +### Example +```fortran +program demo_ligamma + use stdlib_specialfunctions_gamma, only : lig => lower_incomplete_gamma + implicit none + integer :: p + real :: p1, x + + p = 3 + p1 = 2.3 + print *, lig(p, -5.0) + +! -2521.02417 + + print *, lig(p1, 5.0) + +! 1.09715652 +end demo_ligamma +``` ## `upper_incomplete_gamma` - calculate upper incomplete gamma integral @@ -161,7 +252,9 @@ Experimental The upper incomplete gamma function is defined as: -\\Gamma (p, x) = \int_{x}^{\infty }t^{p-1}e^{-t}dt, \; \; p >0,\; x,p\in \mathbb{R} +\Gamma (p, x) = \int_{x}^{\infty }t^{p-1}e^{-t}dt, \; \; p >0,\; x \in \mathbb{R} + +When x < 0, p must be a positive integer. ### Syntax @@ -181,6 +274,21 @@ Elemental function The function returns a real type value with the same kind as argument x. +### Example +```fortran +program demo_uigamma + use stdlib_specialfunctions_gamma, only : uig => upper_incomplete_gamma + implicit none + + print *, uig(3, -5.0) + +!2523.02295 + + print *, uig(2.3, 5.0) + +!6.95552528E-02 +end program demo_uigamma +``` ## `log_lower_incomplete_gamma` - calculate logarithm of the lower incomplete gamma integral @@ -190,7 +298,7 @@ Experimental ### Description -Compute the logarithm of the absolute value of the lower incomplete gamma function. +Compute the natural logarithm of the absolute value of the lower incomplete gamma function. ### Syntax @@ -219,7 +327,7 @@ Experimental ### Description -Compute the logarithm of the absolute value of the upper incomplete gamma function. +Compute the natural logarithm of the absolute value of the upper incomplete gamma function. ### Syntax @@ -248,7 +356,7 @@ Experimental ### Description -The regularized gamma quotient p, also known as normalized incomplete gamma function, is defined as: +The regularized gamma quotient P, also known as normalized incomplete gamma function, is defined as: P(p,x)=\gamma(p,x)/\Gamma(p) @@ -272,6 +380,17 @@ Elemental function The function returns a real type value with the same kind as argument x. +### Example +```fortran +program demo_gamma_p + use stdlib_specialfunctions_gamma, only : rgp => regularized_gamma_p + implicit none + + print *, rgp(3.0, 5.0) + +! 0.875347972 +end program demo_gamma_p +``` ## `regularized_gamma_q` - calculate the gamma quotient Q @@ -304,3 +423,15 @@ Elemental function ### Return value The function returns a real type value with the same kind as argument x. + +### Example +```fortran +program demo_gamma_q + use stdlib_specialfunctions_gamma, only : rgq => regularized_gamma_q + implicit none + + print *, rgq(3.0, 5.0) + +! 0.124652028 +end program demo_gamma_q +``` diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp index 5047aea50..a8542c43c 100644 --- a/src/stdlib_specialfunctions_gamma.fypp +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -2,7 +2,7 @@ #:set WITH_XDP = False #:include "common.fypp" #:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_KINDS_TYPES -module specialfunctions_gamma +module stdlib_specialfunctions_gamma use stdlib_kinds, only : sp, dp, xdp, qp, int8, int16, int32, int64 use stdlib_error, only : error_stop @@ -424,6 +424,8 @@ contains ! ! 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 ! @@ -432,7 +434,7 @@ contains real(${k1}$) :: d integer :: m, i complex(${k2}$) :: zr, zr2, sum, s - real(${k1}$), parameter :: z_limit = 7_${k1}$, zero_k1 = 0.0_${k1}$ + 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) @@ -459,14 +461,11 @@ contains -.1236960214226927445425171034927132488108E+14_${k2}$] ! parameters from above reference + z2 = z - if(z % re > zero_k1) then - - z2 = z - - else + if(z % re < zero_k1) then - z2 = cmplx(abs(z % re), - z % im, kind = ${k1}$) !for Re(z) < 0 + z2 = cmplx(abs(z % re), - z % im, kind = ${k1}$) + 1 end if @@ -508,11 +507,13 @@ contains if(z % re < zero_k1) then - sum = -pi / (sin(pi * z2) * z2 * sum) !Re(Z) < 0 + 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 = sum + res = cmplx(sum, kind = ${k1}$) end function l_gamma_${t1[0]}$${k1}$ #:endfor @@ -865,9 +866,8 @@ contains else - s1 = -x + p * log(abs(x)) - res = gpx(p, x) * exp(s1) - res = (-1) ** nint(p) * res + 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}$ @@ -935,11 +935,16 @@ contains y = one - exp(-x + p * log(x) - s1) * gpx(p, x) res = s1 + log(y) - else if(x <= p) then + 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}$ @@ -1007,10 +1012,9 @@ contains else - y = log_gamma(p) - s1 = -x + p * log(abs(x)) - y - res = gpx(p, x) * exp(s1) - res = (one - (-1) ** nint(p) * res) * exp(y) + + 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}$ @@ -1088,10 +1092,8 @@ contains else - y = log_gamma(p) - s1 = -x + p * log(abs(x)) - y - res = gpx(p, x) * exp(s1) - res = log(one - (-1) ** nint(p) * res) + y + 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}$ @@ -1103,7 +1105,7 @@ contains #: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) & + impure elemental function l_ingamma_up_${t1[0]}$${k1}$${k2}$(p, x) & result(res) ${t1}$, intent(in) :: p @@ -1130,9 +1132,9 @@ contains else y = l_gamma(p, one) - s1 = -x + p * log(abs(x)) - y - res = gpx(p, x) * exp(s1) - res = log(one - (-1) ** p * res) + y + 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}$ @@ -1283,4 +1285,4 @@ contains #:endfor #:endfor -end module specialfunctions_gamma +end module stdlib_specialfunctions_gamma diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index de110ca62..a4250d7ba 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -26,6 +26,7 @@ add_subdirectory(logger) add_subdirectory(optval) add_subdirectory(selection) add_subdirectory(sorting) +add_subdirectory(specialfunctions) add_subdirectory(stats) add_subdirectory(string) add_subdirectory(system) diff --git a/src/tests/Makefile.manual b/src/tests/Makefile.manual index 370f2e3c8..90bedf3b3 100644 --- a/src/tests/Makefile.manual +++ b/src/tests/Makefile.manual @@ -23,6 +23,7 @@ all test clean:: $(MAKE) -f Makefile.manual --directory=selection $@ $(MAKE) -f Makefile.manual --directory=sorting $@ $(MAKE) -f Makefile.manual --directory=quadrature $@ + $(MAKE) -f Makefile.manual --directory=specialfunctions $@ $(MAKE) -f Makefile.manual --directory=stats $@ $(MAKE) -f Makefile.manual --directory=string $@ $(MAKE) -f Makefile.manual --directory=math $@ From de54f6abe54e53d2133d908eed6ecdb89310fc89 Mon Sep 17 00:00:00 2001 From: jim-215-fisher Date: Wed, 2 Feb 2022 14:01:56 -0500 Subject: [PATCH 05/15] add test file --- .../test_specialfunctions_gamma.fypp | 587 ++++++++++++++++++ 1 file changed, 587 insertions(+) create mode 100644 src/tests/specialfunctions/test_specialfunctions_gamma.fypp diff --git a/src/tests/specialfunctions/test_specialfunctions_gamma.fypp b/src/tests/specialfunctions/test_specialfunctions_gamma.fypp new file mode 100644 index 000000000..38525b58f --- /dev/null +++ b/src/tests/specialfunctions/test_specialfunctions_gamma.fypp @@ -0,0 +1,587 @@ +#:set WITH_QP = False +#:set WITH_XDP = False +#:include "common.fypp" +#:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_KINDS_TYPES +module test_specialfunctions_gamma + use testdrive, only : new_unittest, unittest_type, error_type, check + use stdlib_kinds, only: sp, dp, int8, int16, int32, int64 + use stdlib_specialfunctions_gamma, only: gamma, log_gamma, log_factorial, & + lower_incomplete_gamma, & + upper_incomplete_gamma, & + log_lower_incomplete_gamma, & + log_upper_incomplete_gamma, & + regularized_gamma_p, & + regularized_gamma_q + + implicit none + private + + public :: collect_specialfunctions_gamma + + #:for k1, t1 in REAL_KINDS_TYPES + ${t1}$, parameter :: tol_${k1}$ = 1000 * epsilon(1.0_${k1}$) + #:endfor + + + + +contains + + subroutine collect_specialfunctions_gamma(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + + testsuite = [ & + new_unittest("log_factorial_iint8", test_logfact_iint8) & + + #:for k1, t1 in INT_KINDS_TYPES + , new_unittest("log_factorial_${t1[0]}$${k1}$", & + test_logfact_${t1[0]}$${k1}$) & + #:endfor + + #:for k1, t1 in CI_KINDS_TYPES + , new_unittest("gamma_${t1[0]}$${k1}$", & + test_gamma_${t1[0]}$${k1}$) & + , new_unittest("log_gamma_${t1[0]}$${k1}$", & + test_loggamma_${t1[0]}$${k1}$) & + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + , new_unittest("lower_incomplete_gamma_${t1[0]}$${k1}$${k2}$", & + test_lincgamma_${t1[0]}$${k1}$${k2}$) & + , new_unittest("log_lower_incomplete_gamma_${t1[0]}$${k1}$${k2}$", & + test_log_lincgamma_${t1[0]}$${k1}$${k2}$) & + , new_unittest("upper_incomplete_gamma_${t1[0]}$${k1}$${k2}$", & + test_uincgamma_${t1[0]}$${k1}$${k2}$) & + , new_unittest("log_upper_incomplete_gamma_${t1[0]}$${k1}$${k2}$", & + test_log_uincgamma_${t1[0]}$${k1}$${k2}$) & + , new_unittest("regularized_gamma_p_${t1[0]}$${k1}$${k2}$", & + test_gamma_p_${t1[0]}$${k1}$${k2}$) & + , new_unittest("regularized_gamma_q_${t1[0]}$${k1}$${k2}$", & + test_gamma_q_${t1[0]}$${k1}$${k2}$) & + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + , new_unittest("lower_incomplete_gamma_${t1[0]}$${k1}$", & + test_lincgamma_${t1[0]}$${k1}$) & + , new_unittest("log_lower_incomplete_gamma_${t1[0]}$${k1}$", & + test_log_lincgamma_${t1[0]}$${k1}$) & + , new_unittest("upper_incomplete_gamma_${t1[0]}$${k1}$", & + test_uincgamma_${t1[0]}$${k1}$) & + , new_unittest("log_upper_incomplete_gamma_${t1[0]}$${k1}$", & + test_log_uincgamma_${t1[0]}$${k1}$) & + , new_unittest("regularized_gamma_p_${t1[0]}$${k1}$", & + test_gamma_p_${t1[0]}$${k1}$) & + , new_unittest("regularized_gamma_q_${t1[0]}$${k1}$", & + test_gamma_q_${t1[0]}$${k1}$) & + #:endfor + ] + end subroutine collect_specialfunctions_gamma + + + + #:for k1, t1 in INT_KINDS_TYPES + + subroutine test_logfact_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 6 + integer :: i + + #:if k1 == "int8" + + ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & + 5_${k1}$, 100_${k1}$] + real(sp), parameter :: ans(n) = [0.0, 0.0, 0.693147180, 3.17805383, & + 4.78749174, 3.63739376e2] + + #:elif k1 == "int16" + + ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & + 7_${k1}$, 500_${k1}$] + real(sp), parameter :: ans(n) = [0.0, 0.0, 0.693147180, 3.17805383, & + 8.52516136, 2.61133046e3] + + #:elif k1 == "int32" + + ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & + 12_${k1}$, 7000_${k1}$] + real(sp), parameter :: ans(n) = [0.0, 0.0, 0.693147180, 3.17805383, & + 1.99872145e1, 5.49810038e4] + + #:elif k1 == "int64" + + ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & + 20_${k1}$, 90000_${k1}$] + real(sp), parameter :: ans(n) = [0.0, 0.0, 0.693147180, 3.17805383, & + 4.23356165e1, 9.36687468e5] + + #:endif + + do i = 1, n + + call check(error, log_factorial(x(i)), ans(i), "Integer kind " & + //"${k1}$ failed", thr = tol_sp, rel = .true.) + + end do + end subroutine test_logfact_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in CI_KINDS_TYPES + + subroutine test_gamma_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + + #:if k1 == "int8" + + ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 4_${k1}$, 6_${k1}$] + ${t1}$, parameter :: ans(n) = [1_${k1}$, 1_${k1}$, 6_${k1}$, 120_${k1}$] + + #:elif k1 == "int16" + + ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 4_${k1}$, 8_${k1}$] + ${t1}$, parameter :: ans(n) = [1_${k1}$, 1_${k1}$, 6_${k1}$, 5040_${k1}$] + + #:elif k1 == "int32" + + ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 4_${k1}$, 13_${k1}$] + ${t1}$, parameter :: ans(n) = [1_${k1}$, 1_${k1}$, 6_${k1}$, & + 479001600_${k1}$] + + #:elif k1 == "int64" + + ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 4_${k1}$, 21_${k1}$] + ${t1}$, parameter :: ans(n) = [1_${k1}$, 1_${k1}$, 6_${k1}$, & + 2432902008176640000_${k1}$] + #:elif t1[0] == "c" + + ${t1}$, parameter :: x(n) = [(0.25_${k1}$, 0.25_${k1}$), & + (0.5_${k1}$, -0.5_${k1}$), & + (1.0_${k1}$, 1.0_${k1}$), & + (-1.254e1_${k1}$, -9.87_${k1}$)] + + ${t1}$, parameter :: ans(n) = & + [(1.6511332803889208_${k1}$, -1.8378758749947890_${k1}$), & + (0.81816399954174739_${k1}$, 0.76331382871398262_${k1}$),& + (0.49801566811835604_${k1}$, -0.15494982830181069_${k1}$),& + (-2.18767396709283064e-21_${k1}$, 2.77577940846953455e-21_${k1}$)] + #:endif + + + #:if t1[0] == "i" + + do i = 1, n + + call check(error, gamma(x(i)), ans(i), "Integer kind ${k1}$ failed") + + end do + + #:elif t1[0] == "c" + + do i = 1, n + + call check(error, gamma(x(i)), ans(i), "Complex kind ${k1}$ failed",& + thr = tol_${k1}$, rel = .true.) + + end do + + #:endif + end subroutine test_gamma_${t1[0]}$${k1}$ + + + + subroutine test_loggamma_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + + #:if k1 == "int8" + + ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 10_${k1}$, 47_${k1}$] + real(sp), parameter :: ans(n) = [0.0, 0.0, 1.28018274e1, 1.32952575e2] + + #:elif k1 == "int16" + + ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 111_${k1}$, 541_${k1}$] + real(sp), parameter :: ans(n) = [0.0, 0.0, 4.10322777e2, 2.86151221e3] + + #:elif k1 == "int32" + + ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 2021_${k1}$, & + 42031_${k1}$] + real(sp), parameter :: ans(n) = [0.0, 0.0, 1.33586470e4, 4.05433461e5] + + #:elif k1 == "int64" + + ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 2021_${k1}$, & + 42031_${k1}$] + real(sp), parameter :: ans(n) = [0.0, 0.0, 1.33586470e4, 4.05433461e5] + + #:elif t1[0] == "c" + + ${t1}$, parameter :: x(n) = [(0.25_${k1}$, 0.25_${k1}$), & + (0.5_${k1}$, -0.5_${k1}$), & + (1.0_${k1}$, 1.0_${k1}$), & + (-1.254e1_${k1}$, -9.87_${k1}$)] + + ${t1}$, parameter :: ans(n) = & + [(0.90447450949333889_${k1}$, -0.83887024394321282_${k1}$),& + (0.11238724280962311_${k1}$, 0.75072920212205074_${k1}$), & + (-0.65092319930185634_${k1}$, -0.30164032046753320_${k1}$),& + (-4.7091788015763380e1_${k1}$, 1.4804627819235690e1_${k1}$)] + #:endif + + + #:if t1[0] == "i" + + do i = 1, n + + call check(error, log_gamma(x(i)), ans(i), "Integer kind ${k1}$ " & + //"failed", thr = tol_sp, rel = .true.) + + end do + + #:elif t1[0] == "c" + + do i = 1, n + + call check(error, log_gamma(x(i)), ans(i), "Complex kind ${k1}$ " & + //"failed", thr = tol_${k1}$, rel = .true.) + + end do + + #:endif + end subroutine test_loggamma_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + + subroutine test_lincgamma_${t1[0]}$${k1}$${k2}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1_${k1}$, 2_${k1}$, 3_${k1}$, 2_${k1}$] + ${t2}$ :: x(n) = [0.5_${k2}$, 3.5_${k2}$, -5.0_${k2}$, -10.0_${k2}$] + + ${t2}$, parameter :: ans(n) = [0.3934693402873667_${k2}$, & + 0.86411177459956675_${k2}$, & + -2.5210237047438023e3_${k2}$, & + 1.9823919215326045e5_${k2}$] + + do i = 1, n + + call check(error, lower_incomplete_gamma(p(i), x(i)), ans(i), & + "Lower incomplete gamma function with p(kind=${k1}$) and " & + //"x(kind=${k2}$) failed", thr = tol_${k2}$, rel = .true.) + + end do + + end subroutine test_lincgamma_${t1[0]}$${k1}$${k2}$ + + + + subroutine test_log_lincgamma_${t1[0]}$${k1}$${k2}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1_${k1}$, 2_${k1}$, 3_${k1}$, 2_${k1}$] + ${t2}$ :: x(n) = [0.5_${k2}$, 3.5_${k2}$, -5.0_${k2}$, -10.0_${k2}$] + + ${t2}$, parameter :: ans(n) = [-0.93275212956718857_${k2}$, & + -0.14605314979599791_${k2}$, & + 7.8324203300567640_${k2}$, & + 1.2197229621760137e1_${k2}$] + + do i = 1, n + + call check(error, log_lower_incomplete_gamma(p(i), x(i)), ans(i), & + "Logarithm of lower incomplete gamma function with " & + //"p(kind=${k1}$) and x(kind=${k2}$) failed", thr = tol_${k2}$, & + rel = .true.) + + end do + + end subroutine test_log_lincgamma_${t1[0]}$${k1}$${k2}$ + + + + subroutine test_uincgamma_${t1[0]}$${k1}$${k2}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1_${k1}$, 2_${k1}$, 3_${k1}$, 2_${k1}$] + ${t2}$ :: x(n) = [0.5_${k2}$, 3.5_${k2}$, -5.0_${k2}$, -10.0_${k2}$] + + ${t2}$, parameter :: ans(n) = [0.60653065971263342_${k2}$, & + 0.13588822540043325_${k2}$, & + 2.5230237047438022E3_${k2}$, & + -1.9823819215326045e5_${k2}$] + + do i = 1, n + + call check(error, upper_incomplete_gamma(p(i), x(i)), ans(i), & + "Logarithm of upper incomplete gamma function with " & + //"p(kind=${k1}$) and x(kind=${k2}$) failed", thr = tol_${k2}$, & + rel = .true.) + + end do + + end subroutine test_uincgamma_${t1[0]}$${k1}$${k2}$ + + + + subroutine test_log_uincgamma_${t1[0]}$${k1}$${k2}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1_${k1}$, 2_${k1}$, 3_${k1}$, 2_${k1}$] + ${t2}$ :: x(n) = [0.5_${k2}$, 3.5_${k2}$, -5.0_${k2}$, -10.0_${k2}$] + + ${t2}$, parameter :: ans(n) = [-0.5_${k2}$, -1.9959226032237259_${k2}$,& + 7.8332133440562161_${k2}$, & + 1.2197224577336219e1_${k2}$] + + do i = 1, n + + call check(error, log_upper_incomplete_gamma(p(i), x(i)), ans(i), & + "Logarithm of upper incomplete gamma function with " & + //"p(kind=${k1}$) and x(kind=${k2}$) failed", thr = tol_${k2}$, & + rel = .true.) + + end do + end subroutine test_log_uincgamma_${t1[0]}$${k1}$${k2}$ + + + + + subroutine test_gamma_p_${t1[0]}$${k1}$${k2}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1_${k1}$, 1_${k1}$, 3_${k1}$, 3_${k1}$] + ${t2}$ :: x(n) = [0.5_${k2}$, 1.5_${k2}$, 0.5_${k2}$, 3.5_${k2}$] + + ${t2}$, parameter :: ans(n) = [0.39346934028736658_${k2}$, & + 0.77686983985157017_${k2}$, & + 1.4387677966970687e-2_${k2}$, & + 0.67915280113786593_${k2}$] + + do i = 1, n + + call check(error, regularized_gamma_p(p(i), x(i)), ans(i), & + "Regularized gamma P function with p(kind=${k1}$) and " & + //"x(kind=${k2}$) failed", thr = tol_${k2}$, rel = .true.) + + end do + end subroutine test_gamma_p_${t1[0]}$${k1}$${k2}$ + + + + subroutine test_gamma_q_${t1[0]}$${k1}$${k2}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1_${k1}$, 1_${k1}$, 3_${k1}$, 3_${k1}$] + ${t2}$ :: x(n) = [0.5_${k2}$, 1.5_${k2}$, 0.5_${k2}$, 3.5_${k2}$] + + ${t2}$, parameter :: ans(n) = [0.60653065971263342_${k2}$, & + 0.22313016014842983_${k2}$, & + 0.98561232203302931_${k2}$, & + 0.32084719886213407_${k2}$] + + do i = 1, n + + call check(error, regularized_gamma_q(p(i), x(i)), ans(i), & + "Regularized gamma Q function with p(kind=${k1}$) and " & + //"x(kind=${k2}$) failed", thr = tol_${k2}$, rel = .true.) + + end do + end subroutine test_gamma_q_${t1[0]}$${k1}$${k2}$ + + #:endfor + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + + subroutine test_lincgamma_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1.0_${k1}$, 2.0_${k1}$, 3.1_${k1}$, 6.5_${k1}$] + ${t1}$ :: x(n) = [0.5_${k1}$, 3.5_${k1}$, 5.0_${k1}$, 3.2_${k1}$] + + ${t1}$, parameter :: ans(n) = [0.3934693402873667_${k1}$, & + 0.86411177459956675_${k1}$, & + 1.8980559470963281_${k1}$, & + 2.0043549563092636e1_${k1}$] + + do i = 1, n + + call check(error, lower_incomplete_gamma(p(i), x(i)), ans(i), & + "Lower incomplete gamma function with p(kind=${k1}$) and " & + //"x(kind=${k1}$) failed", thr = tol_${k1}$, rel = .true.) + + end do + + end subroutine test_lincgamma_${t1[0]}$${k1}$ + + + + subroutine test_log_lincgamma_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1.0_${k1}$, 2.0_${k1}$, 3.1_${k1}$, 6.5_${k1}$] + ${t1}$ :: x(n) = [0.5_${k1}$, 3.5_${k1}$, 5.0_${k1}$, 3.2_${k1}$] + + ${t1}$, parameter :: ans(n) = [-0.93275212956718857_${k1}$, & + -0.14605314979599791_${k1}$, & + 0.64083017662175706_${k1}$, & + 2.9979073844388951_${k1}$] + + do i = 1, n + + call check(error, log_lower_incomplete_gamma(p(i), x(i)), ans(i), & + "Logarithm of lower incomplete gamma function with " & + //"p(kind=${k1}$) and x(kind=${k1}$) failed", thr = tol_${k1}$, & + rel = .true.) + + end do + + end subroutine test_log_lincgamma_${t1[0]}$${k1}$ + + + + subroutine test_uincgamma_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1.0_${k1}$, 2.0_${k1}$, 3.1_${k1}$, 6.5_${k1}$] + ${t1}$ :: x(n) = [0.5_${k1}$, 3.5_${k1}$, 5.0_${k1}$, 3.2_${k1}$] + + ${t1}$, parameter :: ans(n) = [0.60653065971263342_${k1}$, & + 0.13588822540043325_${k1}$, & + 0.29956433129614910_${k1}$, & + 2.6784172825195172e2_${k1}$] + + do i = 1, n + + call check(error, upper_incomplete_gamma(p(i), x(i)), ans(i), & + "Logarithm of upper incomplete gamma function with " & + //"p(kind=${k1}$) and x(kind=${k1}$) failed", thr = tol_${k1}$, & + rel = .true.) + + end do + + end subroutine test_uincgamma_${t1[0]}$${k1}$ + + + + subroutine test_log_uincgamma_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1.0_${k1}$, 2.0_${k1}$, 3.1_${k1}$, 6.5_${k1}$] + ${t1}$ :: x(n) = [0.5_${k1}$, 3.5_${k1}$, 5.0_${k1}$, 3.2_${k1}$] + + ${t1}$, parameter :: ans(n) = [-0.5_${k1}$, -1.9959226032237259_${k1}$,& + -1.2054260888453405_${k1}$, & + 5.5903962398338761_${k1}$] + + do i = 1, n + + call check(error, log_upper_incomplete_gamma(p(i), x(i)), ans(i), & + "Logarithm of upper incomplete gamma function with " & + //"p(kind=${k1}$) and x(kind=${k1}$) failed", thr = tol_${k1}$, & + rel = .true.) + + end do + end subroutine test_log_uincgamma_${t1[0]}$${k1}$ + + + + subroutine test_gamma_p_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1.3_${k1}$, 1.3_${k1}$, 3.7_${k1}$, 3.7_${k1}$] + ${t1}$ :: x(n) = [0.5_${k1}$, 2.1_${k1}$, 2.6_${k1}$, 5.1_${k1}$] + + ${t1}$, parameter :: ans(n) = [0.26487356764588505_${k1}$, & + 0.81011791338807457_${k1}$, & + 0.32198359288949589_${k1}$, & + 0.79435732817518852_${k1}$] + + do i = 1, n + + call check(error, regularized_gamma_p(p(i), x(i)), ans(i), & + "Regularized gamma P function with p(kind=${k1}$) and " & + //"x(kind=${k1}$) failed", thr = tol_${k1}$, rel = .true.) + + end do + end subroutine test_gamma_p_${t1[0]}$${k1}$ + + + + subroutine test_gamma_q_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + integer, parameter :: n = 4 + integer :: i + ${t1}$ :: p(n) = [1.3_${k1}$, 1.3_${k1}$, 3.7_${k1}$, 3.7_${k1}$] + ${t1}$ :: x(n) = [0.5_${k1}$, 2.1_${k1}$, 2.6_${k1}$, 5.1_${k1}$] + + ${t1}$, parameter :: ans(n) = [0.73512643235411495_${k1}$, & + 0.18988208661192543_${k1}$, & + 0.67801640711050411_${k1}$, & + 0.20564267182481148_${k1}$] + + do i = 1, n + + call check(error, regularized_gamma_q(p(i), x(i)), ans(i), & + "Regularized gamma Q function with p(kind=${k1}$) and " & + //"x(kind=${k1}$) failed", thr = tol_${k1}$, rel = .true.) + + end do + end subroutine test_gamma_q_${t1[0]}$${k1}$ + + #:endfor +end module test_specialfunctions_gamma + + + +program tester + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_specialfunctions_gamma, only : collect_specialfunctions_gamma + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [new_testsuite("Gamma special function", & + collect_specialfunctions_gamma)] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program tester From e464d4b820d7a24292fcbce342e0caf284634660 Mon Sep 17 00:00:00 2001 From: jim-215-fisher Date: Wed, 2 Feb 2022 14:32:22 -0500 Subject: [PATCH 06/15] change kind qp --- src/stdlib_specialfunctions_gamma.fypp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp index a8542c43c..7129fddf4 100644 --- a/src/stdlib_specialfunctions_gamma.fypp +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -3,7 +3,8 @@ #:include "common.fypp" #:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_KINDS_TYPES module stdlib_specialfunctions_gamma - use stdlib_kinds, only : sp, dp, xdp, qp, int8, int16, int32, int64 + 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 From d3ab59bc4fb84dc72069a7c04c4e901dc161b8ce Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Wed, 2 Feb 2022 15:01:07 -0500 Subject: [PATCH 07/15] change to tab line 26 --- src/tests/Makefile.manual | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests/Makefile.manual b/src/tests/Makefile.manual index 90bedf3b3..5f3722d47 100644 --- a/src/tests/Makefile.manual +++ b/src/tests/Makefile.manual @@ -23,7 +23,7 @@ all test clean:: $(MAKE) -f Makefile.manual --directory=selection $@ $(MAKE) -f Makefile.manual --directory=sorting $@ $(MAKE) -f Makefile.manual --directory=quadrature $@ - $(MAKE) -f Makefile.manual --directory=specialfunctions $@ + $(MAKE) -f Makefile.manual --directory=specialfunctions $@ $(MAKE) -f Makefile.manual --directory=stats $@ $(MAKE) -f Makefile.manual --directory=string $@ $(MAKE) -f Makefile.manual --directory=math $@ From 0d3505c4be99c1d1b91b00ca907415ff604a63d5 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 22 Feb 2022 18:53:59 -0500 Subject: [PATCH 08/15] Update doc/specs/stdlib_specialfunctions_gamma.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_specialfunctions_gamma.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_specialfunctions_gamma.md b/doc/specs/stdlib_specialfunctions_gamma.md index da1e17436..565789121 100644 --- a/doc/specs/stdlib_specialfunctions_gamma.md +++ b/doc/specs/stdlib_specialfunctions_gamma.md @@ -6,7 +6,7 @@ title: specialfunctions_gamma [TOC] -## `gamma` - Calculate gamma function with any type of argument +## `gamma` - Calculate the gamma function ### Status From 33ae27e746b327e60e56a07b755e4439ba807ab9 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 22 Feb 2022 18:54:30 -0500 Subject: [PATCH 09/15] Update doc/specs/stdlib_specialfunctions_gamma.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_specialfunctions_gamma.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_specialfunctions_gamma.md b/doc/specs/stdlib_specialfunctions_gamma.md index 565789121..71bcc0582 100644 --- a/doc/specs/stdlib_specialfunctions_gamma.md +++ b/doc/specs/stdlib_specialfunctions_gamma.md @@ -18,7 +18,7 @@ The gamma function is defined as the analytic continuation of a convergent impro \Gamma(z)=\int_{0}^{\infty}x^{z-1}e^{-x}dx, \;\; z\in \mathbb{C} \setminus 0, -1, -2, \cdots -Fortran 2018 standard implements intrinsic gamma function of real type argument in single and double precisions. Here the gamma function is extended to both integer and complex arguments. The values of gamma function with integer arguments are exact. The values of gamma function with complex arguments are approximated in single and double precisions by using Lanczos approxiamtion. +Fortran 2018 standard implements the intrinsic gamma function of real type argument in single and double precisions. Here the gamma function is extended to both integer and complex arguments. The values of the gamma function with integer arguments are exact. The values of the gamma function with complex arguments are approximated in single and double precisions by using Lanczos approximation. ### Syntax From 6d94b71cfebb9577d0fb3274d4846b029e439bd4 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 22 Feb 2022 18:54:43 -0500 Subject: [PATCH 10/15] Update doc/specs/stdlib_specialfunctions_gamma.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_specialfunctions_gamma.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_specialfunctions_gamma.md b/doc/specs/stdlib_specialfunctions_gamma.md index 71bcc0582..8ad005bff 100644 --- a/doc/specs/stdlib_specialfunctions_gamma.md +++ b/doc/specs/stdlib_specialfunctions_gamma.md @@ -77,7 +77,7 @@ program demo_gamma end program demo_gamma ``` -## `log_gamma` - calculate natural logarithm of gamma function with any type of argument +## `log_gamma` - Calculate the natural logarithm of the gamma function ### Status From c434857a15714c0ee7cf7ad392de12b4c572fae6 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 22 Feb 2022 19:31:03 -0500 Subject: [PATCH 11/15] Update doc/specs/stdlib_specialfunctions_gamma.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_specialfunctions_gamma.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_specialfunctions_gamma.md b/doc/specs/stdlib_specialfunctions_gamma.md index 8ad005bff..34b9c68d5 100644 --- a/doc/specs/stdlib_specialfunctions_gamma.md +++ b/doc/specs/stdlib_specialfunctions_gamma.md @@ -145,7 +145,7 @@ program demo_log_gamma end program demo_log_gamma ``` -## `log_factorial` - calculate logarithm of a factorial +## `log_factorial` - calculate the logarithm of a factorial ### Status From 88c00fd23ab6976f56cfdea616f83908b89289d7 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 22 Feb 2022 19:31:31 -0500 Subject: [PATCH 12/15] Update doc/specs/stdlib_specialfunctions_gamma.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_specialfunctions_gamma.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_specialfunctions_gamma.md b/doc/specs/stdlib_specialfunctions_gamma.md index 34b9c68d5..25df93e82 100644 --- a/doc/specs/stdlib_specialfunctions_gamma.md +++ b/doc/specs/stdlib_specialfunctions_gamma.md @@ -242,7 +242,7 @@ program demo_ligamma end demo_ligamma ``` -## `upper_incomplete_gamma` - calculate upper incomplete gamma integral +## `upper_incomplete_gamma` - calculate the upper incomplete gamma integral ### Status From 64a80eff83816c7a518d9bc0fbb655ec72e269ca Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 22 Feb 2022 19:34:12 -0500 Subject: [PATCH 13/15] Update src/tests/specialfunctions/test_specialfunctions_gamma.fypp Co-authored-by: Jeremie Vandenplas --- src/tests/specialfunctions/test_specialfunctions_gamma.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tests/specialfunctions/test_specialfunctions_gamma.fypp b/src/tests/specialfunctions/test_specialfunctions_gamma.fypp index 38525b58f..26421bded 100644 --- a/src/tests/specialfunctions/test_specialfunctions_gamma.fypp +++ b/src/tests/specialfunctions/test_specialfunctions_gamma.fypp @@ -270,8 +270,8 @@ contains type(error_type), allocatable, intent(out) :: error integer, parameter :: n = 4 integer :: i - ${t1}$ :: p(n) = [1_${k1}$, 2_${k1}$, 3_${k1}$, 2_${k1}$] - ${t2}$ :: x(n) = [0.5_${k2}$, 3.5_${k2}$, -5.0_${k2}$, -10.0_${k2}$] + ${t1}$, parameter :: p(n) = [1_${k1}$, 2_${k1}$, 3_${k1}$, 2_${k1}$] + ${t2}$, parameter :: x(n) = [0.5_${k2}$, 3.5_${k2}$, -5.0_${k2}$, -10.0_${k2}$] ${t2}$, parameter :: ans(n) = [0.3934693402873667_${k2}$, & 0.86411177459956675_${k2}$, & From 039dfabbc8fbc7c27706e28c42cb1281bb2c61b0 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 22 Feb 2022 19:38:22 -0500 Subject: [PATCH 14/15] Update doc/specs/stdlib_specialfunctions_gamma.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_specialfunctions_gamma.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_specialfunctions_gamma.md b/doc/specs/stdlib_specialfunctions_gamma.md index 25df93e82..a34dc452a 100644 --- a/doc/specs/stdlib_specialfunctions_gamma.md +++ b/doc/specs/stdlib_specialfunctions_gamma.md @@ -290,7 +290,7 @@ program demo_uigamma end program demo_uigamma ``` -## `log_lower_incomplete_gamma` - calculate logarithm of the lower incomplete gamma integral +## `log_lower_incomplete_gamma` - calculate the natural logarithm of the lower incomplete gamma integral ### Status From d11e0dd0407941da9cf3947a4122ef560d4c8341 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 22 Feb 2022 19:55:51 -0500 Subject: [PATCH 15/15] wording change --- doc/specs/stdlib_specialfunctions_gamma.md | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/doc/specs/stdlib_specialfunctions_gamma.md b/doc/specs/stdlib_specialfunctions_gamma.md index a34dc452a..bd277d4d8 100644 --- a/doc/specs/stdlib_specialfunctions_gamma.md +++ b/doc/specs/stdlib_specialfunctions_gamma.md @@ -39,23 +39,23 @@ The function returns a value with the same type and kind as input argument. ### Example ```fortran program demo_gamma - use iso_fortran_env, only : real64, int64 + use stdlib_kinds, only : dp, int64 use stdlib_specialfunctions_gamma, only : gamma implicit none integer :: i integer(int64) :: n real :: x - real(real64) :: y + real(dp) :: y complex :: z - complex(real64) :: z1 + complex(dp) :: z1 i = 10 n = 15_int64 x = 2.5 - y = 4.3_real64 + y = 4.3_dp z = (2.3, 0.6) - z1 = (-4.2_real64, 3.1_real64) + z1 = (-4.2_dp, 3.1_dp) print *, gamma(i) !integer gives exact result ! 362880 @@ -85,7 +85,7 @@ Experimental ### Description -Due to the different branch cut structures and a different principal branch, natural logarithm of gamma function log_gamma(z) with complex argument is different from the ln(Gamma(z)). The two have the same real part but different imaginary part. +Mathematically, logarithm of gamma function is a special function with complex arguments by itself. Due to the different branch cut structures and a different principal branch, natural logarithm of gamma function log_gamma(z) with complex argument is different from the ln(Gamma(z)). The two have the same real part but different imaginary part. Fortran 2018 standard implements intrinsic log_gamma function of absolute value of real type argument in single and double precision. Here the log_gamma function is extended to both integer and complex arguments. The values of log_gamma function with complex arguments are approximated in single and double precisions by using Stirling's approximation. @@ -109,15 +109,15 @@ The function returns real single precision values for integer input arguments, w ```fortran program demo_log_gamma - use iso_fortran_env, only : real64 + use stdlib_kinds, only : dp use stdlib_specialfunctions_gamma, only : log_gamma implicit none integer :: i real :: x - real(real64) :: y + real(dp) :: y complex :: z - complex(real64) :: z1 + complex(dp) :: z1 i = 10 x = 8.76 @@ -174,7 +174,7 @@ The function returns real type values with single precision. ### Example ```fortran program demo_log_factorial - use iso_fortran_env, only : int64 + use stdlib_kinds, only : int64 use stdlib_specialfunctions_gamma, only : lf => log_factorial implicit none integer :: n