From e8f78efb5537fe7b0497ecb9706b03063f028d31 Mon Sep 17 00:00:00 2001
From: Michael Helton <heltonmc@protonmail.com>
Date: Tue, 25 Oct 2022 12:08:50 -0400
Subject: [PATCH 1/2] dispatch int gamma

---
 src/besselj.jl     | 4 ++--
 src/gamma.jl       | 3 +--
 test/gamma_test.jl | 3 +++
 3 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/besselj.jl b/src/besselj.jl
index d8a1843..976e3c1 100644
--- a/src/besselj.jl
+++ b/src/besselj.jl
@@ -348,10 +348,10 @@ In general, this is most accurate for small arguments and when nu > x.
 function besselj_power_series(v, x::T) where T
     MaxIter = 3000
     S = promote_type(T, Float64)
-    v, x = S(v), S(x)
+    x = S(x)
 
     out = zero(S)
-    a = (x/2)^v / gamma(v + one(S))
+    a = (x/2)^v / gamma(v + 1)
     t2 = (x/2)^2
     for i in 0:MaxIter
         out += a
diff --git a/src/gamma.jl b/src/gamma.jl
index a00fc43..4b38547 100644
--- a/src/gamma.jl
+++ b/src/gamma.jl
@@ -1,5 +1,5 @@
 # Adapted from Cephes Mathematical Library (MIT license https://en.smath.com/view/CephesMathLibrary/license) by Stephen L. Moshier
-gamma(z::Number) = _gamma(float(z))
+gamma(x::Float64) = _gamma(x)
 _gamma(x::Float32) = Float32(_gamma(Float64(x)))
 
 function _gamma(x::Float64)
@@ -50,7 +50,6 @@ function _gamma(x::Float64)
     return z * p / q
 end
 
-
 function gamma(n::Integer)
     n < 0 && throw(DomainError(n, "`n` must not be negative."))
     n == 0 && return Inf*float(n)
diff --git a/test/gamma_test.jl b/test/gamma_test.jl
index 85052ac..8dfd9c7 100644
--- a/test/gamma_test.jl
+++ b/test/gamma_test.jl
@@ -3,3 +3,6 @@ x = rand(10000)*170
 @test SpecialFunctions.gamma.(BigFloat.(-x)) ≈ Bessels.gamma.(-x)
 @test isnan(Bessels.gamma(NaN))
 @test isinf(Bessels.gamma(Inf))
+
+x = [0, 1, 2, 3, 8, 15, 20, 30]
+@test SpecialFunctions.gamma.(x) ≈ Bessels.gamma.(x)

From 5ae131e05648cc5bcf33577abd216c90709d5991 Mon Sep 17 00:00:00 2001
From: Michael Helton <heltonmc@protonmail.com>
Date: Wed, 2 Nov 2022 12:06:18 -0400
Subject: [PATCH 2/2] fix dispatch

---
 src/besselj.jl | 4 ++--
 src/gamma.jl   | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/besselj.jl b/src/besselj.jl
index 976e3c1..d8a1843 100644
--- a/src/besselj.jl
+++ b/src/besselj.jl
@@ -348,10 +348,10 @@ In general, this is most accurate for small arguments and when nu > x.
 function besselj_power_series(v, x::T) where T
     MaxIter = 3000
     S = promote_type(T, Float64)
-    x = S(x)
+    v, x = S(v), S(x)
 
     out = zero(S)
-    a = (x/2)^v / gamma(v + 1)
+    a = (x/2)^v / gamma(v + one(S))
     t2 = (x/2)^2
     for i in 0:MaxIter
         out += a
diff --git a/src/gamma.jl b/src/gamma.jl
index 4b38547..cc85c3e 100644
--- a/src/gamma.jl
+++ b/src/gamma.jl
@@ -1,6 +1,6 @@
 # Adapted from Cephes Mathematical Library (MIT license https://en.smath.com/view/CephesMathLibrary/license) by Stephen L. Moshier
 gamma(x::Float64) = _gamma(x)
-_gamma(x::Float32) = Float32(_gamma(Float64(x)))
+gamma(x::Float32) = Float32(_gamma(Float64(x)))
 
 function _gamma(x::Float64)
     T = Float64
@@ -52,7 +52,7 @@ end
 
 function gamma(n::Integer)
     n < 0 && throw(DomainError(n, "`n` must not be negative."))
-    n == 0 && return Inf*float(n)
+    n == 0 && return Inf*one(n)
     n > 20 && return gamma(float(n))
     @inbounds return Float64(factorial(n-1))
 end