From e3ee922286410ce6962ad6f88882e6185cdd106c Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 27 Mar 2025 09:27:58 +0100 Subject: [PATCH 01/12] gamma: generalize internal calculation to "next" more accurate precision --- src/stdlib_specialfunctions_gamma.fypp | 130 +++++++++++-------------- 1 file changed, 58 insertions(+), 72 deletions(-) diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp index d051c94df..a011d978b 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] + #: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,9 +325,6 @@ contains #:endfor - - - #:for k1, t1 in INT_KINDS_TYPES impure elemental function l_gamma_${t1[0]}$${k1}$(z) result(res) ! @@ -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] + #: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 @@ -557,14 +545,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] + #: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 +671,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. @@ -824,7 +810,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 +847,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 +887,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 +924,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 +956,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 +994,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 +1036,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 +1074,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 +1115,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 +1150,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 +1186,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 +1221,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 From 1d924489b3bdcaafe6a5078e4261d764e52dc3e8 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 27 Mar 2025 12:46:43 +0100 Subject: [PATCH 02/12] fix typo: `z_limit` is `real` --- src/stdlib_specialfunctions_gamma.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp index a011d978b..9efbbe99c 100644 --- a/src/stdlib_specialfunctions_gamma.fypp +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -424,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) From 9b253b19b987e1afc4b24ad747a245cd5562628e Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 28 Mar 2025 08:22:24 +0100 Subject: [PATCH 03/12] generalize tests --- .../test_specialfunctions_gamma.fypp | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index c29cf7a60..96adef261 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -1,10 +1,10 @@ #: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_specialfunctions_gamma, only: gamma, log_gamma, log_factorial, & lower_incomplete_gamma, & upper_incomplete_gamma, & @@ -18,7 +18,7 @@ module test_specialfunctions_gamma public :: collect_specialfunctions_gamma - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES ${t1}$, parameter :: tol_${k1}$ = 1000 * epsilon(1.0_${k1}$) #:endfor @@ -35,7 +35,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 +43,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 +59,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}$", & @@ -124,7 +124,7 @@ contains - #:for k1, t1 in CI_KINDS_TYPES + #:for k1, t1 in CI_KINDS_TYPES[:-1] subroutine test_gamma_${t1[0]}$${k1}$(error) type(error_type), allocatable, intent(out) :: error @@ -262,7 +262,7 @@ contains #:for k1, t1 in INT_KINDS_TYPES - #:for k2, t2 in R_KINDS_TYPES + #:for k2, t2 in REAL_KINDS_TYPES[:-1] subroutine test_lincgamma_${t1[0]}$${k1}$${k2}$(error) type(error_type), allocatable, intent(out) :: error @@ -411,7 +411,7 @@ contains - #:for k1, t1 in R_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES[:-1] subroutine test_lincgamma_${t1[0]}$${k1}$(error) type(error_type), allocatable, intent(out) :: error From a12c18d259f7eae4ff30a4a6a1a7ca4c9181b025 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 28 Mar 2025 08:31:57 +0100 Subject: [PATCH 04/12] temporary: verbose test result for CI error --- .../test_specialfunctions_gamma.fypp | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index 96adef261..7d1cee5cf 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -5,6 +5,7 @@ module test_specialfunctions_gamma use testdrive, only : new_unittest, unittest_type, error_type, check 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, & @@ -113,7 +114,7 @@ contains #:endif do i = 1, n - + call check(error, log_factorial(x(i)), ans(i), "Integer kind " & //"${k1}$ failed", thr = tol_dp, rel = .true.) @@ -130,6 +131,7 @@ contains type(error_type), allocatable, intent(out) :: error integer, parameter :: n = 5 integer :: i + type(state_type) :: err #:if k1 == "int8" @@ -182,13 +184,20 @@ contains #:elif t1[0] == "c" do i = 1, n + + err = state_type(STDLIB_VALUE_ERROR,'Complex ${k1}$ failed: x=',x(i),' gamma=',gamma(x(i)), & + 'expected=',ans(i),' tol=',tol_${k1}$) - call check(error, gamma(x(i)), ans(i), "Complex kind ${k1}$ failed",& + call check(error, gamma(x(i)), ans(i), err%print(),& thr = tol_${k1}$, rel = .true.) - end do + !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}$ From 5022e2554a3869fed2c5f8a8c29ed28d0603e2e6 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 28 Mar 2025 14:14:58 +0100 Subject: [PATCH 05/12] for all non-`sp` kinds, use max available precision --- src/stdlib_specialfunctions_gamma.fypp | 6 +++--- .../test_specialfunctions_gamma.fypp | 18 +++++++++++------- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp index 9efbbe99c..3b27d095d 100644 --- a/src/stdlib_specialfunctions_gamma.fypp +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -219,7 +219,7 @@ contains #! 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] + #: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 @@ -408,7 +408,7 @@ contains #! 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] + #: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) ! @@ -549,7 +549,7 @@ contains #! 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] + #: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) ! diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index 7d1cee5cf..4848113f2 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -185,15 +185,14 @@ contains do i = 1, n - err = state_type(STDLIB_VALUE_ERROR,'Complex ${k1}$ failed: x=',x(i),' gamma=',gamma(x(i)), & - 'expected=',ans(i),' tol=',tol_${k1}$) + err = state_type(STDLIB_VALUE_ERROR,'Complex ${k1}$ failed: x=',x(i), & + ' gamma=',gamma(x(i)), & + 'expected=',ans(i), & + ' tol=',tol_${k1}$) call check(error, gamma(x(i)), ans(i), err%print(),& thr = tol_${k1}$, rel = .true.) - !call check(error, gamma(x(i)), ans(i), "Complex kind ${k1}$ failed",& - ! thr = tol_${k1}$, rel = .true.) - end do #:endif @@ -257,8 +256,13 @@ contains do i = 1, n - call check(error, log_gamma(x(i)), ans(i), "Complex kind ${k1}$ " & - //"failed", thr = tol_${k1}$, rel = .true.) + err = state_type(STDLIB_VALUE_ERROR,'Complex ${k1}$ failed: x=',x(i), & + ' log(gamma)=',log_gamma(x(i)), & + 'expected=',ans(i), & + ' tol=',tol_${k1}$) + + call check(error, log_gamma(x(i)), ans(i), err%print(),& + thr = tol_${k1}$, rel = .true.) end do From 0ba315568227b0f47577b77d8823456cc0ba1791 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 28 Mar 2025 14:16:15 +0100 Subject: [PATCH 06/12] Update test_specialfunctions_gamma.fypp --- test/specialfunctions/test_specialfunctions_gamma.fypp | 1 + 1 file changed, 1 insertion(+) diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index 4848113f2..29cbee91d 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -205,6 +205,7 @@ contains type(error_type), allocatable, intent(out) :: error integer, parameter :: n = 4 integer :: i + type(state_type) :: err #:if k1 == "int8" From 43944cee6e130c894ac2bde2dfda7e112ffc8687 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 28 Mar 2025 14:21:33 +0100 Subject: [PATCH 07/12] Update test_specialfunctions_gamma.fypp --- test/specialfunctions/test_specialfunctions_gamma.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index 29cbee91d..2bd0bc9de 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -20,7 +20,7 @@ module test_specialfunctions_gamma public :: collect_specialfunctions_gamma #:for k1, t1 in REAL_KINDS_TYPES - ${t1}$, parameter :: tol_${k1}$ = 1000 * epsilon(1.0_${k1}$) + ${t1}$, parameter :: tol_${k1}$ = sqrt(epsilon(1.0_${k1}$)) #:endfor contains From e0e998f533a63fac79096bc5592ca34fd2926e6e Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 30 Mar 2025 10:16:15 +0200 Subject: [PATCH 08/12] Update test/specialfunctions/test_specialfunctions_gamma.fypp Co-authored-by: Jeremie Vandenplas --- test/specialfunctions/test_specialfunctions_gamma.fypp | 1 - 1 file changed, 1 deletion(-) diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index 2bd0bc9de..3f6127cbd 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -114,7 +114,6 @@ contains #:endif do i = 1, n - call check(error, log_factorial(x(i)), ans(i), "Integer kind " & //"${k1}$ failed", thr = tol_dp, rel = .true.) From 1dd60534ca853e690fd5f0005317541c8db4769a Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 30 Mar 2025 10:16:22 +0200 Subject: [PATCH 09/12] Update test/specialfunctions/test_specialfunctions_gamma.fypp Co-authored-by: Jeremie Vandenplas --- test/specialfunctions/test_specialfunctions_gamma.fypp | 1 - 1 file changed, 1 deletion(-) diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index 3f6127cbd..7849d2b65 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -193,7 +193,6 @@ contains thr = tol_${k1}$, rel = .true.) end do - #:endif end subroutine test_gamma_${t1[0]}$${k1}$ From 4ab2c810ace685f5579f25f0c48b6a19bc5f1205 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 30 Mar 2025 10:21:31 +0200 Subject: [PATCH 10/12] fix `integer`->`real` warning --- src/stdlib_specialfunctions_gamma.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp index 3b27d095d..fb3a9802a 100644 --- a/src/stdlib_specialfunctions_gamma.fypp +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -718,7 +718,7 @@ contains if(mod(n, 2) == 0) then - a = (1 - p - n / 2) * x + a = (one - p - n / 2) * x else From 220b79f695cb632289347f0e2d85b1e2cf224b07 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 30 Mar 2025 10:33:20 +0200 Subject: [PATCH 11/12] generalize `l_gamma`, `l_factorial` kinds --- src/stdlib_specialfunctions_gamma.fypp | 15 +-- .../test_specialfunctions_gamma.fypp | 113 +++++++----------- 2 files changed, 54 insertions(+), 74 deletions(-) diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp index fb3a9802a..9ee722b2d 100644 --- a/src/stdlib_specialfunctions_gamma.fypp +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -331,7 +331,7 @@ contains ! 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}$ @@ -350,7 +350,7 @@ contains do i = one, z - one - res = res + log(real(i)) + res = res + log(real(i,dp)) end do @@ -512,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") @@ -528,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}$ diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index 7849d2b65..6249e51a3 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -80,42 +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