diff --git a/doc/specs/stdlib_specialfunctions_gamma.md b/doc/specs/stdlib_specialfunctions_gamma.md new file mode 100644 index 000000000..bd277d4d8 --- /dev/null +++ b/doc/specs/stdlib_specialfunctions_gamma.md @@ -0,0 +1,437 @@ +--- +title: specialfunctions_gamma +--- + +# Special functions gamma + +[TOC] + +## `gamma` - Calculate the gamma function + +### Status + +Experimental + +### Description + +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 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 + +`result = [[stdlib_specialfunctions_gamma(module):gamma(interface)]] (x)` + +### Class + +Elemental function + +### Arguments + +`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 stdlib_kinds, only : dp, int64 + use stdlib_specialfunctions_gamma, only : gamma + implicit none + + integer :: i + integer(int64) :: n + real :: x + real(dp) :: y + complex :: z + complex(dp) :: z1 + + i = 10 + n = 15_int64 + x = 2.5 + y = 4.3_dp + z = (2.3, 0.6) + z1 = (-4.2_dp, 3.1_dp) + + print *, gamma(i) !integer gives exact result +! 362880 + + print *, gamma(n) +! 87178291200 + + print *, gamma(x) ! intrinsic function call +! 1.32934034 + + 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 the natural logarithm of the gamma function + +### Status + +Experimental + +### Description + +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. + +### Syntax + +`result = [[stdlib_specialfunctions_gamma(module):log_gamma(interface)]] (x)` + +### Class + +Elemental function + +### Arguments + +`x`: Shall be a positive integer or a complex type number. + +### Return value + +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_kinds, only : dp + use stdlib_specialfunctions_gamma, only : log_gamma + implicit none + + integer :: i + real :: x + real(dp) :: y + complex :: z + complex(dp) :: z1 + + i = 10 + x = 8.76 + y = x + z = (5.345, -3.467) + z1 = z + print *, log_gamma(i) !default single precision output +!12.8018274 + + print *, log_gamma(x) !intrinsic function call + +!10.0942659 + + print *, log_gamma(y) !intrinsic function call + +!10.094265528673880 + + print *, log_gamma(z) !same kind as input + +!(2.56165648, -5.73382425) + + print *, log_gamma(z1) + +!(2.5616575105114614, -5.7338247782852498) +end program demo_log_gamma +``` + +## `log_factorial` - calculate the logarithm of a factorial + +### Status + +Experimental + +### Description + +Compute the natural logarithm of factorial, log(n!) + +### Syntax + +`result = [[stdlib_specialfunctions_gamma(module):log_factorial(interface)]] (x)` + +### Class + +Elemental function + +### Arguments + +`x`: Shall be a positive integer type number. + +### Return value + +The function returns real type values with single precision. + +### Example +```fortran +program demo_log_factorial + use stdlib_kinds, 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 + +### 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\in \mathbb{R} + +When x < 0, p must be positive integer. + +### 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. + +### 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 the 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 \in \mathbb{R} + +When x < 0, p must be a positive integer. + +### 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. + +### 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 the natural logarithm of the lower incomplete gamma integral + +### Status + +Experimental + +### Description + +Compute the natural 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 natural 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. + +### 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 + +### 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. + +### 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/CMakeLists.txt b/src/CMakeLists.txt index dce8e3054..40a5940a6 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/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp new file mode 100644 index 000000000..7129fddf4 --- /dev/null +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -0,0 +1,1289 @@ +#:set WITH_QP = False +#:set WITH_XDP = False +#:include "common.fypp" +#:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_KINDS_TYPES +module stdlib_specialfunctions_gamma + use iso_fortran_env, only : qp => real128 + use stdlib_kinds, only : sp, dp, int8, int16, int32, int64 + use stdlib_error, only : error_stop + + implicit none + private + + integer(int8), parameter :: max_fact_int8 = 6_int8 + integer(int16), parameter :: max_fact_int16 = 8_int16 + integer(int32), parameter :: max_fact_int32 = 13_int32 + integer(int64), parameter :: max_fact_int64 = 21_int64 + + #:for k1, t1 in REAL_KINDS_TYPES + ${t1}$, parameter :: tol_${k1}$ = epsilon(1.0_${k1}$) + #:endfor + real(qp), parameter :: tol_qp = epsilon(1.0_qp) + + + + public :: gamma, log_gamma, log_factorial + public :: lower_incomplete_gamma, log_lower_incomplete_gamma + public :: upper_incomplete_gamma, log_upper_incomplete_gamma + public :: regularized_gamma_p, regularized_gamma_q + + + + interface gamma + !! Gamma function for integer and complex numbers + !! + #:for k1, t1 in CI_KINDS_TYPES + module procedure gamma_${t1[0]}$${k1}$ + #:endfor + end interface gamma + + + + interface log_gamma + !! Logarithm of gamma function + !! + #:for k1, t1 in CI_KINDS_TYPES + module procedure l_gamma_${t1[0]}$${k1}$ + #:endfor + end interface log_gamma + + + + interface log_factorial + !! Logarithm of factorial n!, integer variable + !! + #:for k1, t1 in INT_KINDS_TYPES + module procedure l_factorial_${t1[0]}$${k1}$ + #:endfor + end interface log_factorial + + + + interface lower_incomplete_gamma + !! Lower incomplete gamma function + !! + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure ingamma_low_${t1[0]}$${k1}$${k2}$ + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + module procedure ingamma_low_${t1[0]}$${k1}$ + #:endfor + end interface lower_incomplete_gamma + + + + interface log_lower_incomplete_gamma + !! Logarithm of lower incomplete gamma function + !! + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure l_ingamma_low_${t1[0]}$${k1}$${k2}$ + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + module procedure l_ingamma_low_${t1[0]}$${k1}$ + #:endfor + end interface log_lower_incomplete_gamma + + + + interface upper_incomplete_gamma + !! Upper incomplete gamma function + !! + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure ingamma_up_${t1[0]}$${k1}$${k2}$ + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + module procedure ingamma_up_${t1[0]}$${k1}$ + #:endfor + end interface upper_incomplete_gamma + + + + interface log_upper_incomplete_gamma + !! Logarithm of upper incomplete gamma function + !! + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure l_ingamma_up_${t1[0]}$${k1}$${k2}$ + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + module procedure l_ingamma_up_${t1[0]}$${k1}$ + #:endfor + end interface log_upper_incomplete_gamma + + + + interface regularized_gamma_p + !! Regularized (normalized) lower incomplete gamma function, P + !! + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure regamma_p_${t1[0]}$${k1}$${k2}$ + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + module procedure regamma_p_${t1[0]}$${k1}$ + #:endfor + end interface regularized_gamma_p + + + + interface regularized_gamma_q + !! Regularized (normalized) upper incomplete gamma function, Q + !! + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure regamma_q_${t1[0]}$${k1}$${k2}$ + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + module procedure regamma_q_${t1[0]}$${k1}$ + #:endfor + end interface regularized_gamma_q + + + + interface gpx + ! Incomplete gamma G function. + ! Internal use only + ! + #:for k1, t1 in REAL_KINDS_TYPES + module procedure gpx_${t1[0]}$${k1}$ !for real p and x + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure gpx_${t1[0]}$${k1}$${k2}$ !for integer p and real x + #:endfor + #:endfor + end interface gpx + + + + interface l_gamma + ! Logarithm of gamma with integer argument for designated output kind. + ! Internal use only + ! + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + module procedure l_gamma_${t1[0]}$${k1}$${k2}$ + #:endfor + #:endfor + end interface l_gamma + + + + + +contains + + #:for k1, t1 in INT_KINDS_TYPES + impure elemental function gamma_${t1[0]}$${k1}$(z) result(res) + ${t1}$, intent(in) :: z + ${t1}$ :: res, i + ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$ + + if(z <= zero) call error_stop("Error(gamma): Gamma function argument" & + //" must be positive integer.") + + if(z > max_fact_${k1}$) call error_stop("Error(gamma): Gamma function" & + //" integer argument is greater than the upper limit from which an"& + //" integer overflow will be generated. Suggest switch to high " & + //" precision or convert to real data type") + + res = one + + do i = one, z - one + + res = res * i + + end do + + end function gamma_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in CMPLX_KINDS_TYPES + #:if k1 == "sp" + #:set k2 = "dp" + #:elif k1 == "dp" + #:set k2 = "qp" + #:endif + #:set t2 = "real({})".format(k2) + + impure elemental function gamma_${t1[0]}$${k1}$(z) result(res) + ${t1}$, intent(in) :: z + ${t1}$ :: res + integer :: i + + real(${k1}$), parameter :: zero_k1 = 0.0_${k1}$ + ${t2}$, parameter :: zero = 0.0_${k2}$, half = 0.5_${k2}$, & + one = 1.0_${k2}$, pi = acos(- one), sqpi = sqrt(pi) + complex(${k2}$) :: y, x, sum + + #:if k1 == "sp" + #! for single precision input, using double precision for calculation + + integer, parameter :: n = 10 + ${t2}$, parameter :: r = 10.900511_${k2}$ + ${t2}$, parameter :: d(0 : n) = [2.48574089138753566e-5_${k2}$, & + 1.05142378581721974_${k2}$, & + -3.45687097222016235_${k2}$, & + 4.51227709466894824_${k2}$, & + -2.98285225323576656_${k2}$, & + 1.05639711577126713_${k2}$, & + -1.95428773191645870e-1_${k2}$, & + 1.70970543404441224e-2_${k2}$, & + -5.71926117404305781e-4_${k2}$, & + 4.63399473359905637e-6_${k2}$, & + -2.71994908488607704e-9_${k2}$] + ! parameters from above referenced source. + + #:elif k1 == "dp" + #! for double precision input, using quadruple precision for calculation + + integer, parameter :: n = 24 + ${t2}$, parameter :: r = 25.617904_${k2}$ + ${t2}$, parameter :: d(0 : n)= & + [1.0087261714899910504854136977047144166e-11_${k2}$, & + 1.6339627701280724777912729825256860624_${k2}$, & + -1.4205787702221583745972794018472259342e+1_${k2}$, & + 5.6689501646428786119793943350900908698e+1_${k2}$, & + -1.3766376824252176069406853670529834070e+2_${k2}$, & + 2.2739972766608392140035874845640820558e+2_${k2}$, & + -2.7058382145757164380300118233258834430e+2_${k2}$, & + 2.39614374587263042692333711131832094166e+2_${k2}$, & + -1.6090450559507517723393498276315290189e+2_${k2}$, & + 8.27378183187161305711485619113605553100e+1_${k2}$, & + -3.2678977082742592701862249152153110206e+1_${k2}$, & + 9.89018079175824824537131521501652931756_${k2}$, & + -2.2762136356329318377213053650799013041_${k2}$, & + 3.93265017303573867227590563182750070164e-1_${k2}$, & + -5.0051054352146209116457193223422284239e-2_${k2}$, & + 4.57142601898244576789629257292603538238e-3_${k2}$, & + -2.8922592124650765614787233510990416584e-4_${k2}$, & + 1.20833375377219592849746118012697473202e-5_${k2}$, & + -3.1220812187551248389268359432609135033e-7_${k2}$, & + 4.55117045361638520378367871355819524460e-9_${k2}$, & + -3.2757632817493581828033170342853173968e-11_${k2}$, & + 9.49784279240135747819870224486376897253e-14_${k2}$, & + -7.9480594917454410117072562195702526836e-17_${k2}$, & + 1.04692819439870077791406760109955648941e-20_${k2}$, & + -5.8990280044857540075384586350723191533e-26_${k2}$] + ! parameters from above referenced source. + + #:endif + + + + if(abs(z % im) < tol_${k1}$) then + + res = cmplx(gamma(z % re), kind = ${k1}$) + return + + end if + + if(z % re > zero_k1) then + + 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 CMPLX_KINDS_TYPES + #:if k1 == "sp" + #:set k2 = "dp" + #:elif k1 == "dp" + #:set k2 = "qp" + #:endif + #:set t2 = "real({})".format(k2) + impure elemental function l_gamma_${t1[0]}$${k1}$(z) result (res) + ! + ! log_gamma function for any complex number, excluding negative whole number + ! "Computation of special functions", Shanjie Zhang & Jianmin Jin, 1996, p.48 + ! "Computing the principal branch of log-gamma", D.E.G. Hare, + ! J. of Algorithms, 25(2), 1997 p. 221–236 + ! + ! Fortran 90 program by Jim-215-Fisher + ! + ${t1}$, intent(in) :: z + ${t1}$ :: res, z1, z2 + real(${k1}$) :: d + integer :: m, i + complex(${k2}$) :: zr, zr2, sum, s + real(${k1}$), parameter :: z_limit = 10_${k1}$, zero_k1 = 0.0_${k1}$ + integer, parameter :: n = 20 + ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$, & + pi = acos(-one), ln2pi = log(2 * pi) + ${t2}$, parameter :: a(n) = [ & + .8333333333333333333333333333333333333333E-1_${k2}$,& + -.2777777777777777777777777777777777777778E-2_${k2}$,& + .7936507936507936507936507936507936507937E-3_${k2}$,& + -.5952380952380952380952380952380952380952E-3_${k2}$,& + .8417508417508417508417508417508417508418E-3_${k2}$,& + -.1917526917526917526917526917526917526918E-2_${k2}$,& + .6410256410256410256410256410256410256410E-2_${k2}$,& + -.2955065359477124183006535947712418300654E-1_${k2}$,& + .1796443723688305731649384900158893966944E+0_${k2}$,& + -.1392432216905901116427432216905901116427E+1_${k2}$,& + .1340286404416839199447895100069013112491E+2_${k2}$,& + -.1568482846260020173063651324520889738281E+3_${k2}$,& + .2193103333333333333333333333333333333333E+4_${k2}$,& + -.3610877125372498935717326521924223073648E+5_${k2}$,& + .6914722688513130671083952507756734675533E+6_${k2}$,& + -.1523822153940741619228336495888678051866E+8_${k2}$,& + .3829007513914141414141414141414141414141E+9_${k2}$,& + -.1088226603578439108901514916552510537473E+11_${k2}$,& + .3473202837650022522522522522522522522523E+12_${k2}$,& + -.1236960214226927445425171034927132488108E+14_${k2}$] + ! parameters from above reference + + z2 = z + + if(z % re < zero_k1) then + + z2 = cmplx(abs(z % re), - z % im, kind = ${k1}$) + 1 + + end if + + d = hypot(z2 % re, z2 % im) + z1 = z2 + m = 0 + + if(d <= z_limit) then !for small |z| + + m = ceiling(z_limit - d) + z1 = z2 + m + + end if + + zr = one / z1 + zr2 = zr * zr + + sum = (((a(20) * zr2 + a(19)) * zr2 + a(18)) * zr2 + a(17)) * zr2 + sum = (((sum + a(16)) * zr2 + a(15)) * zr2 + a(14)) * zr2 + sum = (((sum + a(13)) * zr2 + a(12)) * zr2 + a(11)) * zr2 + sum = (((sum + a(10)) * zr2 + a(9)) * zr2 + a(8)) * zr2 + sum = (((sum + a(7)) * zr2 + a(6)) * zr2 + a(5)) * zr2 + sum = (((sum + a(4)) * zr2 + a(3)) * zr2 + a(2)) * zr2 + sum = (sum + a(1)) * zr + ln2pi / 2 - z1 + (z1 - 0.5_${k2}$) * log(z1) + + if(m /= 0) then + + s = cmplx(zero, zero, kind = ${k2}$) + + do i = 1, m + + s = s + log(cmplx(z1, kind = ${k2}$) - i) + + end do + + sum = sum - s + + end if + + if(z % re < zero_k1) then + + sum = log(pi) - log(sin(pi * z)) - sum + m = ceiling((2 * z % re - 3) / 4) + sum % im = sum % im + 2 * pi * m * sign(1.0_${k1}$, z % im) + + end if + + res = cmplx(sum, kind = ${k1}$) + end function l_gamma_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in INT_KINDS_TYPES + impure elemental function l_factorial_${t1[0]}$${k1}$(n) result(res) + ! + ! Log(n!) + ! + ${t1}$, intent(in) :: n + real :: res + ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$ + real, parameter :: zero_k2 = 0.0 + + if(n < zero) call error_stop("Error(l_factorial): Logarithm of" & + //" factorial function argument must be non-negative") + + select case(n) + + case (zero) + + res = zero_k2 + + case (one) + + res = zero_k2 + + case (two : ) + + res = l_gamma(n + 1, 1.0D0) + + end select + end function l_factorial_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + #:if k1 == "sp" + #:set k2 = "dp" + #:elif k1 == "dp" + #:set k2 = "qp" + #:endif + #:set t2 = "real({})".format(k2) + + impure elemental function gpx_${t1[0]}$${k1}$(p, x) result(res) + ! + ! Approximation of incomplete gamma G function with real argument p. + ! + ! Based on Rémy Abergel and Lionel Moisan "Algorithm 1006, Fast and + ! Accurate Evaluation of a Generalized Incomplete Gamma Function", ACM + ! Transactions on Mathematical Software, March 2020. + ! + ! Fortran 90 program by Jim-215-Fisher + ! + ${t1}$, intent(in) :: p, x + integer :: n, m + + ${t2}$ :: res, p_lim, a, b, g, c, d, y, ss + ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$ + ${t2}$, parameter :: dm = tiny(1.0_${k2}$) * 10 ** 6 + ${t1}$, parameter :: zero_k1 = 0.0_${k1}$ + + if(p <= zero_k1) call error_stop("Error(gpx): Incomplete gamma" & + //" function must have a positive parameter p") + + if(x < -9.0_${k1}$) then + + p_lim = 5.0_${k1}$ * (sqrt(abs(x)) - 1.0_${k1}$) + + elseif(x >= -9.0_${k1}$ .and. x <= zero_k1) then + + p_lim = zero_k1 + + else + + p_lim = x + + endif + + if(x < zero_k1 .and. p < p_lim .and. abs(anint(p) - p) > tol_${k1}$) & + call error_stop("Error(gpx): Incomplete gamma function with " & + //"negative x must come with a whole number p not too small") + + if(p >= p_lim) then !use modified Lentz method of continued fraction + !for eq. (15) in the above reference. + a = one + b = p + g = a / b + c = a / dm + d = one / b + n = 2 + + do + + if(mod(n, 2) == 0) then + a = (one - p - n / 2) * x + else + a = (n / 2) * x + end if + + b = p - one + n + d = d * a + b + + if(d == zero) d = dm + + c = b + a / c + + if(c == zero) c = dm + + d = one / d + y = c * d + g = g * y + n = n + 1 + + if(abs(y - one) < tol_${k2}$) exit + + end do + + else if(x >= zero_k1) then !use modified Lentz method of continued + !fraction for eq. (16) in the reference. + a = one + b = x + one - p + g = a / b + c = a / dm + d = one / b + n = 2 + + do + + a = (n - 1) * (1 + p - n) + b = b + 2 + d = d * a + b + + if(d == zero) d = dm + + c = b + a / c + + if(c == zero) c = dm + + d = one / d + y = c * d + g = g * y + n = n + 1 + + if(abs(y - one) < tol_${k2}$) exit + + end do + + else !Algorithm 2 in the reference + + m = nint(ss) + a = - x + c = one / a + d = p - one + b = c * (a - d) + n = 1 + + do + + c = d * (d - one) / (a * a) + d = d - 2 + y = c * (a - d) + b = b + y + n = n + 1 + + if(n > int((p - 2) / 2) .or. y < b * tol_${k2}$) exit + + end do + + if(y >= b * tol_${k2}$ .and. mod(m , 2) /= 0) b = b + d * c / a + + g = ((-1) ** m * exp(-a + log_gamma(p) - (p - 1) * log(a)) + b) / a + end if + + res = g + end function gpx_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + impure elemental function gpx_${t1[0]}$${k1}$${k2}$(p, x) result(res) + ! + ! Approximation of incomplete gamma G function with integer argument p. + ! + ! Based on Rémy Abergel and Lionel Moisan "Algorithm 1006, Fast and + ! Accurate Evaluation of a Generalized Incomplete Gamma Function", ACM + ! Transactions on Mathematical Software, March 2020. + ! + ${t1}$, intent(in) :: p + ${t2}$, intent(in) :: x + ${t2}$ :: res, p_lim, a, b, g, c, d, y + integer :: n, m + ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$ + ${t2}$, parameter :: dm = tiny(1.0_${k2}$) * 10 ** 6 + ${t1}$, parameter :: zero_k1 = 0_${k1}$, two = 2_${k1}$ + + if(p <= zero_k1) call error_stop("Error(gpx): Incomplete gamma " & + //"function must have a positive parameter p") + + if(x < -9.0_${k2}$) then + + p_lim = 5.0_${k2}$ * (sqrt(abs(x)) - 1.0_${k2}$) + + else if(x >= -9.0_${k2}$ .and. x <= zero) then + + p_lim = zero + + else + + p_lim = x + + end if + + if(real(p, ${k2}$) >= p_lim) then + + a = one + b = p + g = a / b + c = a / dm + d = one / b + n = 2 + + do + + if(mod(n, 2) == 0) then + + a = (1 - p - n / 2) * x + + else + + a = (n / 2) * x + + end if + + b = p - 1 + n + d = d * a + b + + if(d == zero) d = dm + + c = b + a / c + + if(c == zero) c = dm + + d = one / d + y = c * d + g = g * y + n = n + 1 + + if(abs(y - one) < tol_${k2}$) exit + + end do + + else if(x >= zero) then + + a = one + b = x + 1 - p + g = a / b + c = a / dm + d = one / b + n = 2 + + do + + a = -(n - 1) * (n - 1 - p) + b = b + 2 + d = d * a + b + + if(d == zero) d = dm + + c = b + a / c + + if(c == zero) c = dm + + d = one / d + y = c * d + g = g * y + n = n + 1 + + if(abs(y - one) < tol_${k2}$) exit + + end do + + else + + a = -x + c = one / a + d = p - 1 + b = c * (a - d) + n = 1 + + do + + c = d * (d - one) / (a * a) + d = d - 2 + y = c * ( a - d) + b = b + y + n = n + 1 + + if(int(n, ${k1}$) > (p - two) / two .or. y < b * tol_${k2}$) exit + + end do + + if(y >= b * tol_${k2}$ .and. mod(p, two) /= zero_k1) & + b = b + d * c / a + + g = ((-1) ** p * exp(-a + l_gamma(p, one) - (p - 1) * log(a)) & + + b ) / a + + end if + + res = g + end function gpx_${t1[0]}$${k1}$${k2}$ + + #:endfor + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function ingamma_low_${t1[0]}$${k1}$(p, x) result(res) + ! + ! Approximation of lower incomplete gamma function with real p. + ! + ${t1}$, intent(in) :: p, x + ${t1}$ :: res, s1, y + ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ + + if(x == zero) then + + res = zero + + else if(x > p) then + + s1 = log_gamma(p) + y = one - exp(-x + p * log(x) - s1) * gpx(p, x) + res = exp(s1 + log(y)) + + else if(x <= p .and. x > zero) then + + s1 = -x + p * log(x) + res = gpx(p, x) * exp(s1) + + else + + call error_stop("Error(Logarithm of upper incomplete gamma " & + //"function): negative x must be with integer p") + + end if + end function ingamma_low_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + impure elemental function ingamma_low_${t1[0]}$${k1}$${k2}$(p, x) & + result(res) + ! + ! Approximation of lower incomplete gamma function with integer p. + ! + ${t1}$, intent(in) :: p + ${t2}$, intent(in) :: x + ${t2}$ :: res, s1, y + ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$ + + if(x == zero) then + + res = zero + + else if(x > real(p, ${k2}$)) then + + s1 = l_gamma(p, one) + y = one - exp(-x + p * log(x) - s1) * gpx(p, x) + res = exp(s1 + log(y)) + + else if(x <= real(p, ${k2}$) .and. x > zero) then + + s1 = -x + p * log(x) + res = gpx(p, x) * exp(s1) + + else + + s1 = -x + p * log(abs(x)) + res = gpx(p, x) * exp(s1) + res = (-1) ** p * res + + end if + end function ingamma_low_${t1[0]}$${k1}$${k2}$ + + #:endfor + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function l_ingamma_low_${t1[0]}$${k1}$(p, x) result(res) + + ${t1}$, intent(in) :: p, x + ${t1}$ :: res, s1, y + ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ + + + if(x == zero) then + + res = zero + + else if(x > p) then + + s1 = log_gamma(p) + y = one - exp(-x + p * log(x) - s1) * gpx(p, x) + res = s1 + log(y) + + else if(x <= p .and. x > zero) then + + s1 = -x + p * log(abs(x)) + res = log(abs(gpx(p, x))) + s1 + + else + + call error_stop("Error(Logarithm of upper incomplete gamma " & + //"function): negative x must be with integer p") + + end if + end function l_ingamma_low_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + impure elemental function l_ingamma_low_${t1[0]}$${k1}$${k2}$(p, x) & + result(res) + + ${t1}$, intent(in) :: p + ${t2}$, intent(in) :: x + ${t2}$ :: res, s1, y + ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$ + + if(x == zero) then + + res = zero + + else if(x > real(p, ${k2}$)) then + + s1 = l_gamma(p, one) + y = one - exp(-x + p * log(x) - s1) * gpx(p, x) + res = s1 + log(y) + + else if(x <= real(p, ${k2}$)) then + + s1 = -x + p * log(abs(x)) + res = log(abs(gpx(p, x))) + s1 + + end if + end function l_ingamma_low_${t1[0]}$${k1}$${k2}$ + + #:endfor + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function ingamma_up_${t1[0]}$${k1}$(p, x) result(res) + ! + ! Approximation of upper incomplete gamma function with real p. + ! + ${t1}$, intent(in) :: p, x + ${t1}$ :: res, s1, y + ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ + + if(x == zero) then + + res = gamma(p) + + else if(x > p) then + + s1 = -x + p * log(x) + res = gpx(p, x) * exp(s1) + + else if(x <= p .and. x > zero) then + + y = log_gamma(p) + s1 = -x + p * log(x) - y + res = (one - gpx(p, x) * exp(s1)) * exp(y) + + else + + + call error_stop("Error(Logarithm of upper incomplete gamma " & + //"function): negative x must be with integer p") + + end if + end function ingamma_up_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + impure elemental function ingamma_up_${t1[0]}$${k1}$${k2}$(p, x) & + result(res) + ! + ! Approximation of upper incomplete gamma function with integer p. + ! + ${t1}$, intent(in) :: p + ${t2}$, intent(in) :: x + ${t2}$ :: res, s1, y + ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$ + + if(x == zero) then + + res = gamma(real(p, ${k2}$)) + + else if(x > real(p, ${k2}$)) then + + s1 = -x + p * log(x) + res = gpx(p, x) * exp(s1) + + else if(x <= real(p, ${k2}$) .and. x > zero) then + + y = l_gamma(p, one) + s1 = -x + p * log(x) - y + res = gpx(p, x) * exp(s1) + res = (one - res) * exp(y) + + else + + y = l_gamma(p, one) + s1 = -x + p * log(abs(x)) - y + res = gpx(p, x) * exp(s1) + res = (one - (-1) ** p * res) * exp(y) + + end if + end function ingamma_up_${t1[0]}$${k1}$${k2}$ + + #:endfor + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function l_ingamma_up_${t1[0]}$${k1}$(p, x) result(res) + + ${t1}$, intent(in) :: p, x + ${t1}$ :: res, s1, y + ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ + + + if(x == zero) then + + res = log_gamma(p) + + else if(x > p) then + + s1 = -x + p * log(x) + res = log(gpx(p, x)) + s1 + + else if(x <= p .and. x > zero) then + + y= log_gamma(p) + s1 = -x + p * log(x) - y + res = gpx(p, x) * exp(s1) + res = log(one - res) + y + + else + + call error_stop("Error(Logarithm of upper incomplete gamma " & + //"function): negative x must be with integer p") + + end if + end function l_ingamma_up_${t1[0]}$${k1}$ + + #:endfor + + + + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + impure elemental function l_ingamma_up_${t1[0]}$${k1}$${k2}$(p, x) & + result(res) + + ${t1}$, intent(in) :: p + ${t2}$, intent(in) :: x + ${t2}$ :: res, s1, y + ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$ + + if(x == zero) then + + res = l_gamma(p, one) + + else if(x > real(p, ${k2}$)) then + + s1 = -x + p * log(x) + res = log(gpx(p, x)) + s1 + + else if(x <= real(p, ${k2}$) .and. x > zero) then + + y = l_gamma(p, one) + s1 = -x + p * log(x) - y + res = gpx(p, x) * exp(s1) + res = log(one - res) + y + + else + + y = l_gamma(p, one) + s1 = -x + p * log(abs(x)) + log(gpx(p, x)) + res = (-1) ** p * exp(s1) + res = log(abs(exp(y) - res)) + + end if + end function l_ingamma_up_${t1[0]}$${k1}$${k2}$ + + #:endfor + #:endfor + + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function regamma_p_${t1[0]}$${k1}$(p, x) result(res) + ! + ! Approximation of regularized incomplete gamma function P(p,x) for real p + ! + ${t1}$, intent(in) :: p, x + ${t1}$ :: res, s1 + ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ + + if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_p" & + //" function is not defined at x < 0") + + + if(x == zero) then + + res = zero + + else if(x > p) then + + s1 = -x + p * log(x) - log_gamma(p) + res = one - exp(s1 + log(gpx(p,x))) + + else if(x <= p) then + + s1 = -x + p * log(abs(x)) - log_gamma(p) + res = exp(log(gpx(p, x)) + s1) + + end if + end function regamma_p_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + impure elemental function regamma_p_${t1[0]}$${k1}$${k2}$(p, x) result(res) + ! + ! Approximation of regularized incomplete gamma function P(p,x) for integer p + ! + ${t1}$, intent(in) :: p + ${t2}$, intent(in) :: x + ${t2}$ :: res, s1 + ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$ + + if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_p" & + //" function is not defined at x < 0") + + + if(x == zero) then + + res = zero + + else if(x > real(p, ${k2}$)) then + + s1 = -x + p * log(x) - l_gamma(p, one) + res = one - exp(s1 + log(gpx(p,x))) + + else if(x <= real(p, ${k2}$)) then + + s1 = -x + p * log(abs(x)) - l_gamma(p, one) + res = exp(log(gpx(p, x)) + s1) + + end if + end function regamma_p_${t1[0]}$${k1}$${k2}$ + + #:endfor + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function regamma_q_${t1[0]}$${k1}$(p, x) result(res) + ! + ! Approximation of regularized incomplete gamma function Q(p,x) for real p + ! + ${t1}$, intent(in) :: p, x + ${t1}$ :: res, s1 + ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$ + + if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_q" & + //" function is not defined at x < 0") + + + if(x == zero) then + + res = one + + else if(x > p) then + + s1 = -x + p * log(x) - log_gamma(p) + res = exp(s1 + log(gpx(p,x))) + + else if(x <= p) then + + s1 = -x + p * log(abs(x)) - log_gamma(p) + res = one - exp(log(gpx(p, x)) + s1) + + end if + end function regamma_q_${t1[0]}$${k1}$ + + #:endfor + + + + #:for k1, t1 in INT_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES + impure elemental function regamma_q_${t1[0]}$${k1}$${k2}$(p, x) result(res) + ! + ! Approximation of regularized incomplet gamma function Q(p,x) for integer p + ! + ${t1}$, intent(in) :: p + ${t2}$, intent(in) :: x + ${t2}$ :: res, s1 + ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$ + + if(x < zero) call error_stop("Error(regamma_q): Regularized gamma_q" & + //" function is not defined at x < 0") + + + if(x == zero) then + + res = one + + else if(x > real(p, ${k2}$)) then + + s1 = -x + p * log(x) - l_gamma(p, one) + res = exp(log(gpx(p,x)) + s1) + + elseif(x <= real(p, ${k2}$)) then + + s1 = -x + p * log(abs(x)) - l_gamma(p, one) + res = one - exp(s1 + log(gpx(p,x))) + + end if + end function regamma_q_${t1[0]}$${k1}$${k2}$ + + #:endfor + #:endfor + +end module stdlib_specialfunctions_gamma 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/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/specialfunctions/test_specialfunctions_gamma.fypp b/src/tests/specialfunctions/test_specialfunctions_gamma.fypp new file mode 100644 index 000000000..26421bded --- /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}$, 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}$, & + -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 diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 651286708..ff9d45063 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -1,4 +1,4 @@ -### Pre-process: .fpp -> .f90 via Fypp +#### Pre-process: .fpp -> .f90 via Fypp # Create a list of the files to be preprocessed set(fppFiles