Merge pull request #42 from JuliaMath/improve_compilation

Reduce time to first bessely

Author: heltonmc

src/Bessels.jl

@@ -57,4 +57,6 @@ include("Polynomials/besselj_polys.jl")
 include("asymptotics.jl")
 include("gamma.jl")
 
+precompile(besselj, (Float64, Float64))
+
 end

src/asymptotics.jl

@@ -200,9 +200,11 @@ nu and x must be real where nu and x can be positive or negative.
 """
 besselj(nu::Real, x::Real) = _besselj(nu, float(x))
 
+_besselj(nu, x::Float32) = Float32(_besselj(nu, Float64(x)))
+
 _besselj(nu, x::Float16) = Float16(_besselj(nu, Float32(x)))
 
-function _besselj(nu, x::T) where T <: Union{Float32, Float64}
+function _besselj(nu, x::T) where T <: Float64
     isinteger(nu) && return _besselj(Int(nu), x)
     abs_nu = abs(nu)
     abs_x = abs(x)
@@ -228,7 +230,7 @@ function _besselj(nu, x::T) where T <: Union{Float32, Float64}
     end
 end
 
-function _besselj(nu::Integer, x::T) where T <: Union{Float32, Float64}
+function _besselj(nu::Integer, x::T) where T <: Float64
     abs_nu = abs(nu)
     abs_x = abs(x)
     sg = iseven(abs_nu) ? 1 : -1
@@ -318,7 +320,7 @@ function besselj_power_series(v, x::T) where T
 end
 
 besselj_series_cutoff(v, x::Float64) = (x < 7.0) || v > (2 + x*(0.109 + 0.062x))
-besselj_series_cutoff(v, x::Float32) = (x < 20.0) || v > (14.4 + x*(-0.455 + 0.027x))
+#besselj_series_cutoff(v, x::Float32) = (x < 20.0) || v > (14.4 + x*(-0.455 + 0.027x))
 
 #=
 # this needs a better way to sum these as it produces large errors

src/besselj.jl

@@ -237,9 +237,11 @@ nu and x must be real where nu and x can be positive or negative.
 """
 bessely(nu::Real, x::Real) = _bessely(nu, float(x))
 
+_bessely(nu, x::Float32) = Float32(_bessely(nu, Float64(x)))
+
 _bessely(nu, x::Float16) = Float16(_bessely(nu, Float32(x)))
 
-function _bessely(nu, x::T) where T <: Union{Float32, Float64}
+function _bessely(nu, x::T) where T <: Float64
     isnan(nu) || isnan(x) && return NaN
     isinteger(nu) && return _bessely(Int(nu), x)
     abs_nu = abs(nu)
@@ -264,7 +266,8 @@ function _bessely(nu, x::T) where T <: Union{Float32, Float64}
         end
     end
 end
-function _bessely(nu::Integer, x::T) where T <: Union{Float32, Float64}
+
+function _bessely(nu::Integer, x::T) where T <: Float64
     abs_nu = abs(nu)
     abs_x = abs(x)
     sg = iseven(abs_nu) ? 1 : -1

src/bessely.jl

@@ -24,9 +24,11 @@ solution to the radial part of the Helmholz equation in spherical coordinates.
 """
 sphericalbesselj(nu::Real, x::Real) = _sphericalbesselj(nu, float(x))
 
+_sphericalbesselj(nu, x::Float32) = Float32(_sphericalbesselj(nu, Float64(x)))
+
 _sphericalbesselj(nu, x::Float16) = Float16(_sphericalbesselj(nu, Float32(x)))
 
-function _sphericalbesselj(nu::Real, x::T) where T <: Union{Float32, Float64}
+function _sphericalbesselj(nu::Real, x::T) where T <: Float64
     x < zero(T) && return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
     if ~isfinite(x)
         isnan(x) && return x
@@ -110,9 +112,11 @@ known as a spherical Neumann function.
 """
 sphericalbessely(nu::Real, x::Real) = _sphericalbessely(nu, float(x))
 
+_sphericalbessely(nu, x::Float32) = Float32(_sphericalbessely(nu, Float64(x)))
+
 _sphericalbessely(nu, x::Float16) = Float16(_sphericalbessely(nu, Float32(x)))
 
-function _sphericalbessely(nu::Real, x::T) where T <: Union{Float32, Float64}
+function _sphericalbessely(nu::Real, x::T) where T <: Float64
     x < zero(T) && return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
     if ~isfinite(x)
         isnan(x) && return x

src/sphericalbessel.jl

@@ -110,17 +110,20 @@ for v in nu, xx in x
     @test isapprox(Bessels.besseljy_positive_args(v, xx)[1], Float64(sf), rtol=5e-11)
 end
 
-# test Float16
+# test Float16 and Float32
 @test besselj(10, Float16(1.0)) isa Float16
+@test besselj(10.2f0, 1.0f0) isa Float32
 
 ## 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 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)
 
-
 # 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
 # values taken from https://keisan.casio.com/exec/system/1180573474 which match mathematica

test/besselj_test.jl

@@ -91,6 +91,7 @@ end
 
 # test Float16
 @test bessely(10, Float16(1.0)) isa Float16
+@test bessely(10.2f0, 1.0f0) isa Float32
 
 # test limits for small arguments see https://github.com/JuliaMath/Bessels.jl/issues/35
 @test bessely(185.0, 1.01) == -Inf

test/bessely_test.jl

@@ -32,9 +32,12 @@ x = 1e-15
 @test isnan(Bessels.sphericalbessely(1.4, NaN))
 @test isnan(Bessels.sphericalbessely(4.0, NaN))
 
-# test Float16 types
+# test Float16, Float32 types
 @test Bessels.sphericalbesselj(Float16(1.4), Float16(1.2)) isa Float16
 @test Bessels.sphericalbessely(Float16(1.4), Float16(1.2)) isa Float16
+@test Bessels.sphericalbesselj(1.4f0, 1.2f0) isa Float32
+@test Bessels.sphericalbessely(1.4f0, 1.2f0) isa Float32
+
 
 for x in 0.5:1.0:100.0, v in [0, 1, 5.5, 8.2, 10]
     @test isapprox(Bessels.sphericalbesselj(v, x), SpecialFunctions.sphericalbesselj(v, x), rtol=1e-12)