Merge pull request #77 from heltonmc/airy

Improved Airy function implementations for real arguments

Project.toml

@@ -1,7 +1,6 @@
 name = "Bessels"
 uuid = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
-authors = ["Michael Helton <[email protected]> and contributors"]
-version = "0.2.8"
+version = "0.3-DEV"
 
 [compat]
 julia = "1.6"

src/airy.jl → src/Airy/cairy.jl

@@ -1,7 +1,7 @@
 #                           Airy functions
 #
-#                       airyai(z), airybi(nu, z)
-#                   airyaiprime(z), airybiprime(nu, z)
+#                       airyai(z), airybi(z)
+#                   airyaiprime(z), airybiprime(z)
 #
 #    A numerical routine to compute the airy functions and their derivatives in the entire complex plane.
 #    These routines are based on the methods reported in [1] which use a combination of the power series
@@ -48,10 +48,9 @@ See also: [`airyaiprime`](@ref), [`airybi`](@ref)
 """
 airyai(z::Number) = _airyai(float(z))
 
-_airyai(x::Float16) = Float16(_airyai(Float32(x)))
 _airyai(z::ComplexF16) = ComplexF16(_airyai(ComplexF32(z)))
 
-function _airyai(z::ComplexOrReal{T}) where T <: Union{Float32, Float64}
+function _airyai(z::Complex{T}) where T <: Union{Float32, Float64}
     if ~isfinite(z)
         if abs(angle(z)) < 2*T(π)/3
             return exp(-z)
@@ -87,7 +86,7 @@ end
 """
     airyaiprime(z)
 
