Refactor code into separate modules

src/AiryFunctions/AiryFunctions.jl

@@ -0,0 +1,20 @@
+module AiryFunctions
+
+using ..GammaFunctions
+using ..Bessels: ComplexOrReal
+using ..Math
+
+export airyai
+export airyaix
+export airyaiprime
+export airyaiprimex
+export airybi
+export airybix
+export airybiprime
+export airybiprimex
+
+include("airy.jl")
+include("cairy.jl")
+include("airy_polys.jl")
+
+end

src/Airy/cairy.jl → src/AiryFunctions/cairy.jl

@@ -21,6 +21,8 @@
 # [1] Jentschura, Ulrich David, and E. Lötstedt. "Numerical calculation of Bessel, Hankel and Airy functions." 
 #     Computer Physics Communications 183.3 (2012): 506-519.
 
+using ..BesselFunctions: besseli_power_series_inu_inum1, besselk_down_recurrence, besselk_ratio_knu_knup1
+
 """
     airyai(z)
 

src/BesselFunctions/BesselFunctions.jl

@@ -0,0 +1,61 @@
+module BesselFunctions
+
+using ..Bessels: ComplexOrReal
+using ..GammaFunctions
+using ..Math
+
+export besselj0
+export besselj1
+export besselj
+export besselj!
+
+export bessely0
+export bessely1
+export bessely
+export bessely!
+
+export sphericalbesselj
+export sphericalbessely
+export sphericalbesseli
+export sphericalbesselk
+
+export besseli
+export besselix
+export besseli0
+export besseli0x
+export besseli1
+export besseli1x
+export besseli!
+
+export besselk
+export besselkx
+export besselk0
+export besselk0x
+export besselk1
+export besselk1x
+export besselk!
+
+export besselh
+export hankelh1
+export hankelh2
+
+include("besseli.jl")
+include("besselj.jl")
+include("besselk.jl")
+include("bessely.jl")
+include("hankel.jl")
+include("sphericalbessel.jl")
+include("modifiedsphericalbessel.jl")
+
+
+include("constants.jl")
+include("U_polynomials.jl")
+include("recurrence.jl")
+include("Polynomials/besselj_polys.jl")
+include("asymptotics.jl")
+
+include("Float128/constants.jl")
+include("Float128/besselj.jl")
+include("Float128/bessely.jl")
+
+end

src/bessely.jl → src/BesselFunctions/bessely.jl

@@ -508,17 +508,6 @@ function bessely_chebyshev_low_orders(v, x)
     return clenshaw_chebyshev(v1, a), clenshaw_chebyshev(v2, a)
 end
 
-# uses the Clenshaw algorithm to recursively evaluate a linear combination of Chebyshev polynomials
-function clenshaw_chebyshev(x, c)
-    x2 = 2x
-    c0 = c[end-1]
-    c1 = c[end]
-    for i in length(c)-2:-1:1
-        c0, c1 = c[i] - c1, c0 + c1 * x2
-    end
-    return c0 + c1 * x
-end
-
 # to generate the Chebyshev weights
 #=
 using ArbNumerics, FastChebInterp

src/Bessels.jl

@@ -10,6 +10,8 @@ export bessely
 
 export sphericalbesselj
 export sphericalbessely
+export sphericalbesseli
+export sphericalbesselk
 
 export besseli
 export besselix
@@ -38,32 +40,18 @@ export airybix
 export airybiprime
 export airybiprimex
 
-const ComplexOrReal{T} = Union{T,Complex{T}}
+export gamma
 
-include("besseli.jl")
-include("besselj.jl")
-include("besselk.jl")
-include("bessely.jl")
-include("hankel.jl")
-include("Airy/airy.jl")
-include("Airy/cairy.jl")
-include("sphericalbessel.jl")
-include("modifiedsphericalbessel.jl")
+const ComplexOrReal{T} = Union{T,Complex{T}}
 
-include("Float128/besseli.jl")
-include("Float128/besselj.jl")
-include("Float128/besselk.jl")
-include("Float128/bessely.jl")
-include("Float128/constants.jl")
+include("Math/Math.jl")
+include("GammaFunctions/GammaFunctions.jl")
+include("BesselFunctions/BesselFunctions.jl")
+include("AiryFunctions/AiryFunctions.jl")
 
-include("constants.jl")
-include("math_constants.jl")
-include("U_polynomials.jl")
-include("recurrence.jl")
-include("misc.jl")
-include("Polynomials/besselj_polys.jl")
-include("asymptotics.jl")
-include("gamma.jl")
+using .GammaFunctions
+using .BesselFunctions
+using .AiryFunctions
 
 precompile(besselj, (Float64, Float64))
 

src/GammaFunctions/GammaFunctions.jl

@@ -0,0 +1,9 @@
+module GammaFunctions
+
+using ..Math: SQ2PI
+
+export gamma
+
+include("gamma.jl")
+
+end

src/misc.jl → src/Math/Math.jl

@@ -1,3 +1,12 @@
+module Math
+
+# constains miscelaneous math functions and math constants
+
+export cos_sum, sin_sum
+export clenshaw_chebyshev
+
+include("math_constants.jl")
+
 # function to more accurately compute cos(x + xn)
 # see https://github.com/heltonmc/Bessels.jl/pull/13
 # written by @oscardssmith
@@ -18,6 +27,7 @@ function cos_sum(x, xn)
         return Base.Math.sin_kernel(y)
     end
 end
+
 # function to more accurately compute sin(x + xn)
 function sin_sum(x, xn)
     s = x + xn
@@ -36,6 +46,7 @@ function sin_sum(x, xn)
         return -Base.Math.cos_kernel(y)
     end
 end
+
 @inline function reduce_pi02_med(x::Float64)
     pio2_1 = 1.57079632673412561417e+00
 
@@ -45,3 +56,16 @@ end
     y = r-w
     return unsafe_trunc(Int, fn), Base.Math.DoubleFloat64(y, (r-y)-w)
 end
+
+# uses the Clenshaw algorithm to recursively evaluate a linear combination of Chebyshev polynomials
+function clenshaw_chebyshev(x, c)
+    x2 = 2x
+    c0 = c[end-1]
+    c1 = c[end]
+    for i in length(c)-2:-1:1
+        c0, c1 = c[i] - c1, c0 + c1 * x2
+    end
+    return c0 + c1 * x
+end
+
+end

src/math_constants.jl → src/Math/math_constants.jl

@@ -1,3 +1,9 @@
+export ONEOSQPI, TWOOPI, PIO2, SQPIO2, SQ1O2PI, SQ2OPI
+export THPIO4, SQ2PI, PIPOW3O2, PIO4, SQ2O2
+
+export GAMMA_TWO_THIRDS, GAMMA_ONE_THIRD, GAMMA_ONE_SIXTH, GAMMA_FIVE_SIXTHS
+
+
 const ONEOSQPI(::Type{BigFloat}) = big"5.6418958354775628694807945156077258584405E-1"
 const TWOOPI(::Type{BigFloat}) = big"6.3661977236758134307553505349005744813784E-1"
 const PIO2(::Type{BigFloat}) = big"1.570796326794896619231321691639751442098584699687552910487472296153908203143099"

test/besselj_test.jl

@@ -102,8 +102,8 @@ for v in nu, xx in x
     xx *= v
     sf = SpecialFunctions.besselj(BigFloat(v), BigFloat(xx))
     @test isapprox(besselj(v, xx), Float64(sf), rtol=5e-14)
-    @test isapprox(Bessels.besseljy_positive_args(v, xx)[1], Float64(sf), rtol=5e-14)
-    @test isapprox(Bessels.besselj(Float32(v), Float32(xx)), Float32(sf))
+    @test isapprox(Bessels.BesselFunctions.besseljy_positive_args(v, xx)[1], Float64(sf), rtol=5e-14)
+    @test isapprox(besselj(Float32(v), Float32(xx)), Float32(sf))
 end
 
 # test half orders (SpecialFunctions does not give big float precision)
@@ -115,8 +115,8 @@ for v in nu, xx in x
     xx *= v
     sf = SpecialFunctions.besselj(v, xx)
     @test isapprox(besselj(v, xx), sf, rtol=1e-12)
-    @test isapprox(Bessels.besseljy_positive_args(v, xx)[1], sf, rtol=1e-12)
-    @test isapprox(Bessels.besselj(Float32(v), Float32(xx)), Float32(sf))
+    @test isapprox(Bessels.BesselFunctions.besseljy_positive_args(v, xx)[1], sf, rtol=1e-12)
+    @test isapprox(besselj(Float32(v), Float32(xx)), Float32(sf))
 end
 
 ## test large orders
@@ -126,7 +126,7 @@ for v in nu, xx in x
     xx *= v
     sf = SpecialFunctions.besselj(v, xx)
     @test isapprox(besselj(v, xx), sf, rtol=5e-11)
-    @test isapprox(Bessels.besseljy_positive_args(v, xx)[1], Float64(sf), rtol=5e-11)
+    @test isapprox(Bessels.BesselFunctions.besseljy_positive_args(v, xx)[1], Float64(sf), rtol=5e-11)
 end
 
 # test nu_range
@@ -143,13 +143,13 @@ end
 
 ## test large arguments
 @test isapprox(besselj(10.0, 150.0), SpecialFunctions.besselj(10.0, 150.0), rtol=1e-12)
-@test isapprox(Bessels.besseljy_large_argument(15.0, 100.0)[1], SpecialFunctions.besselj(15.0, 100.0), rtol=1e-12)
-@test isapprox(Bessels.besseljy_large_argument(15.0, 45.0)[1], SpecialFunctions.besselj(15.0, 45.0), rtol=1e-12)
-@test isapprox(Bessels.besseljy_large_argument(15.0, 25.5)[1], SpecialFunctions.besselj(15.0, 25.5), rtol=1e-12)
+@test isapprox(Bessels.BesselFunctions.besseljy_large_argument(15.0, 100.0)[1], SpecialFunctions.besselj(15.0, 100.0), rtol=1e-12)
+@test isapprox(Bessels.BesselFunctions.besseljy_large_argument(15.0, 45.0)[1], SpecialFunctions.besselj(15.0, 45.0), rtol=1e-12)
+@test isapprox(Bessels.BesselFunctions.besseljy_large_argument(15.0, 25.5)[1], SpecialFunctions.besselj(15.0, 25.5), rtol=1e-12)
 
 # test BigFloat for single point
-@test isapprox(Bessels.besseljy_debye(big"2000", big"1500.0")[1], SpecialFunctions.besselj(big"2000", big"1500"), rtol=5e-20)
-@test isapprox(Bessels.besseljy_large_argument(big"20", big"1500.0")[1], SpecialFunctions.besselj(big"20", big"1500"), rtol=5e-20)
+@test isapprox(Bessels.BesselFunctions.besseljy_debye(big"2000", big"1500.0")[1], SpecialFunctions.besselj(big"2000", big"1500"), rtol=5e-20)
+@test isapprox(Bessels.BesselFunctions.besseljy_large_argument(big"20", big"1500.0")[1], SpecialFunctions.besselj(big"20", big"1500"), rtol=5e-20)
 
 # need to test accuracy of negative orders and negative arguments and all combinations within
 # SpecialFunctions.jl doesn't provide these so will hand check against hard values

test/bessely_test.jl

@@ -82,7 +82,7 @@ for v in nu, xx in x
     xx *= v
     sf = SpecialFunctions.bessely(BigFloat(v), BigFloat(xx))
     @test isapprox(bessely(v, xx), Float64(sf), rtol=2e-13)
-    @test isapprox(Bessels.besseljy_positive_args(v, xx)[2], Float64(sf), rtol=5e-12)
+    @test isapprox(Bessels.BesselFunctions.besseljy_positive_args(v, xx)[2], Float64(sf), rtol=5e-12)
     @test isapprox(bessely(Float32(v), Float32(xx)), Float32(sf))
 end
 
@@ -95,7 +95,7 @@ for v in nu, xx in x
     xx *= v
     sf = SpecialFunctions.bessely(v, xx)
     @test isapprox(bessely(v, xx), sf, rtol=5e-12)
-    @test isapprox(Bessels.besseljy_positive_args(v, xx)[2], SpecialFunctions.bessely(v, xx), rtol=5e-12)
+    @test isapprox(Bessels.BesselFunctions.besseljy_positive_args(v, xx)[2], SpecialFunctions.bessely(v, xx), rtol=5e-12)
     @test isapprox(bessely(Float32(v), Float32(xx)), Float32(sf))
 end
 

test/gamma_test.jl

@@ -1,16 +1,16 @@
 for (T, max, rtol) in ((Float16, 13, 1.0), (Float32, 43, 1.0), (Float64, 170, 7))
-    v = rand(T, 10000)*max
+    local v = rand(T, 10000)*max
     for x in v
-        @test isapprox(T(SpecialFunctions.gamma(widen(x))), Bessels.gamma(x), rtol=rtol*eps(T))
+        @test isapprox(T(SpecialFunctions.gamma(widen(x))), gamma(x), rtol=rtol*eps(T))
         if isinteger(x) && x != 0
-            @test_throws DomainError Bessels.gamma(-x)
+            @test_throws DomainError gamma(-x)
         else
-            @test isapprox(T(SpecialFunctions.gamma(widen(-x))), Bessels.gamma(-x), atol=nextfloat(T(0.),2), rtol=rtol*eps(T))
+            @test isapprox(T(SpecialFunctions.gamma(widen(-x))), gamma(-x), atol=nextfloat(T(0.),2), rtol=rtol*eps(T))
         end
     end
-    @test isnan(Bessels.gamma(T(NaN)))
-    @test isinf(Bessels.gamma(T(Inf)))
+    @test isnan(gamma(T(NaN)))
+    @test isinf(gamma(T(Inf)))
 end
 
 x = [0, 1, 2, 3, 8, 15, 20, 30]
-@test SpecialFunctions.gamma.(x) ≈ Bessels.gamma.(x)
+@test SpecialFunctions.gamma.(x) ≈ gamma.(x)