fix: special gamma functions

Author: wassup05

src/stdlib_specialfunctions_gamma.fypp

@@ -575,9 +575,9 @@ contains
     ! Fortran 90 program by Jim-215-Fisher
     !
         ${t1}$, intent(in) :: p, x
-        integer :: n, m
+        integer :: n
 
-        ${t2}$ :: res, p_lim, a, b, g, c, d, y, ss
+        ${t2}$ :: res, p_lim, a, b, g, c, d, y
         ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$
         ${t2}$, parameter :: dm = tiny(1.0_${k2}$) * 10 ** 6
         ${t1}$, parameter :: zero_k1 = 0.0_${k1}$
@@ -603,6 +603,9 @@ contains
             call error_stop("Error(gpx): Incomplete gamma function with "      &
             //"negative x must come with a whole number p not too small")
 
+        if(x < zero_k1) call error_stop("Error(gpx): Incomplete gamma"         &
+            // " function with negative x must have an integer parameter p")
+
         if(p >= p_lim) then     !use modified Lentz method of continued fraction
                                 !for eq. (15) in the above reference.
             a = one
@@ -670,7 +673,6 @@ contains
 
         else                            !Algorithm 2 in the reference
 
-            m = nint(ss)
             a = - x
             c = one / a
             d = p - one
@@ -689,9 +691,9 @@ contains
 
             end do
 
-            if(y >= b * tol_${k2}$ .and. mod(m , 2) /= 0) b = b + d * c / a
+            if(y >= b * tol_${k2}$ .and. mod(int(p), 2) /= 0) b = b + d * c / a
 
-            g = ((-1) ** m * exp(-a + log_gamma(p) - (p - 1) * log(a)) + b) / a
+            g = ((-1) ** p * exp(-a + log_gamma(p) - (p - 1) * log(a)) + b) / a
         end if
 
         res = g