Float64 improvements

Author: oscardssmith

src/gamma.jl

@@ -1,12 +1,12 @@
 # Float64 version adapted from Cephes Mathematical Library (MIT license https://en.smath.com/view/CephesMathLibrary/license) by Stephen L. Moshier
-function gamma(x::Float64)
+function gamma(_x::Float64)
     T = Float64
+    x = _x
     if x < 0
-        (isinteger(x) || x == -Inf) && throw(DomainError(x, "NaN result for non-NaN input."))
-        xp1 = abs(x) + 1.0
-        return π / (sinpi(xp1) * _gamma(xp1))
+        s = sinpi(_x)
+        s == 0 && throw(DomainError(x, "NaN result for non-NaN input."))
+        x = 1.0 - x
     end
-    isfinite(x) || return x
     if x > 11.5
         w = inv(x)
         s = (
@@ -17,34 +17,39 @@ function gamma(x::Float64)
         w = muladd(w, evalpoly(w, s), 1.0)
         # avoid overflow
         v = x ^ muladd(0.5, x, -0.25)
-        y = v * (v / exp(x))
 
-        return SQ2PI(T) * y * w
+        res = SQ2PI(T) * v * (v / exp(x)) * w
+	
+    	if _x < 0
+	    return π / (res * s)
+	else
+	    return res
+	end
     end
     P = (
-        1.000000000000000000009e0, 8.378004301573126728826e-1, 3.629515436640239168939e-1, 1.113062816019361559013e-1,
-        2.385363243461108252554e-2, 4.092666828394035500949e-3, 4.542931960608009155600e-4, 4.212760487471622013093e-5
+	1.000000000000000000009e0, 8.378004301573126728826e-1, 3.629515436640239168939e-1, 1.113062816019361559013e-1,
+	2.385363243461108252554e-2, 4.092666828394035500949e-3, 4.542931960608009155600e-4, 4.212760487471622013093e-5
     )
     Q = (
-        9.999999999999999999908e-1, 4.150160950588455434583e-1, -2.243510905670329164562e-1, -4.633887671244534213831e-2,
-        2.773706565840072979165e-2, -7.955933682494738320586e-4, -1.237799246653152231188e-3, 2.346584059160635244282e-4,
-        -1.397148517476170440917e-5
+	9.999999999999999999908e-1, 4.150160950588455434583e-1, -2.243510905670329164562e-1, -4.633887671244534213831e-2,
+	2.773706565840072979165e-2, -7.955933682494738320586e-4, -1.237799246653152231188e-3, 2.346584059160635244282e-4,
+	-1.397148517476170440917e-5
     )
 
     z = 1.0
     while x >= 3.0
-        x -= 1.0
-        z *= x
+	x -= 1.0
+	z *= x
     end
     while x < 2.0
-        z /= x
-        x += 1.0
+	z /= x
+	x += 1.0
     end
 
     x -= 2.0
     p = evalpoly(x, P)
     q = evalpoly(x, Q)
-    return z * p / q
+    retrun _x < 0 ? π * q / (s * z * p) : z * p / q
 end
 
 function gamma(_x::Float32)
@@ -89,7 +94,6 @@ function gamma(_x::Float16)
     return Float16(_x < 0 ? Float32(π)*den / (s*z*num) : z * num / den)
 end
 
-_gamma(x) = gamma(x) #easier than fixing this in other places.
 function gamma(n::Integer)
     n < 0 && throw(DomainError(n, "`n` must not be negative."))
     n == 0 && return Inf*one(n)