fix inf and zero cases

Author: heltonmc

src/airy.jl

@@ -24,7 +24,7 @@ function _airyai(x::T) where T <: Union{Float32, Float64}
     z = 2 * x_abs^(T(3)/2) / 3
 
     if x > zero(T)
-        return sqrt(x / 3) * besselk(one(T)/3, z) / π
+        return isinf(z) ? zero(T) : sqrt(x / 3) * besselk(one(T)/3, z) / π
     elseif x < zero(T)
         Jv, Yv = besseljy_positive_args(one(T)/3, z)
         Jmv = (Jv - sqrt(T(3)) * Yv) / 2
@@ -48,7 +48,7 @@ function _airyaiprime(x::T) where T <: Union{Float32, Float64}
     z = 2 * x_abs^(T(3)/2) / 3
 
     if x > zero(T)
-        return -x * besselk(T(2)/3, z) / (sqrt(T(3)) * π)
+        return isinf(z) ? zero(T) : -x * besselk(T(2)/3, z) / (sqrt(T(3)) * π)
     elseif x < zero(T)
         Jv, Yv = besseljy_positive_args(T(2)/3, z)
         Jmv = -(Jv + sqrt(T(3))*Yv) / 2
@@ -72,7 +72,7 @@ function _airybi(x::T) where T <: Union{Float32, Float64}
     z = 2 * x_abs^(T(3)/2) / 3
 
     if x > zero(T)
-        return sqrt(x / 3) * (besseli(-one(T)/3, z) + besseli(one(T)/3, z))
+        return isinf(z) ? z : sqrt(x / 3) * (besseli(-one(T)/3, z) + besseli(one(T)/3, z))
     elseif x < zero(T)
         Jv, Yv = besseljy_positive_args(one(T)/3, z)
         Jmv = (Jv - sqrt(T(3)) * Yv) / 2
@@ -86,13 +86,17 @@ end
     airybiprime(x)
 Derivative of the Airy function of the second kind ``\\operatorname{Bi}'(x)``.
 """
-function airybiprime(x::T) where T
+airybiprime(x::Real) = _airybiprime(float(x))
+
+_airybiprime(x::Float16) = Float16(_airybiprime(Float32(x)))
+
+function _airybiprime(x::T) where T <: Union{Float32, Float64}
     isnan(x) && return x
     x_abs = abs(x)
     z = 2 * x_abs^(T(3)/2) / 3
 
     if x > zero(T)
-        return x * (besseli(T(2)/3, z) + besseli(-T(2)/3, z)) / sqrt(T(3))
+        return isinf(z) ? z : x * (besseli(T(2)/3, z) + besseli(-T(2)/3, z)) / sqrt(T(3))
     elseif x < zero(T)
         Jv, Yv = besseljy_positive_args(T(2)/3, z)
         Jmv = -(Jv + sqrt(T(3))*Yv) / 2

src/besseli.jl

@@ -207,6 +207,7 @@ Modified Bessel function of the first kind of order nu, ``I_{nu}(x)`` for positi
 function besseli_positive_args(nu, x::T) where T <: Union{Float32, Float64}
     iszero(nu) && return besseli0(x)
     isone(nu) && return besseli1(x)
+    isinf(x) && return T(Inf)
 
     # use large argument expansion if x >> nu
     besseli_large_argument_cutoff(nu, x) && return besseli_large_argument(nu, x)

src/besselk.jl

@@ -217,6 +217,7 @@ Modified Bessel function of the second kind of order nu, ``K_{nu}(x)`` valid for
 """
 function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64}
     iszero(x) && return T(Inf)
+    isinf(x) && return zero(T)
 
     # dispatch to avoid uniform expansion when nu = 0 
     iszero(nu) && return besselk0(x)

test/airy_test.jl

@@ -10,15 +10,15 @@ for x in [0.0; rand(10000)*100]
     @test isapprox(airyaiprime(-x), SpecialFunctions.airyaiprime(-x), rtol=1e-9)
 
     @test isapprox(airybi(x), SpecialFunctions.airybi(x), rtol=1e-12)
-    @test isapprox(airybi(-x), SpecialFunctions.airybi(-x), rtol=1e-9)
+    @test isapprox(airybi(-x), SpecialFunctions.airybi(-x), rtol=5e-8)
 
     @test isapprox(airybiprime(x), SpecialFunctions.airybiprime(x), rtol=1e-12)
     @test isapprox(airybiprime(-x), SpecialFunctions.airybiprime(-x), rtol=1e-9)
 end
 
 # test Inf
 @test iszero(airyai(Inf))
-@test iszer(airyaiprime(Inf))
+@test iszero(airyaiprime(Inf))
 @test isinf(airybi(Inf))
 @test isinf(airybiprime(Inf))
 

test/runtests.jl

@@ -2,10 +2,10 @@ using Bessels
 using Test
 import SpecialFunctions
 
-#@time @testset "besseli" begin include("besseli_test.jl") end
-#@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 "hankel" begin include("hankel_test.jl") end
-#@time @testset "gamma" begin include("gamma_test.jl") end
+@time @testset "besseli" begin include("besseli_test.jl") end
+@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 "hankel" begin include("hankel_test.jl") end
+@time @testset "gamma" begin include("gamma_test.jl") end
 @time @testset "airy" begin include("airy_test.jl") end