-Returns the derivative of the Airy function of the first kind, ``\\operatorname{Ai}'(z)``.
+Returns the derivative of the Airy function of the first kind, 1`\\operatorname{Ai}'(z)``.
 Routine supports single and double precision (e.g., `Float32`,  `Float64`, `ComplexF64`) for real and complex arguments.
 
 # Examples
@@ -106,7 +105,6 @@ See also: [`airyai`](@ref), [`airybi`](@ref)
 """
 airyaiprime(z::Number) = _airyaiprime(float(z))
 
-_airyaiprime(x::Float16) = Float16(_airyaiprime(Float32(x)))
 _airyaiprime(z::ComplexF16) = ComplexF16(_airyaiprime(ComplexF32(z)))
 
 function _airyaiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64}
@@ -164,7 +162,6 @@ See also: [`airybiprime`](@ref), [`airyai`](@ref)
 """
 airybi(z::Number) = _airybi(float(z))
 
-_airybi(x::Float16) = Float16(_airybi(Float32(x)))
 _airybi(z::ComplexF16) = ComplexF16(_airybi(ComplexF32(z)))
 
 function _airybi(z::ComplexOrReal{T}) where T <: Union{Float32, Float64}
@@ -224,7 +221,6 @@ See also: [`airybi`](@ref), [`airyai`](@ref)
 """
 airybiprime(z::Number) = _airybiprime(float(z))
 
-_airybiprime(x::Float16) = Float16(_airybiprime(Float32(x)))
 _airybiprime(z::ComplexF16) = ComplexF16(_airybiprime(ComplexF32(z)))
 
 function _airybiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64}

src/Bessels.jl

@@ -30,9 +30,13 @@ export hankelh1
 export hankelh2
 
 export airyai
+export airyaix
 export airyaiprime
+export airyaiprimex
 export airybi
+export airybix
 export airybiprime
+export airybiprimex
 
 const ComplexOrReal{T} = Union{T,Complex{T}}
 
@@ -41,7 +45,8 @@ include("besselj.jl")
 include("besselk.jl")
 include("bessely.jl")
 include("hankel.jl")
-include("airy.jl")
+include("Airy/airy.jl")
+include("Airy/cairy.jl")
 include("sphericalbessel.jl")
 include("modifiedsphericalbessel.jl")
 

src/math_constants.jl

@@ -16,6 +16,8 @@ const SQPIO2(::Type{Float64}) = 1.25331413731550025
 const SQ1O2PI(::Type{Float64}) = 0.3989422804014327
 const SQ2PI(::Type{Float64}) = 2.5066282746310007
 const SQ2O2(::Type{Float64}) = 0.7071067811865476
+const PIPOW3O2(::Type{Float64}) = 5.568327996831708 # pi^(3/2)
+const ONEOSQPI(::Type{Float64}) = 0.5641895835477563
 
 const PIO4(::Type{Float32}) = 0.78539816339744830962f0
 const TWOOPI(::Type{Float32}) = 0.636619772367581343075535f0

test/Manifest.toml

@@ -2,6 +2,7 @@
 
 [[ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+version = "1.1.1"
 
 [[Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
@@ -30,6 +31,7 @@ version = "3.41.0"
 [[CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+version = "0.5.2+0"
 
 [[Dates]]
 deps = ["Printf"]
@@ -50,8 +52,12 @@ uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 version = "0.8.6"
 
 [[Downloads]]
-deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
 uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+version = "1.6.0"
+
+[[FileWatching]]
+uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
 [[InteractiveUtils]]
 deps = ["Markdown"]
@@ -77,10 +83,12 @@ version = "1.3.0"
 [[LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+version = "0.6.3"
 
 [[LibCURL_jll]]
 deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
 uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+version = "7.84.0+0"
 
 [[LibGit2]]
 deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
@@ -89,6 +97,7 @@ uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 [[LibSSH2_jll]]
 deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
 uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+version = "1.10.2+0"
 
 [[Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -113,23 +122,28 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 [[MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+version = "2.28.0+0"
 
 [[Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
 [[MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+version = "2022.2.1"
 
 [[NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+version = "1.2.0"
 
 [[OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.20+0"
 
 [[OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
+version = "0.8.1+0"
 
 [[OpenSpecFun_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
@@ -140,6 +154,7 @@ version = "0.5.5+0"
 [[Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+version = "1.8.0"
 
 [[Preferences]]
 deps = ["TOML"]
@@ -161,6 +176,7 @@ uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+version = "0.7.0"
 
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -189,10 +205,12 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 [[TOML]]
 deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+version = "1.0.0"
 
 [[Tar]]
 deps = ["ArgTools", "SHA"]
 uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+version = "1.10.1"
 
 [[Test]]
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
@@ -208,15 +226,19 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 [[Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+version = "1.2.12+3"
 
 [[libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+version = "5.1.1+0"
 
 [[nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+version = "1.48.0+0"
 
 [[p7zip_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
+version = "17.4.0+0"

test/airy_test.jl

@@ -1,34 +1,51 @@
 # test the airy functions
 # these are prone to some cancellation as they are indirectly calculated from relations to bessel functions
 # this is amplified for negative arguments
-# the implementations provided by SpecialFunctions.jl suffer from same inaccuracies
-for x in [0.0; 1e-17:0.1:100.0]
-    @test isapprox(airyai(x), SpecialFunctions.airyai(x), rtol=2e-13)
-    @test isapprox(airyai(-x), SpecialFunctions.airyai(-x), rtol=3e-12)
 
-    @test isapprox(airyaiprime(x), SpecialFunctions.airyaiprime(x), rtol=2e-13)
-    @test isapprox(airyaiprime(-x), SpecialFunctions.airyaiprime(-x), rtol=5e-12)
-
-    @test isapprox(airybi(x), SpecialFunctions.airybi(x), rtol=2e-13)
-    @test isapprox(airybi(-x), SpecialFunctions.airybi(-x), rtol=5e-12)
-
-    @test isapprox(airybiprime(x), SpecialFunctions.airybiprime(x), rtol=2e-13)
-    @test isapprox(airybiprime(-x), SpecialFunctions.airybiprime(-x), rtol=5e-12)
+# test Float64 for positive arguments
+for line in eachline("data/airy/airy_positive_args.csv")
+    T = Float64
+    local x, aix, aiprimex, bix, biprimex
+    x, aix, aiprimex, bix, biprimex = parse.(Float64, split(line))
+    @test isapprox(airyaix(x), aix, rtol=4.4e-16)
+    @test isapprox(airybix(x), bix, rtol=5.5e-16)
+    @test isapprox(airyaiprimex(x), aiprimex, rtol=4.9e-16)
+    @test isapprox(airybiprimex(x), biprimex, rtol=5.9e-16)
+    if x < 105.0
+        x_big = BigFloat(x)
+        scale = exp(-2 * sqrt(x_big) * x_big / 3)
+        if x < 2.0
+            tol = 4e-16
+        else
+            tol = 2.4e-16 * sqrt(x) * x
+        end
+        @test isapprox(airyai(x), T(aix * scale), rtol=tol)
+        @test isapprox(airybi(x), T(bix / scale), rtol=tol)
+        @test isapprox(airyaiprime(x), T(aiprimex * scale), rtol=tol)
+        @test isapprox(airybiprime(x), T(biprimex / scale), rtol=tol)
+    end
 end
 
-# Float32
-for x in [0.0; 0.5:0.5:30.0]
-    @test isapprox(airyai(x), SpecialFunctions.airyai(x), rtol=2e-13)
-    @test isapprox(airyai(-x), SpecialFunctions.airyai(-x), rtol=3e-12)
-
-    @test isapprox(airyaiprime(x), SpecialFunctions.airyaiprime(x), rtol=2e-13)
-    @test isapprox(airyaiprime(-x), SpecialFunctions.airyaiprime(-x), rtol=5e-12)
-
-    @test isapprox(airybi(x), SpecialFunctions.airybi(x), rtol=2e-13)
-    @test isapprox(airybi(-x), SpecialFunctions.airybi(-x), rtol=5e-12)
-
-    @test isapprox(airybiprime(x), SpecialFunctions.airybiprime(x), rtol=2e-13)
-    @test isapprox(airybiprime(-x), SpecialFunctions.airybiprime(-x), rtol=5e-12)
+# test Float64 for negative arguments
+for line in eachline("data/airy/airy_negative_args.csv")
+    local x, aix, aiprimex, bix, biprimex
+    x, ai, aiprime, bi, biprime = parse.(Float64, split(line))
+    if x >= -9.5
+        tol = 2.4e-16
+        @test isapprox(airyai(x), ai, rtol=tol)
+        @test isapprox(airybi(x), bi, rtol=tol)
+        @test isapprox(airyaiprime(x), aiprime, rtol=tol)
+        @test isapprox(airybiprime(x), biprime, rtol=tol)
+    elseif x >= -1e8
+        tol = 0.8e-16 * abs(x)^(5/4)
+        tol2 = 0.8e-16 * abs(x)^(7/4)
+        @test isapprox(airyai(x), ai, atol=tol)
+        @test isapprox(airybi(x), bi, atol=tol)
+        @test isapprox(airybix(x), bi, atol=tol)
+        @test isapprox(airyaiprime(x), aiprime, atol=tol2)
+        @test isapprox(airybiprime(x), biprime, atol=tol2)
+        @test isapprox(airybiprimex(x), biprime, atol=tol2)
+    end
 end
 
 for x in [0.0, 0.01, 0.5, 1.0, 2.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 50.0], a in 0:pi/12:2pi
@@ -39,12 +56,26 @@ for x in [0.0, 0.01, 0.5, 1.0, 2.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 50.0], a
     @test isapprox(airybiprime(z), SpecialFunctions.airybiprime(z), rtol=1e-11)
 end
 
+
 # test Inf
+# results match mathematica using Limit[AiryAi[x], x -> Infinity]
+# for scaled arguments Limit[AiryAi[x] * Exp[2 * Sqrt[x] * x / 3], x -> Infinity]
+@test airyai(Inf) === 0.0
+@test airyaiprime(Inf) === 0.0
+@test airybi(Inf) === Inf
+@test airybiprime(Inf) === Inf
 
-@test iszero(airyai(Inf))
-@test iszero(airyaiprime(Inf))
-@test isinf(airybi(Inf))
-@test isinf(airybiprime(Inf))
+@test airyaix(Inf) === 0.0
+@test airyaiprimex(Inf) === -Inf
+@test airybix(Inf) === 0.0
+@test airybiprimex(Inf) === Inf
+
+# test negative infinite arguments
+@test_throws DomainError airyai(-Inf)
+# @test airyaiprime(-Inf) === Nan # value is indeterminate
+@test_throws DomainError airybi(-Inf)
+@test_throws DomainError airybix(-Inf)
+# @test airybiprime(-Inf) === Nan # value is indeterminate
 
 @test airyai(Inf + 0.0im) === exp(-(Inf + 0.0im))
 @test airyaiprime(Inf + 0.0im) === -exp(-(Inf + 0.0im))
@@ -55,13 +86,17 @@ end
 @test airybiprime(Inf + 0.0*im) === exp((Inf + 0.0*im))
 @test airybiprime(-Inf + 0.0*im) === -1 / (-Inf + 0.0*im)
 
+# test NaNs
+for f in (:airyai, :airyaix, :airyaiprime, :airyaiprimex, :airybi, :airybix, :airybiprime, :airybiprimex)
+    @test @eval isnan($f(NaN))
+end
 
-# test Float16 types
-@test airyai(Float16(1.2)) isa Float16
-@test airyai(ComplexF16(1.2)) isa ComplexF16
-@test airyaiprime(Float16(1.9)) isa Float16
-@test airyaiprime(ComplexF16(1.2)) isa ComplexF16
-@test airybi(Float16(1.2)) isa Float16
-@test airybi(ComplexF16(1.2)) isa ComplexF16
-@test airybiprime(Float16(1.9)) isa Float16
-@test airybiprime(ComplexF16(1.2)) isa ComplexF16
+# test Float16 and Float32 types
+for f in (:airyai, :airyaix, :airyaiprime, :airyaiprimex, :airybi, :airybix, :airybiprime, :airybiprimex), T in (:Float16, :Float32)
+    @test @eval $f($T(1.2)) isa $T
+end
+
+# test ComplexF16 and ComplexF32 types
+for f in (:airyai, :airyaiprime, :airybi, :airybiprime), T in (:ComplexF16, :ComplexF32)
+    @test @eval $f($T(1.2 + 1.1im)) isa $T
+end