diff --git a/example/specialfunctions_gamma/example_gamma.f90 b/example/specialfunctions_gamma/example_gamma.f90 index 6df5f42e4..921283d71 100644 --- a/example/specialfunctions_gamma/example_gamma.f90 +++ b/example/specialfunctions_gamma/example_gamma.f90 @@ -1,5 +1,5 @@ program example_gamma - use stdlib_kinds, only: dp, int64 + use stdlib_kinds, only: sp, dp, int64 use stdlib_specialfunctions_gamma, only: gamma implicit none @@ -7,15 +7,13 @@ program example_gamma integer(int64) :: n real :: x real(dp) :: y - complex :: z - complex(dp) :: z1 + complex(sp) :: z 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 @@ -32,6 +30,4 @@ program example_gamma print *, gamma(z) ! (0.988054395, 0.383354813) - print *, gamma(z1) -! (-2.78916032990983999E-005, 9.83164600163221218E-006) end program example_gamma diff --git a/example/specialfunctions_gamma/example_log_gamma.f90 b/example/specialfunctions_gamma/example_log_gamma.f90 index 70b093670..08971c8a4 100644 --- a/example/specialfunctions_gamma/example_log_gamma.f90 +++ b/example/specialfunctions_gamma/example_log_gamma.f90 @@ -1,19 +1,18 @@ program example_log_gamma - use stdlib_kinds, only: dp + use stdlib_kinds, only: sp, dp use stdlib_specialfunctions_gamma, only: log_gamma implicit none integer :: i real :: x real(dp) :: y - complex :: z - complex(dp) :: z1 + complex(sp) :: z i = 10 x = 8.76 y = x z = (5.345, -3.467) - z1 = z + print *, log_gamma(i) !default single precision output !12.8018274 @@ -29,7 +28,4 @@ program example_log_gamma !(2.56165648, -5.73382425) - print *, log_gamma(z1) - -!(2.5616575105114614, -5.7338247782852498) end program example_log_gamma diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp index d051c94df..9ee722b2d 100644 --- a/src/stdlib_specialfunctions_gamma.fypp +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -1,28 +1,24 @@ #:include "common.fypp" -#:set R_KINDS_TYPES = [KT for KT in REAL_KINDS_TYPES if KT[0] in ["sp","dp"]] -#:set C_KINDS_TYPES = [KT for KT in CMPLX_KINDS_TYPES if KT[0] in ["sp","dp"]] -#:set CI_KINDS_TYPES = INT_KINDS_TYPES + C_KINDS_TYPES +#:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_KINDS_TYPES +#:set IDX_CMPLX_KINDS_TYPES = [(i, CMPLX_KINDS[i], CMPLX_TYPES[i], CMPLX_INIT[i]) for i in range(len(CMPLX_KINDS))] +#:set IDX_REAL_KINDS_TYPES = [(i, REAL_KINDS[i], REAL_TYPES[i], REAL_INIT[i]) for i in range(len(REAL_KINDS))] module stdlib_specialfunctions_gamma - use iso_fortran_env, only : qp => real128 use ieee_arithmetic, only: ieee_value, ieee_quiet_nan - use stdlib_kinds, only : sp, dp, int8, int16, int32, int64 + 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 = 6_int8 + integer(int8), parameter :: max_fact_int8 = 6_int8 integer(int16), parameter :: max_fact_int16 = 8_int16 integer(int32), parameter :: max_fact_int32 = 13_int32 integer(int64), parameter :: max_fact_int64 = 21_int64 - #:for k1, t1 in R_KINDS_TYPES + #: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 @@ -33,7 +29,7 @@ module stdlib_specialfunctions_gamma interface gamma !! Gamma function for integer and complex numbers !! - #:for k1, t1 in CI_KINDS_TYPES + #:for k1, t1 in CI_KINDS_TYPES[:-1] module procedure gamma_${t1[0]}$${k1}$ #:endfor end interface gamma @@ -43,7 +39,7 @@ module stdlib_specialfunctions_gamma interface log_gamma !! Logarithm of gamma function !! - #:for k1, t1 in CI_KINDS_TYPES + #:for k1, t1 in CI_KINDS_TYPES[:-1] module procedure l_gamma_${t1[0]}$${k1}$ #:endfor end interface log_gamma @@ -64,12 +60,12 @@ module stdlib_specialfunctions_gamma !! Lower incomplete gamma function !! #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] module procedure ingamma_low_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] module procedure ingamma_low_${t1[0]}$${k1}$ #:endfor end interface lower_incomplete_gamma @@ -80,12 +76,12 @@ module stdlib_specialfunctions_gamma !! Logarithm of lower incomplete gamma function !! #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] module procedure l_ingamma_low_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] module procedure l_ingamma_low_${t1[0]}$${k1}$ #:endfor end interface log_lower_incomplete_gamma @@ -96,12 +92,12 @@ module stdlib_specialfunctions_gamma !! Upper incomplete gamma function !! #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] module procedure ingamma_up_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] module procedure ingamma_up_${t1[0]}$${k1}$ #:endfor end interface upper_incomplete_gamma @@ -112,12 +108,12 @@ module stdlib_specialfunctions_gamma !! Logarithm of upper incomplete gamma function !! #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] module procedure l_ingamma_up_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] module procedure l_ingamma_up_${t1[0]}$${k1}$ #:endfor end interface log_upper_incomplete_gamma @@ -128,12 +124,12 @@ module stdlib_specialfunctions_gamma !! Regularized (normalized) lower incomplete gamma function, P !! #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] module procedure regamma_p_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] module procedure regamma_p_${t1[0]}$${k1}$ #:endfor end interface regularized_gamma_p @@ -144,12 +140,12 @@ module stdlib_specialfunctions_gamma !! Regularized (normalized) upper incomplete gamma function, Q !! #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] module procedure regamma_q_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] module procedure regamma_q_${t1[0]}$${k1}$ #:endfor end interface regularized_gamma_q @@ -160,12 +156,12 @@ module stdlib_specialfunctions_gamma ! Incomplete gamma G function. ! Internal use only ! - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] module procedure gpx_${t1[0]}$${k1}$ !for real p and x #:endfor #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] module procedure gpx_${t1[0]}$${k1}$${k2}$ !for integer p and real x #:endfor #:endfor @@ -178,7 +174,7 @@ module stdlib_specialfunctions_gamma ! Internal use only ! #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] module procedure l_gamma_${t1[0]}$${k1}$${k2}$ #:endfor #:endfor @@ -219,14 +215,12 @@ contains - #:for k1, t1 in C_KINDS_TYPES - #:if k1 == "sp" - #:set k2 = "dp" - #:elif k1 == "dp" - #:set k2 = "qp" - #:endif - #:set t2 = "real({})".format(k2) - + #! Because the KIND lists are sorted by increasing accuracy, + #! gamma will use the next available more accurate KIND for the + #! internal more accurate solver. + #:for i, k1, t1, i1 in IDX_CMPLX_KINDS_TYPES[:-1] + #:set k2 = CMPLX_KINDS[i + 1] if k1 == "sp" else CMPLX_KINDS[-1] + #:set t2 = "real({})".format(k2) impure elemental function gamma_${t1[0]}$${k1}$(z) result(res) ${t1}$, intent(in) :: z ${t1}$ :: res @@ -255,8 +249,8 @@ contains -2.71994908488607704e-9_${k2}$] ! parameters from above referenced source. - #:elif k1 == "dp" - #! for double precision input, using quadruple precision for calculation + #:else + #! for double or extended precision input, using quadruple precision for calculation integer, parameter :: n = 24 ${t2}$, parameter :: r = 25.617904_${k2}$ @@ -290,8 +284,6 @@ contains #:endif - - if(abs(z % im) < tol_${k1}$) then res = cmplx(gamma(z % re), kind = ${k1}$) @@ -333,16 +325,13 @@ contains #: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 + real(dp) :: res ${t1}$ :: i ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$ @@ -361,7 +350,7 @@ contains do i = one, z - one - res = res + log(real(i)) + res = res + log(real(i,dp)) end do @@ -374,7 +363,7 @@ contains #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] impure elemental function l_gamma_${t1[0]}$${k1}$${k2}$(z, x) result(res) ! @@ -415,13 +404,12 @@ contains - #:for k1, t1 in C_KINDS_TYPES - #:if k1 == "sp" - #:set k2 = "dp" - #:elif k1 == "dp" - #:set k2 = "qp" - #:endif - #:set t2 = "real({})".format(k2) + #! Because the KIND lists are sorted by increasing accuracy, + #! gamma will use the next available more accurate KIND for the + #! internal more accurate solver. + #:for i, k1, t1, i1 in IDX_CMPLX_KINDS_TYPES[:-1] + #:set k2 = CMPLX_KINDS[i + 1] if k1 == "sp" else CMPLX_KINDS[-1] + #: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 @@ -436,7 +424,7 @@ contains real(${k1}$) :: d integer :: m, i complex(${k2}$) :: zr, zr2, sum, s - real(${k1}$), parameter :: z_limit = 10_${k1}$, zero_k1 = 0.0_${k1}$ + real(${k1}$), parameter :: z_limit = 10.0_${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) @@ -524,14 +512,15 @@ contains #:for k1, t1 in INT_KINDS_TYPES + #:set k2, t2 = REAL_KINDS[-2], REAL_TYPES[-2] impure elemental function l_factorial_${t1[0]}$${k1}$(n) result(res) ! ! Log(n!) ! ${t1}$, intent(in) :: n - real(dp) :: res + ${t2}$ :: res ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$ - real(dp), parameter :: zero_dp = 0.0_dp + ${t2}$, parameter :: zero_${k2}$ = 0.0_${k2}$, one_${k2}$ = 1.0_${k2}$ if(n < zero) call error_stop("Error(l_factorial): Logarithm of" & //" factorial function argument must be non-negative") @@ -540,15 +529,15 @@ contains case (zero) - res = zero_dp + res = zero_${k2}$ case (one) - res = zero_dp + res = zero_${k2}$ case (two : ) - res = l_gamma(n + 1, 1.0_dp) + res = l_gamma(n + 1, one_${k2}$) end select end function l_factorial_${t1[0]}$${k1}$ @@ -557,14 +546,12 @@ contains - #:for k1, t1 in R_KINDS_TYPES - #:if k1 == "sp" - #:set k2 = "dp" - #:elif k1 == "dp" - #:set k2 = "qp" - #:endif - #:set t2 = "real({})".format(k2) - + #! Because the KIND lists are sorted by increasing accuracy, + #! gamma will use the next available more accurate KIND for the + #! internal more accurate solver. + #:for i, k1, t1, i1 in IDX_REAL_KINDS_TYPES[:-1] + #:set k2 = REAL_KINDS[i + 1] if k1 == "sp" else REAL_KINDS[-1] + #:set t2 = REAL_TYPES[i + 1] impure elemental function gpx_${t1[0]}$${k1}$(p, x) result(res) ! ! Approximation of incomplete gamma G function with real argument p. @@ -685,7 +672,7 @@ contains #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] impure elemental function gpx_${t1[0]}$${k1}$${k2}$(p, x) result(res) ! ! Approximation of incomplete gamma G function with integer argument p. @@ -732,7 +719,7 @@ contains if(mod(n, 2) == 0) then - a = (1 - p - n / 2) * x + a = (one - p - n / 2) * x else @@ -824,7 +811,7 @@ contains - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] impure elemental function ingamma_low_${t1[0]}$${k1}$(p, x) result(res) ! ! Approximation of lower incomplete gamma function with real p. @@ -861,7 +848,7 @@ contains #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] impure elemental function ingamma_low_${t1[0]}$${k1}$${k2}$(p, x) & result(res) ! @@ -901,7 +888,7 @@ contains - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] impure elemental function l_ingamma_low_${t1[0]}$${k1}$(p, x) result(res) ${t1}$, intent(in) :: p, x @@ -938,7 +925,7 @@ contains #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] impure elemental function l_ingamma_low_${t1[0]}$${k1}$${k2}$(p, x) & result(res) @@ -970,7 +957,7 @@ contains - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] impure elemental function ingamma_up_${t1[0]}$${k1}$(p, x) result(res) ! ! Approximation of upper incomplete gamma function with real p. @@ -1008,7 +995,7 @@ contains #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] impure elemental function ingamma_up_${t1[0]}$${k1}$${k2}$(p, x) & result(res) ! @@ -1050,7 +1037,7 @@ contains - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] impure elemental function l_ingamma_up_${t1[0]}$${k1}$(p, x) result(res) ${t1}$, intent(in) :: p, x @@ -1088,7 +1075,7 @@ contains #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] impure elemental function l_ingamma_up_${t1[0]}$${k1}$${k2}$(p, x) & result(res) @@ -1129,7 +1116,7 @@ contains - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] impure elemental function regamma_p_${t1[0]}$${k1}$(p, x) result(res) ! ! Approximation of regularized incomplete gamma function P(p,x) for real p @@ -1164,7 +1151,7 @@ contains #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] 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 @@ -1200,7 +1187,7 @@ contains - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] impure elemental function regamma_q_${t1[0]}$${k1}$(p, x) result(res) ! ! Approximation of regularized incomplete gamma function Q(p,x) for real p @@ -1235,7 +1222,7 @@ contains #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] 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 diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index c29cf7a60..6249e51a3 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -1,10 +1,11 @@ #:include "common.fypp" -#:set R_KINDS_TYPES = [KT for KT in REAL_KINDS_TYPES if KT[0] in ["sp","dp"]] -#:set C_KINDS_TYPES = [KT for KT in CMPLX_KINDS_TYPES if KT[0] in ["sp","dp"]] -#:set CI_KINDS_TYPES = INT_KINDS_TYPES + C_KINDS_TYPES +#:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_KINDS_TYPES +#:set IDX_CMPLX_KINDS_TYPES = [(i, CMPLX_KINDS[i], CMPLX_TYPES[i], CMPLX_INIT[i]) for i in range(len(CMPLX_KINDS))] +#:set IDX_REAL_KINDS_TYPES = [(i, REAL_KINDS[i], REAL_TYPES[i], REAL_INIT[i]) for i in range(len(REAL_KINDS))] 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_kinds, only: sp, dp, xdp, qp, int8, int16, int32, int64 + use stdlib_error, only: state_type, STDLIB_VALUE_ERROR use stdlib_specialfunctions_gamma, only: gamma, log_gamma, log_factorial, & lower_incomplete_gamma, & upper_incomplete_gamma, & @@ -18,8 +19,8 @@ module test_specialfunctions_gamma public :: collect_specialfunctions_gamma - #:for k1, t1 in R_KINDS_TYPES - ${t1}$, parameter :: tol_${k1}$ = 1000 * epsilon(1.0_${k1}$) + #:for k1, t1 in REAL_KINDS_TYPES + ${t1}$, parameter :: tol_${k1}$ = sqrt(epsilon(1.0_${k1}$)) #:endfor contains @@ -35,7 +36,7 @@ contains test_logfact_${t1[0]}$${k1}$) & #:endfor - #:for k1, t1 in CI_KINDS_TYPES + #:for k1, t1 in CI_KINDS_TYPES[:-1] , new_unittest("gamma_${t1[0]}$${k1}$", & test_gamma_${t1[0]}$${k1}$) & , new_unittest("log_gamma_${t1[0]}$${k1}$", & @@ -43,7 +44,7 @@ contains #:endfor #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] , 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}$", & @@ -59,7 +60,7 @@ contains #:endfor #:endfor - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] , new_unittest("lower_incomplete_gamma_${t1[0]}$${k1}$", & test_lincgamma_${t1[0]}$${k1}$) & , new_unittest("log_lower_incomplete_gamma_${t1[0]}$${k1}$", & @@ -79,43 +80,32 @@ contains #:for k1, t1 in INT_KINDS_TYPES - + #:set k2, t2 = REAL_KINDS[-2], REAL_TYPES[-2] subroutine test_logfact_${t1[0]}$${k1}$(error) - type(error_type), allocatable, intent(out) :: error - integer, parameter :: n = 6 + type(error_type), allocatable, intent(out) :: error integer :: i - - #:if k1 == "int8" - - ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & - 5_${k1}$, 100_${k1}$] - real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, & - 4.78749174278204_dp, 3.637393755555e2_dp] - - #:elif k1 == "int16" - - ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & - 7_${k1}$, 500_${k1}$] - real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, & - 8.52516136106541_dp, 2.611330458460e3_dp] - #:elif k1 == "int32" - - ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & - 12_${k1}$, 7000_${k1}$] - real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, & - 1.998721449566e1_dp, 5.498100377941e4_dp] - #:elif k1 == "int64" - - ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & - 20_${k1}$, 90000_${k1}$] - real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, & - 4.233561646075e1_dp, 9.366874681600e5_dp] - #:endif + + integer, parameter :: xtest(*) = [0,1,2,4,5,7,12,20,100,500,7000,90000] + ${t2}$, parameter :: res(*) = [0.0_${k2}$, & + 0.0_${k2}$, & + 0.69314718055994_${k2}$, & + 3.17805383034794_${k2}$, & + 4.78749174278204_${k2}$, & + 8.52516136106541_${k2}$, & + 1.998721449566e1_${k2}$, & + 4.233561646075e1_${k2}$, & + 3.637393755555e2_${k2}$, & + 2.611330458460e3_${k2}$, & + 5.498100377941e4_${k2}$, & + 9.366874681600e5_${k2}$] + + ${t1}$, parameter :: x(*) = pack(xtest, xtest