move gamma and add tests

Author: heltonmc

src/Bessels.jl

@@ -40,6 +40,8 @@ include("recurrence.jl")
 include("misc.jl")
 include("Polynomials/besselj_polys.jl")
 include("asymptotics.jl")
+include("gamma.jl")
+
 #include("hankel.jl")
 
 end

src/besselj.jl

@@ -236,71 +236,3 @@ function log_besselj_small_arguments_orders(v, x::T) where T
     logout = -loggamma(v + 1) + fma(v, log(x/2), log(out))
     return exp(logout)
 end
-
-# For 0.0 <= x < 171.5
-# Mean ULP = 0.55
-# Max ULP = 2.4
-# Adapted from Cephes Mathematical Library (MIT license https://en.smath.com/view/CephesMathLibrary/license) by Stephen L. Moshier
-function gamma(x)
-    if x < 0
-        isinteger(x) && throw(DomainError(x, "NaN result for non-NaN input."))
-        xp1 = abs(x) + 1.0
-        return π / sinpi(xp1) / _gamma(xp1)
-    else
-        return _gamma(x)
-    end
-end
-function _gamma(x)
-    if x > 11.5
-        return large_gamma(x)
-    elseif x <= 11.5
-        return small_gamma(x)
-    elseif isnan(x)
-        return x
-    end
-end
-function large_gamma(x)
-    isinf(x) && return x
-    T = Float64
-    w = inv(x)
-    s = (
-        8.333333333333331800504e-2, 3.472222222230075327854e-3, -2.681327161876304418288e-3, -2.294719747873185405699e-4,
-        7.840334842744753003862e-4, 6.989332260623193171870e-5, -5.950237554056330156018e-4, -2.363848809501759061727e-5,
-        7.147391378143610789273e-4
-    )
-    w = w * evalpoly(w, s) + one(T)
-    # lose precision on following block
-    y = exp((x))
-    # avoid overflow
-    v = x^(0.5 * x - 0.25)
-    y = v * (v / y)
-
-    return SQ2PI(T) * y * w
-end
-function small_gamma(x)
-    T = Float64
-    P = (
-        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
-    )
-
-    z = one(T)
-    while x >= 3.0
-        x -= one(T)
-        z *= x
-    end
-    while x < 2.0
-        z /= x
-        x += one(T)
-    end
-
-    x -= T(2)
-    p = evalpoly(x, P)
-    q = evalpoly(x, Q)
-    return z * p / q
-end

src/gamma.jl

@@ -0,0 +1,64 @@
+# Adapted from Cephes Mathematical Library (MIT license https://en.smath.com/view/CephesMathLibrary/license) by Stephen L. Moshier
+function gamma(x)
+    if x < 0
+        isinteger(x) && throw(DomainError(x, "NaN result for non-NaN input."))
+        xp1 = abs(x) + 1.0
+        return π / sinpi(xp1) / _gamma(xp1)
+    else
+        return _gamma(x)
+    end
+end
+function _gamma(x)
+    if x > 11.5
+        return large_gamma(x)
+    elseif x <= 11.5
+        return small_gamma(x)
+    elseif isnan(x)
+        return x
+    end
+end
+function large_gamma(x)
+    isinf(x) && return x
+    T = Float64
+    w = inv(x)
+    s = (
+        8.333333333333331800504e-2, 3.472222222230075327854e-3, -2.681327161876304418288e-3, -2.294719747873185405699e-4,
+        7.840334842744753003862e-4, 6.989332260623193171870e-5, -5.950237554056330156018e-4, -2.363848809501759061727e-5,
+        7.147391378143610789273e-4
+    )
+    w = w * evalpoly(w, s) + one(T)
+    # lose precision on following block
+    y = exp((x))
+    # avoid overflow
+    v = x^(0.5 * x - 0.25)
+    y = v * (v / y)
+
+    return SQ2PI(T) * y * w
+end
+function small_gamma(x)
+    T = Float64
+    P = (
+        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
+    )
+
+    z = one(T)
+    while x >= 3.0
+        x -= one(T)
+        z *= x
+    end
+    while x < 2.0
+        z /= x
+        x += one(T)
+    end
+
+    x -= T(2)
+    p = evalpoly(x, P)
+    q = evalpoly(x, Q)
+    return z * p / q
+end

test/bessely_test.jl

@@ -68,8 +68,8 @@ y1_32 = bessely1.(Float32.(x))
 ## Tests for bessely
 
 ## test all numbers and orders for 0<nu<100
-x = [0.05, 0.1, 0.2, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.92, 0.95, 0.97, 0.99, 1.0, 1.01, 1.05, 1.1, 1.2, 1.4, 2.0]
-nu = [2, 4, 6, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100]
+x = [0.05, 0.1, 0.2, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.91, 0.92, 0.93, 0.95, 0.96, 0.97, 0.98, 0.99, 0.995, 0.999, 1.0, 1.001, 1.01, 1.05, 1.1, 1.2, 1.4, 1.6, 1.8, 1.9, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 10.0]
+nu = [2, 4, 6, 10, 15, 20, 25, 30, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 110, 125, 150, 175, 200]
 for v in nu, xx in x
     xx *= v
     sf = SpecialFunctions.bessely(BigFloat(v), BigFloat(xx))

test/gamma_test.jl

@@ -0,0 +1,3 @@
+x = rand(10000)*170
+@test SpecialFunctions.gamma.(BigFloat.(x)) ≈ Bessels.gamma.(x)
+@test SpecialFunctions.gamma.(BigFloat.(-x)) ≈ Bessels.gamma.(-x)

test/runtests.jl

@@ -6,3 +6,4 @@ import SpecialFunctions
 @time @testset "besselk" begin include("besselk_test.jl") end
 @time @testset "besselj" begin include("besselj_test.jl") end
 @time @testset "bessely" begin include("bessely_test.jl") end
+@time @testset "gamma" begin include("gamma_test.jl") end