Merge pull request #876 from cyrilgandon/fix-doubling-precision

fix(factorial): result should be of the default kind, not double default kind

Author: jvdp1

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,15 +539,15 @@ contains
 
         case (zero)
 
-            res = zero_k2
+            res = zero_dp
 
         case (one)
 
-            res = zero_k2
+            res = zero_dp
 
         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}$

test/specialfunctions/test_specialfunctions_gamma.fypp

@@ -89,36 +89,33 @@ 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,    &
-                                         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(sp), 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(sp), 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(sp), 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
 
             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}$