Merge pull request #730 from banana-bred/master

PR related to #729, gamma function of pure imaginary returns correct value

Author: jvdp1

src/stdlib_specialfunctions_gamma.fypp

@@ -298,14 +298,14 @@ contains
 
         end if
 
-        if(z % re > zero_k1) then
+        if(z % re < zero_k1) then
 
-            y = z - one
+            x = cmplx(abs(z % re), - z % im, kind = ${k1}$)
+            y = x - one
 
         else
 
-            x = cmplx(abs(z % re), - z % im, kind = ${k1}$)
-            y = x - one
+            y = z - one
 
         end if
 

test/specialfunctions/test_specialfunctions_gamma.fypp

@@ -96,21 +96,21 @@ contains
                                          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,    &
@@ -134,42 +134,46 @@ contains
 
     subroutine test_gamma_${t1[0]}$${k1}$(error)
         type(error_type), allocatable, intent(out) :: error
-        integer, parameter :: n = 4
+        integer, parameter :: n = 5
         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}$]
+        ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 4_${k1}$, 5_${k1}$, 6_${k1}$]
+        ${t1}$, parameter :: ans(n) = [1_${k1}$, 1_${k1}$, 6_${k1}$, 24_${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}$]
+        ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 4_${k1}$, 5_${k1}$, 8_${k1}$]
+        ${t1}$, parameter :: ans(n) = [1_${k1}$, 1_${k1}$, 6_${k1}$, 24_${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}$,           &
+        ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 4_${k1}$, 8_${k1}$, 13_${k1}$]
+        ${t1}$, parameter :: ans(n) = [1_${k1}$, 1_${k1}$, 6_${k1}$, 5040_${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}$,           &
+        ${t1}$, parameter :: x(n) = [1_${k1}$, 2_${k1}$, 4_${k1}$, 13_${k1}$, 21_${k1}$]
+        ${t1}$, parameter :: ans(n) = [1_${k1}$, 1_${k1}$, 6_${k1}$, 479001600_${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}$)]
+                                     (-1.254e1_${k1}$, -9.87_${k1}$),          &
+                                     (0.0_${k1}$, 1.0_${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}$)]
+             (-2.18767396709283064e-21_${k1}$, 2.77577940846953455e-21_${k1}$),&
+                    (-0.15494982830181067_${k1}$, -0.49801566811835607_${k1}$) &
+             ]
         #:endif