From c14c9d80f7b467349623ca4a25ea4d1e70b3f82d Mon Sep 17 00:00:00 2001 From: Cyril Gandon Date: Mon, 23 Sep 2024 09:44:43 +0200 Subject: [PATCH 1/4] fix(factorial): result should be of the default kind, not double default kind Fixes https://github.com/fortran-lang/stdlib/issues/875 --- 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 e91f14954..c8dc4c444 100644 --- a/src/stdlib_specialfunctions_gamma.fypp +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -547,7 +547,7 @@ contains case (two : ) - res = l_gamma(n + 1, 1.0D0) + res = l_gamma(n + 1, 1.0_dp) end select end function l_factorial_${t1[0]}$${k1}$ From a65cc2f2da275bb94a04bae3c5957d16cf018504 Mon Sep 17 00:00:00 2001 From: Cyril Gandon Date: Fri, 27 Sep 2024 08:27:59 +0200 Subject: [PATCH 2/4] Promote real to real_dp as it should be --- src/stdlib_specialfunctions_gamma.fypp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp index c8dc4c444..8fff966c2 100644 --- a/src/stdlib_specialfunctions_gamma.fypp +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -528,9 +528,9 @@ contains ! Log(n!) ! ${t1}$, intent(in) :: n - real :: res + real(dp) :: res ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$ - real, parameter :: zero_k2 = 0.0 + real(dp), parameter :: zero_dp = 0.0_dp if(n < zero) call error_stop("Error(l_factorial): Logarithm of" & //" factorial function argument must be non-negative") @@ -539,11 +539,11 @@ contains case (zero) - res = zero_k2 + res = zero_dp case (one) - res = zero_k2 + res = zero_dp case (two : ) From 2856a7383d73f1cc790d37796484198b894c712a Mon Sep 17 00:00:00 2001 From: Cyril Gandon Date: Fri, 27 Sep 2024 16:50:21 +0200 Subject: [PATCH 3/4] fix test --- test/specialfunctions/test_specialfunctions_gamma.fypp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index 62ee4c1f9..cd8b3bef1 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -89,28 +89,28 @@ contains ${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, & + real(dp), 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, & + real(dp), 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, & + real(dp), 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, & + real(dp), parameter :: ans(n) = [0.0, 0.0, 0.693147180, 3.17805383, & 4.23356165e1, 9.36687468e5] #:endif @@ -118,7 +118,7 @@ contains do i = 1, n call check(error, log_factorial(x(i)), ans(i), "Integer kind " & - //"${k1}$ failed", thr = tol_sp, rel = .true.) + //"${k1}$ failed", thr = tol_dp, rel = .true.) end do end subroutine test_logfact_${t1[0]}$${k1}$ From 6d8076f800dcdc393057b71e92bf55e6e6ef88fa Mon Sep 17 00:00:00 2001 From: Cyril Gandon Date: Mon, 30 Sep 2024 09:47:34 +0200 Subject: [PATCH 4/4] Add extra digits to fix new precision for test of log_factorial --- .../test_specialfunctions_gamma.fypp | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index cd8b3bef1..c29cf7a60 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -89,30 +89,27 @@ contains ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & 5_${k1}$, 100_${k1}$] - real(dp), parameter :: ans(n) = [0.0, 0.0, 0.693147180, 3.17805383, & - 4.78749174, 3.63739376e2] + 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, 0.0, 0.693147180, 3.17805383, & - 8.52516136, 2.61133046e3] - + 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, 0.0, 0.693147180, 3.17805383, & - 1.99872145e1, 5.49810038e4] - + 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, 0.0, 0.693147180, 3.17805383, & - 4.23356165e1, 9.36687468e5] - + real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, & + 4.233561646075e1_dp, 9.366874681600e5_dp] #:endif do i = 1, n