From 33abbaf3c952bc16d144e6cb485b1f55ea1dce78 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 21 Sep 2022 11:56:26 -0400 Subject: [PATCH 01/14] init --- src/airy.jl | 101 ++++++++++++++++++++++++++++++++++++++++++++++ src/besseli.jl | 41 ++++++++++++++++++- src/gamma.jl | 106 ++++++++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 242 insertions(+), 6 deletions(-) diff --git a/src/airy.jl b/src/airy.jl index 1ed969f..494bad6 100644 --- a/src/airy.jl +++ b/src/airy.jl @@ -105,3 +105,104 @@ function _airybiprime(x::T) where T <: Union{Float32, Float64} return T(0.4482883573538264) end end + +GAMMA_TWO_THIRDS(::Type{Float64}) = 1.3541179394264005 +GAMMA_ONE_THIRD(::Type{Float64}) = 2.6789385347077475 +GAMMA_ONE_THIRD(::Type{T}) where T <: AbstractFloat = T(big"2.678938534707747633655692940974677644128689377957301100950428327590417610167733") +GAMMA_TWO_THIRDS(::Type{T}) where T <: AbstractFloat = T(big"1.354117939426400416945288028154513785519327266056793698394022467963782965401746") + + + +const ComplexOrReal{T} = Union{T,Complex{T}} +# returns both (airyai, airyaiprime) using the power series definition for small arguments +function airyai_prime_power_series(x::ComplexOrReal{T}) where T + S = eltype(x) + MaxIter = 3000 + ai1 = zero(S) + ai2 = zero(S) + aip1 = zero(S) + aip2 = zero(S) + x2 = x*x + x3 = x2*x + t = one(S) / GAMMA_TWO_THIRDS(T) + t2 = 3*x / GAMMA_ONE_THIRD(T) + t3 = one(S) / GAMMA_ONE_THIRD(T) + t4 = 3*x2 / (2*GAMMA_TWO_THIRDS(T)) + + for i in 0:MaxIter + ai1 += t + ai2 += t2 + aip1 += t3 + aip2 += t4 + @show ai1, ai2 + abs(t2) < 1e-10*eps(T) * abs(ai2) && break + t *= x3 * inv(9*(i + one(T))*(i + T(2)/3)) + t2 *= x3 * inv(9*(i + one(T))*(i + T(4)/3)) + t3 *= x3 * inv(9*(i + one(T))*(i + T(1)/3)) + t4 *= x3 * inv(9*(i + one(T))*(i + T(5)/3)) + end + return (ai1*3^(-T(2)/3) - ai2*3^(-T(4)/3), -aip1*3^(-T(1)/3) + aip2*3^(-T(5)/3)) +end + + + +function d(x::ComplexOrReal{T}) where T + S = eltype(x) + MaxIter = 3000 + ai = zero(S) + ai2 = zero(S) + + x2 = x*x + x3 = x2*x + t = one(S) / GAMMA_TWO_THIRDS(T) + t2 = 3*x / GAMMA_ONE_THIRD(T) + a = 1 / cbrt(T(9)) + + for i in 0:MaxIter + ai += a * muladd(a, -t2, t) + @show t, t2 + abs(t) < eps(T) * abs(ai) && break + t *= x3 * inv(9*(i + one(T))*(i + T(2)/3)) + t2 *= x3 * inv(9*(i + one(T))*(i + T(4)/3)) + end + return ai +end + +function d(x::ComplexOrReal{T}) where T + S = eltype(x) + MaxIter = 3000 + ai1 = zero(S) + ai2 = zero(S) + + x2 = x*x + x3 = x2*x + t = one(S) / GAMMA_TWO_THIRDS(T) + t2 = 3*x / GAMMA_ONE_THIRD(T) + + + for i in 0:MaxIter + ai1 += t + ai2 += t2 + abs(t) < eps(T) * abs(ai1) && break + t *= x3 * inv(9*(i + one(T))*(i + T(2)/3)) + t2 *= x3 * inv(9*(i + one(T))*(i + T(4)/3)) + end + return (ai1*3^(-T(2)/3) - ai2*3^(-T(4)/3)) +end + + +function airy_large_arg_a(x::T) where T + S = eltype(x) + MaxIter = 3000 + xsqr = sqrt(x) + + out = zero(S) + t = gamma(one(T) / 6) * gamma(T(5) / 6) / 4 + a = 4*xsqr*x + for i in 0:MaxIter + out += t + abs(t) < eps(T) * abs(out) && break + t *= -3*(i + one(T)/6) * (i + T(5)/6) / (a*(i + one(T))) + end + return out * exp(-a / 6) / (pi^(3/2) * sqrt(xsqr)) +end \ No newline at end of file diff --git a/src/besseli.jl b/src/besseli.jl index ce9c967..87b1cb2 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -164,7 +164,7 @@ end Modified Bessel function of the second kind of order nu, ``I_{nu}(x)``. """ -besseli(nu::Real, x::Real) = _besseli(nu, float(x)) +besseli(nu::Real, x) = _besseli(nu, float(x)) _besseli(nu, x::Float16) = Float16(_besseli(nu, Float32(x))) @@ -394,3 +394,42 @@ function steed(n, x::T) where T return h end =# +function _besseli(nu, x::Complex{T}) where T <: Union{Float32, Float64} + isinteger(nu) && return _besseli(Int(nu), x) + ~isfinite(x) && return x + abs_nu = abs(nu) + abs_x = abs(x) + + if nu >= 0 + if x >= 0 + return besseli_positive_args(abs_nu, abs_x) + else + return cispi(abs_nu) * besseli_positive_args(abs_nu, abs_x) + end + else + return throw(DomainError(nu, "nu must be positive if x is complex")) + #= + if x >= 0 + return besseli_positive_args(abs_nu, abs_x) + 2 / π * sinpi(abs_nu) * besselk_positive_args(abs_nu, abs_x) + else + Iv = besseli_positive_args(abs_nu, abs_x) + Kv = besselk_positive_args(abs_nu, abs_x) + return cispi(abs_nu) * Iv + 2 / π * sinpi(abs_nu) * (cispi(-abs_nu) * Kv - im * π * Iv) + end + =# + end +end +function _besseli(nu::Integer, x::Complex{T}) where T <: Union{Float32, Float64} + ~isfinite(x) && return x + abs_nu = abs(nu) + abs_x = abs(x) + sg = iseven(abs_nu) ? 1 : -1 + + if x >= 0 + return besseli_positive_args(abs_nu, abs_x) + else + return sg * besseli_positive_args(abs_nu, abs_x) + end +end + +Base.eps(::Type{Complex{T}}) where {T<:AbstractFloat} = eps(T) \ No newline at end of file diff --git a/src/gamma.jl b/src/gamma.jl index 488ad00..a0668bc 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -1,16 +1,18 @@ # Adapted from Cephes Mathematical Library (MIT license https://en.smath.com/view/CephesMathLibrary/license) by Stephen L. Moshier -function gamma(x) +gamma(z::Number) = _gamma(float(z)) +_gamma(x::Float32) = Float32(_gamma(Float64(x))) + +function _gamma(x::Float64) if x < 0 isinteger(x) && throw(DomainError(x, "NaN result for non-NaN input.")) xp1 = abs(x) + 1.0 - return π / sinpi(xp1) / _gamma(xp1) + return π / sinpi(xp1) / _gammax(xp1) else - return _gamma(x) + return _gammax(x) end end # only have a Float64 implementations -gamma(x::Float32) = Float32(gamma(Float64(x))) -function _gamma(x) +function _gammax(x) if x > 11.5 return large_gamma(x) elseif x <= 11.5 @@ -64,3 +66,97 @@ function small_gamma(x) q = evalpoly(x, Q) return z * p / q end + +_gamma(z::Complex) = exp(loggamma(z)) + +loggamma(z::Complex) = _loggamma(float(z)) + +_loggamma(x::Complex{Float32}) = ComplexF32(gamma(ComplexF64(x))) + +# copied from https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/src/gamma.jl + +# Compute the logΓ(z) function using a combination of the asymptotic series, +# the Taylor series around z=1 and z=2, the reflection formula, and the shift formula. +# Many details of these techniques are discussed in D. E. G. Hare, +# "Computing the principal branch of log-Gamma," J. Algorithms 25, pp. 221-236 (1997), +# and similar techniques are used (in a somewhat different way) by the +# SciPy loggamma function. The key identities are also described +# at http://functions.wolfram.com/GammaBetaErf/LogGamma/ +function _loggamma(z::Complex{Float64}) + x, y = reim(z) + yabs = abs(y) + if !isfinite(x) || !isfinite(y) # Inf or NaN + if isinf(x) && isfinite(y) + return Complex(x, x > 0 ? (y == 0 ? y : copysign(Inf, y)) : copysign(Inf, -y)) + elseif isfinite(x) && isinf(y) + return Complex(-Inf, y) + else + return Complex(NaN, NaN) + end + elseif x > 7 || yabs > 7 # use the Stirling asymptotic series for sufficiently large x or |y| + return loggamma_asymptotic(z) + elseif x < 0.1 # use reflection formula to transform to x > 0 + if x == 0 && y == 0 # return Inf with the correct imaginary part for z == 0 + return Complex(Inf, signbit(x) ? copysign(Float64(π), -y) : -y) + end + # the 2pi * floor(...) stuff is to choose the correct branch cut for log(sinpi(z)) + return Complex(Float64(logπ), copysign(Float64(twoπ), y) * floor(0.5*x+0.25)) - + log(sinpi(z)) - loggamma(1-z) + elseif abs(x - 1) + yabs < 0.1 + # taylor series at z=1 + # ... coefficients are [-eulergamma; [(-1)^k * zeta(k)/k for k in 2:15]] + w = Complex(x - 1, y) + return w * @evalpoly(w, -5.7721566490153286060651188e-01,8.2246703342411321823620794e-01, + -4.0068563438653142846657956e-01,2.705808084277845478790009e-01, + -2.0738555102867398526627303e-01,1.6955717699740818995241986e-01, + -1.4404989676884611811997107e-01,1.2550966952474304242233559e-01, + -1.1133426586956469049087244e-01,1.000994575127818085337147e-01, + -9.0954017145829042232609344e-02,8.3353840546109004024886499e-02, + -7.6932516411352191472827157e-02,7.1432946295361336059232779e-02, + -6.6668705882420468032903454e-02) + elseif abs(x - 2) + yabs < 0.1 + # taylor series at z=2 + # ... coefficients are [1-eulergamma; [(-1)^k * (zeta(k)-1)/k for k in 2:12]] + w = Complex(x - 2, y) + return w * @evalpoly(w, 4.2278433509846713939348812e-01,3.2246703342411321823620794e-01, + -6.7352301053198095133246196e-02,2.0580808427784547879000897e-02, + -7.3855510286739852662729527e-03,2.8905103307415232857531201e-03, + -1.1927539117032609771139825e-03,5.0966952474304242233558822e-04, + -2.2315475845357937976132853e-04,9.9457512781808533714662972e-05, + -4.4926236738133141700224489e-05,2.0507212775670691553131246e-05) + end + # use recurrence relation loggamma(z) = loggamma(z+1) - log(z) to shift to x > 7 for asymptotic series + shiftprod = Complex(x,yabs) + x += 1 + sb = false # == signbit(imag(shiftprod)) == signbit(yabs) + # To use log(product of shifts) rather than sum(logs of shifts), + # we need to keep track of the number of + to - sign flips in + # imag(shiftprod), as described in Hare (1997), proposition 2.2. + signflips = 0 + while x <= 7 + shiftprod *= Complex(x,yabs) + sb′ = signbit(imag(shiftprod)) + signflips += sb′ & (sb′ != sb) + sb = sb′ + x += 1 + end + shift = log(shiftprod) + if signbit(y) # if y is negative, conjugate the shift + shift = Complex(real(shift), signflips*-6.2831853071795864769252842 - imag(shift)) + else + shift = Complex(real(shift), imag(shift) + signflips*6.2831853071795864769252842) + end + return loggamma_asymptotic(Complex(x,y)) - shift +end + +# asymptotic series for log(gamma(z)), valid for sufficiently large real(z) or |imag(z)| +function loggamma_asymptotic(z::Complex{Float64}) + zinv = inv(z) + t = zinv*zinv + # coefficients are bernoulli[2:n+1] .// (2*(1:n).*(2*(1:n) - 1)) + return (z-0.5)*log(z) - z + 9.1893853320467274178032927e-01 + # <-- log(2pi)/2 + zinv*@evalpoly(t, 8.3333333333333333333333368e-02,-2.7777777777777777777777776e-03, + 7.9365079365079365079365075e-04,-5.9523809523809523809523806e-04, + 8.4175084175084175084175104e-04,-1.9175269175269175269175262e-03, + 6.4102564102564102564102561e-03,-2.9550653594771241830065352e-02) +end \ No newline at end of file From cf4ae3035c44086851c4c5f6ec585ef70b2ceae9 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Fri, 23 Sep 2022 16:16:57 -0400 Subject: [PATCH 02/14] add airy --- src/airy.jl | 262 ++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 225 insertions(+), 37 deletions(-) diff --git a/src/airy.jl b/src/airy.jl index 494bad6..f6d02f1 100644 --- a/src/airy.jl +++ b/src/airy.jl @@ -115,83 +115,212 @@ GAMMA_TWO_THIRDS(::Type{T}) where T <: AbstractFloat = T(big"1.35411793942640041 const ComplexOrReal{T} = Union{T,Complex{T}} # returns both (airyai, airyaiprime) using the power series definition for small arguments -function airyai_prime_power_series(x::ComplexOrReal{T}) where T + +function airyai2(z::ComplexOrReal{T}) where T + x, y = real(z), imag(z) + zabs = abs(z) + + if zabs > 8 + a = airy_large_arg_a(z) + if x < zero(T) && abs(y) < 5.5 + b = airy_large_arg_b(z) + y >= zero(T) ? (return a + im*b) : (return a - im*b) + end + return a + elseif x < 2 && abs(y) > -1.4*(x + 5.5) + return airyai_power_series(z) + elseif x > zero(T) + zz = 2 * z^(T(3)/2) / 3 + + return sqrt(z / 3) * besselk_continued_fraction_shift(one(T)/3, zz) / π + else + return exp(pi*im/3) * airyai2(-z*exp(pi*im/3)) + exp(-pi*im/3) * airyai2(-z*exp(-pi*im/3)) + end +end + +function airyaiprime2(z::ComplexOrReal{T}) where T + x, y = real(z), imag(z) + zabs = abs(z) + + if zabs > 8 + c = airy_large_arg_c(z) + if x < zero(T) && abs(y) < 5.5 + d = airy_large_arg_d(z) + y >= zero(T) ? (return c + im*d) : (return c - im*d) + end + return c + elseif x < 2 && abs(y) > -1.4*(x + 5.5) + return airyaiprime_power_series(z) + elseif x > zero(T) + zz = 2 * z^(T(3)/2) / 3 + return -z * besselk_continued_fraction_shift(T(2)/3, zz) / (π * sqrt(T(3))) + else + return -(exp(2pi*im/3)*airyaiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*airyaiprime(-z*exp(-pi*im/3))) + end +end + + +function airybi2(z::ComplexOrReal{T}) where T + x, y = real(z), imag(z) + zabs = abs(z) + if zabs > 9 + b = airy_large_arg_b(z) + if abs(y) <= 1.7*(x - 6) + return 2*b + else + a = airy_large_arg_a(z) + if x < zero(T) + if abs(y) < 5 + if y < zero(T) + return b - im*a + else + return b + im*a + end + end + end + if y < zero(T) + return 2*b - im*a + else + return 2*b + im*a + end + end + elseif x < 2 && abs(y) > -1.4*(x + 5.5) + return airybi_power_series(z) + elseif x > 2 && abs(y) > 3 + zz = 2 * z^(T(3)/2) / 3 + shift = 20 + order = one(T)/3 + inu, inum1 = besseli_power_series_inu_inum1(order + shift, zz) + inu, inum1 = besselk_down_recurrence(zz, inum1, inu, order + shift - 1, order) + + inu2, inum2 = besseli_power_series_inu_inum1(-order + shift, zz) + inu2, inum2 = besselk_down_recurrence(zz, inum2, inu2, -order + shift - 1, -order) + return sqrt(z/3) * (inu + inu2) + else + return exp(pi*im/3) * airybi2(-z*exp(pi*im/3)) + exp(-pi*im/3) * airybi2(-z*exp(-pi*im/3)) + end + +end + +function airybiprime2(z::ComplexOrReal{T}) where T + x, y = real(z), imag(z) + zabs = abs(z) + if zabs > 9 + d = airy_large_arg_d(z) + if abs(y) <= 1.7*(x - 6) + return 2*d + else + c = airy_large_arg_c(z) + if x < zero(T) + if abs(y) < 5 + if y < zero(T) + return d - im*c + else + return d + im*c + end + end + end + if y < zero(T) + return 2*d - im*c + else + return 2*d + im*c + end + end + elseif x < 0 && abs(y) < -1.4*(x + 5.5) + return -(exp(2pi*im/3)*airybiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*airybiprime(-z*exp(-pi*im/3))) + elseif x > 2 && abs(y) > 3 + zz = 2 * z^(T(3)/2) / 3 + shift = 20 + order = T(2)/3 + inu, inum1 = besseli_power_series_inu_inum1(order + shift, zz) + inu, inum1 = besselk_down_recurrence(zz, inum1, inu, order + shift - 1, order) + + inu2, inum2 = besseli_power_series_inu_inum1(-order + shift, zz) + inu2, inum2 = besselk_down_recurrence(zz, inum2, inu2, -order + shift - 1, -order) + return z / sqrt(3) * (inu + inu2) + else + return airybiprime_power_series(z) + end +end + +function airyai_power_series(x::ComplexOrReal{T}) where T S = eltype(x) MaxIter = 3000 ai1 = zero(S) ai2 = zero(S) - aip1 = zero(S) - aip2 = zero(S) x2 = x*x x3 = x2*x t = one(S) / GAMMA_TWO_THIRDS(T) t2 = 3*x / GAMMA_ONE_THIRD(T) - t3 = one(S) / GAMMA_ONE_THIRD(T) - t4 = 3*x2 / (2*GAMMA_TWO_THIRDS(T)) - + for i in 0:MaxIter ai1 += t ai2 += t2 - aip1 += t3 - aip2 += t4 - @show ai1, ai2 - abs(t2) < 1e-10*eps(T) * abs(ai2) && break + abs(t) < eps(T) * abs(ai1) && break t *= x3 * inv(9*(i + one(T))*(i + T(2)/3)) t2 *= x3 * inv(9*(i + one(T))*(i + T(4)/3)) - t3 *= x3 * inv(9*(i + one(T))*(i + T(1)/3)) - t4 *= x3 * inv(9*(i + one(T))*(i + T(5)/3)) end - return (ai1*3^(-T(2)/3) - ai2*3^(-T(4)/3), -aip1*3^(-T(1)/3) + aip2*3^(-T(5)/3)) + return (ai1*3^(-T(2)/3) - ai2*3^(-T(4)/3)) end - - - -function d(x::ComplexOrReal{T}) where T +function airyaiprime_power_series(x::ComplexOrReal{T}) where T S = eltype(x) MaxIter = 3000 - ai = zero(S) + ai1 = zero(S) ai2 = zero(S) + x2 = x*x + x3 = x2*x + t = one(S) / GAMMA_ONE_THIRD(T) + t2 = 3*x2 / (2*GAMMA_TWO_THIRDS(T)) + for i in 0:MaxIter + ai1 += t + ai2 += t2 + abs(t) < eps(T) * abs(ai1) && break + t *= x3 * inv(9*(i + one(T))*(i + T(1)/3)) + t2 *= x3 * inv(9*(i + one(T))*(i + T(5)/3)) + end + return -ai1*3^(-T(1)/3) + ai2*3^(-T(5)/3) +end +function airybi_power_series(x::ComplexOrReal{T}) where T + S = eltype(x) + MaxIter = 3000 + ai1 = zero(S) + ai2 = zero(S) x2 = x*x x3 = x2*x t = one(S) / GAMMA_TWO_THIRDS(T) t2 = 3*x / GAMMA_ONE_THIRD(T) - a = 1 / cbrt(T(9)) - + for i in 0:MaxIter - ai += a * muladd(a, -t2, t) - @show t, t2 - abs(t) < eps(T) * abs(ai) && break + ai1 += t + ai2 += t2 + abs(t) < eps(T) * abs(ai1) && break t *= x3 * inv(9*(i + one(T))*(i + T(2)/3)) t2 *= x3 * inv(9*(i + one(T))*(i + T(4)/3)) end - return ai + return (ai1*3^(-T(1)/6) + ai2*3^(-T(5)/6)) end - -function d(x::ComplexOrReal{T}) where T +function airybiprime_power_series(x::ComplexOrReal{T}) where T S = eltype(x) MaxIter = 3000 ai1 = zero(S) ai2 = zero(S) - x2 = x*x x3 = x2*x - t = one(S) / GAMMA_TWO_THIRDS(T) - t2 = 3*x / GAMMA_ONE_THIRD(T) + t = one(S) / GAMMA_ONE_THIRD(T) + t2 = 3*x2 / (2*GAMMA_TWO_THIRDS(T)) - for i in 0:MaxIter ai1 += t ai2 += t2 abs(t) < eps(T) * abs(ai1) && break - t *= x3 * inv(9*(i + one(T))*(i + T(2)/3)) - t2 *= x3 * inv(9*(i + one(T))*(i + T(4)/3)) + t *= x3 * inv(9*(i + one(T))*(i + T(1)/3)) + t2 *= x3 * inv(9*(i + one(T))*(i + T(5)/3)) end - return (ai1*3^(-T(2)/3) - ai2*3^(-T(4)/3)) + return (ai1*3^(T(1)/6) + ai2*3^(-T(7)/6)) end - -function airy_large_arg_a(x::T) where T +function airy_large_arg_a(x::ComplexOrReal{T}) where T S = eltype(x) MaxIter = 3000 xsqr = sqrt(x) @@ -201,8 +330,67 @@ function airy_large_arg_a(x::T) where T a = 4*xsqr*x for i in 0:MaxIter out += t - abs(t) < eps(T) * abs(out) && break + abs(t) < eps(T)*50 * abs(out) && break t *= -3*(i + one(T)/6) * (i + T(5)/6) / (a*(i + one(T))) end return out * exp(-a / 6) / (pi^(3/2) * sqrt(xsqr)) -end \ No newline at end of file +end + +function airy_large_arg_b(x::ComplexOrReal{T}) where T + S = eltype(x) + MaxIter = 3000 + xsqr = sqrt(x) + + out = zero(S) + t = gamma(one(T) / 6) * gamma(T(5) / 6) / 4 + a = 4*xsqr*x + for i in 0:MaxIter + out += t + abs(t) < eps(T)*50 * abs(out) && break + t *= 3*(i + one(T)/6) * (i + T(5)/6) / (a*(i + one(T))) + end + return out * exp(a / 6) / (pi^(3/2) * sqrt(xsqr)) +end + +function airy_large_arg_c(x::ComplexOrReal{T}) where T + S = eltype(x) + MaxIter = 3000 + xsqr = sqrt(x) + + out = zero(S) + t = gamma(-one(T) / 6) * gamma(T(7) / 6) / 4 + a = 4*xsqr*x + for i in 0:MaxIter + out += t + abs(t) < eps(T)*50* abs(out) && break + t *= -3*(i - one(T)/6) * (i + T(7)/6) / (a*(i + one(T))) + end + return out * exp(-a / 6) * sqrt(xsqr) / pi^(3/2) +end + +function airy_large_arg_d(x::ComplexOrReal{T}) where T + S = eltype(x) + MaxIter = 3000 + xsqr = sqrt(x) + + out = zero(S) + t = gamma(-one(T) / 6) * gamma(T(7) / 6) / 4 + a = 4*xsqr*x + for i in 0:MaxIter + out += t + abs(t) < eps(T)*50 * abs(out) && break + t *= 3*(i - one(T)/6) * (i + T(7)/6) / (a*(i + one(T))) + end + return -out * exp(a / 6) * sqrt(xsqr) / pi^(3/2) +end + +# calculates besselk from the power series of besseli using the continued fraction and wronskian +# this shift the order higher first to avoid cancellation in the power series of besseli along the imaginary axis +# for real arguments this is not needed because besseli can be computed stably over the entire real axis +function besselk_continued_fraction_shift(nu, x) + shift = 20 + inu, inum1 = besseli_power_series_inu_inum1(nu + shift, x) + inu, inum1 = besselk_down_recurrence(x, inum1, inu, nu + shift - 1, nu) + H_knu = besselk_ratio_knu_knup1(nu-1, x) + return 1 / (x * (inum1 + inu / H_knu)) +end From ac98ea4d96eff259fffdfe5ad9bce46bd64751cc Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Sat, 24 Sep 2022 16:35:59 -0400 Subject: [PATCH 03/14] add complex support --- src/airy.jl | 360 ++++++++++++++++++++++++---------------------- src/besseli.jl | 9 +- src/besselk.jl | 27 ++-- src/recurrence.jl | 3 - 4 files changed, 209 insertions(+), 190 deletions(-) diff --git a/src/airy.jl b/src/airy.jl index f6d02f1..ac816f9 100644 --- a/src/airy.jl +++ b/src/airy.jl @@ -10,101 +10,7 @@ # In the future, these could be replaced with more custom routines as they depend on a single variable. # -""" - airyai(x) -Airy function of the first kind ``\\operatorname{Ai}(x)``. -""" -airyai(x::Real) = _airyai(float(x)) -_airyai(x::Float16) = Float16(_airyai(Float32(x))) - -function _airyai(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 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 - return isinf(z) ? throw(DomainError(x, "airyai(x) is not defined.")) : sqrt(x_abs) * (Jmv + Jv) / 3 - elseif iszero(x) - return T(0.3550280538878172) - end -end - -""" - airyaiprime(x) -Derivative of the Airy function of the first kind ``\\operatorname{Ai}'(x)``. -""" -airyaiprime(x::Real) = _airyaiprime(float(x)) - -_airyaiprime(x::Float16) = Float16(_airyaiprime(Float32(x))) - -function _airyaiprime(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 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 - return isinf(z) ? throw(DomainError(x, "airyaiprime(x) is not defined.")) : x_abs * (Jv - Jmv) / 3 - elseif iszero(x) - return T(-0.2588194037928068) - end -end - -""" - airybi(x) -Airy function of the second kind ``\\operatorname{Bi}(x)``. -""" -airybi(x::Real) = _airybi(float(x)) - -_airybi(x::Float16) = Float16(_airybi(Float32(x))) - -function _airybi(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 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 - return isinf(z) ? throw(DomainError(x, "airybi(x) is not defined.")) : sqrt(x_abs/3) * (Jmv - Jv) - elseif iszero(x) - return T(0.6149266274460007) - end -end - -""" - airybiprime(x) -Derivative of the Airy function of the second kind ``\\operatorname{Bi}'(x)``. -""" -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 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 - return isinf(z) ? throw(DomainError(x, "airybiprime(x) is not defined.")) : x_abs * (Jv + Jmv) / sqrt(T(3)) - elseif iszero(x) - return T(0.4482883573538264) - end -end GAMMA_TWO_THIRDS(::Type{Float64}) = 1.3541179394264005 GAMMA_ONE_THIRD(::Type{Float64}) = 2.6789385347077475 @@ -116,25 +22,34 @@ GAMMA_TWO_THIRDS(::Type{T}) where T <: AbstractFloat = T(big"1.35411793942640041 const ComplexOrReal{T} = Union{T,Complex{T}} # returns both (airyai, airyaiprime) using the power series definition for small arguments +airyai(z::Number) = _airyai(float(x)) + +_airyai(x::T) where T <: ComplexOrReal{Float16} = T(_airyai(Float32(x))) + function airyai2(z::ComplexOrReal{T}) where T x, y = real(z), imag(z) zabs = abs(z) - if zabs > 8 - a = airy_large_arg_a(z) - if x < zero(T) && abs(y) < 5.5 - b = airy_large_arg_b(z) - y >= zero(T) ? (return a + im*b) : (return a - im*b) - end - return a - elseif x < 2 && abs(y) > -1.4*(x + 5.5) - return airyai_power_series(z) - elseif x > zero(T) + airy_large_argument_cutoff(z) && return airyai_large_argument(z) + airyai_power_series_cutoff(x, y) && return airyai_power_series(z) + + if x > zero(T) + # use relation to besselk (http://dlmf.nist.gov/9.6.E1) zz = 2 * z^(T(3)/2) / 3 - return sqrt(z / 3) * besselk_continued_fraction_shift(one(T)/3, zz) / π else - return exp(pi*im/3) * airyai2(-z*exp(pi*im/3)) + exp(-pi*im/3) * airyai2(-z*exp(-pi*im/3)) + # z is close to the negative real axis + # for imag(z) == 0 use reflection to compute in terms of bessel functions of first kind (http://dlmf.nist.gov/9.6.E5) + # use computation for real numbers then convert to input type for stability + # for imag(z) != 0 use rotation identity (http://dlmf.nist.gov/9.2.E14) + if iszero(y) + xx = 2 * x^(T(3)/2) / 3 + Jv, Yv = besseljy_positive_args(one(T)/3, xx) + Jmv = (Jv - sqrt(T(3)) * Yv) / 2 + return convert(eltype(z), sqrt(abs(x)) * (Jmv + Jv) / 3) + else + return exp(pi*im/3) * airyai2(-z*exp(pi*im/3)) + exp(-pi*im/3) * airyai2(-z*exp(-pi*im/3)) + end end end @@ -142,51 +57,38 @@ function airyaiprime2(z::ComplexOrReal{T}) where T x, y = real(z), imag(z) zabs = abs(z) - if zabs > 8 - c = airy_large_arg_c(z) - if x < zero(T) && abs(y) < 5.5 - d = airy_large_arg_d(z) - y >= zero(T) ? (return c + im*d) : (return c - im*d) - end - return c - elseif x < 2 && abs(y) > -1.4*(x + 5.5) - return airyaiprime_power_series(z) - elseif x > zero(T) + airy_large_argument_cutoff(z) && return airyaiprime_large_argument(z) + airyai_power_series_cutoff(x, y) && return airyaiprime_power_series(z) + + if x > zero(T) + # use relation to besselk (http://dlmf.nist.gov/9.6.E2) zz = 2 * z^(T(3)/2) / 3 return -z * besselk_continued_fraction_shift(T(2)/3, zz) / (π * sqrt(T(3))) else - return -(exp(2pi*im/3)*airyaiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*airyaiprime(-z*exp(-pi*im/3))) + # z is close to the negative real axis + # for imag(z) == 0 use reflection to compute in terms of bessel functions of first kind (http://dlmf.nist.gov/9.6.E5) + # use computation for real numbers then convert to input type for stability + # for imag(z) != 0 use rotation identity (http://dlmf.nist.gov/9.2.E14) + if iszero(y) + xx = 2 * x^(T(3)/2) / 3 + Jv, Yv = besseljy_positive_args(T(2)/3, xx) + Jmv = -(Jv + sqrt(T(3))*Yv) / 2 + return convert(eltype(z), abs(x) * (Jv - Jmv) / 3 + else + return -(exp(2pi*im/3)*airyaiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*airyaiprime(-z*exp(-pi*im/3))) + end end end - function airybi2(z::ComplexOrReal{T}) where T x, y = real(z), imag(z) zabs = abs(z) - if zabs > 9 - b = airy_large_arg_b(z) - if abs(y) <= 1.7*(x - 6) - return 2*b - else - a = airy_large_arg_a(z) - if x < zero(T) - if abs(y) < 5 - if y < zero(T) - return b - im*a - else - return b + im*a - end - end - end - if y < zero(T) - return 2*b - im*a - else - return 2*b + im*a - end - end - elseif x < 2 && abs(y) > -1.4*(x + 5.5) - return airybi_power_series(z) - elseif x > 2 && abs(y) > 3 + + airy_large_argument_cutoff(z) && return airybi_large_argument(z) + + airybi_power_series_cutoff(x, y) && return airybi_power_series(z) + + if x > zero(T) zz = 2 * z^(T(3)/2) / 3 shift = 20 order = one(T)/3 @@ -197,38 +99,25 @@ function airybi2(z::ComplexOrReal{T}) where T inu2, inum2 = besselk_down_recurrence(zz, inum2, inu2, -order + shift - 1, -order) return sqrt(z/3) * (inu + inu2) else - return exp(pi*im/3) * airybi2(-z*exp(pi*im/3)) + exp(-pi*im/3) * airybi2(-z*exp(-pi*im/3)) - end - + if iszero(y) + xx = 2 * x^(T(3)/2) / 3 + Jv, Yv = besseljy_positive_args(one(T)/3, xx) + Jmv = (Jv - sqrt(T(3)) * Yv) / 2 + return convert(eltype(z), sqrt(abs(x)/3) * (Jmv - Jv)) + else + return exp(pi*im/3) * airybi2(-z*exp(pi*im/3)) + exp(-pi*im/3) * airybi2(-z*exp(-pi*im/3)) + end end function airybiprime2(z::ComplexOrReal{T}) where T x, y = real(z), imag(z) zabs = abs(z) - if zabs > 9 - d = airy_large_arg_d(z) - if abs(y) <= 1.7*(x - 6) - return 2*d - else - c = airy_large_arg_c(z) - if x < zero(T) - if abs(y) < 5 - if y < zero(T) - return d - im*c - else - return d + im*c - end - end - end - if y < zero(T) - return 2*d - im*c - else - return 2*d + im*c - end - end - elseif x < 0 && abs(y) < -1.4*(x + 5.5) - return -(exp(2pi*im/3)*airybiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*airybiprime(-z*exp(-pi*im/3))) - elseif x > 2 && abs(y) > 3 + + airy_large_argument_cutoff(z) && return airybiprime_large_argument(z) + + airybi_power_series_cutoff(x, y) && return airybiprime_power_series(z) + + if x > zero(T) zz = 2 * z^(T(3)/2) / 3 shift = 20 order = T(2)/3 @@ -238,11 +127,25 @@ function airybiprime2(z::ComplexOrReal{T}) where T inu2, inum2 = besseli_power_series_inu_inum1(-order + shift, zz) inu2, inum2 = besselk_down_recurrence(zz, inum2, inu2, -order + shift - 1, -order) return z / sqrt(3) * (inu + inu2) - else - return airybiprime_power_series(z) - end + else + if iszero(y) + xx = 2 * x^(T(3)/2) / 3 + Jv, Yv = besseljy_positive_args(T(2)/3, xx) + Jmv = -(Jv + sqrt(T(3))*Yv) / 2 + return convert(eltype(z), abs(x) * (Jv + Jmv) / sqrt(T(3))) + else + return -(exp(2pi*im/3)*airybiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*airybiprime(-z*exp(-pi*im/3))) + end end +##### +##### Power series for airyai(x) +##### + +# cutoffs for power series valid for both airyai and airyaiprime +airyai_power_series_cutoff(x::T, y::T) where T <: Float64 = x < 2 && abs(y) > -1.4*(x + 5.5) +airyai_power_series_cutoff(x::T, y::T) where T <: Float32 = x < 4.5f0 && abs(y) > -1.4f0*(x + 9.5f0) + function airyai_power_series(x::ComplexOrReal{T}) where T S = eltype(x) MaxIter = 3000 @@ -262,6 +165,11 @@ function airyai_power_series(x::ComplexOrReal{T}) where T end return (ai1*3^(-T(2)/3) - ai2*3^(-T(4)/3)) end + +##### +##### Power series for airyaiprime(x) +##### + function airyaiprime_power_series(x::ComplexOrReal{T}) where T S = eltype(x) MaxIter = 3000 @@ -281,6 +189,18 @@ function airyaiprime_power_series(x::ComplexOrReal{T}) where T end return -ai1*3^(-T(1)/3) + ai2*3^(-T(5)/3) end + +##### +##### Power series for airybi(x) +##### + +# cutoffs for power series valid for both airybi and airybiprime +# has a more complicated validity as it works well close to positive real line and for small negative arguments also works for angle(z) ~ 2pi/3 +# the statements are somewhat complicated but we absolutely want to hit this branch when we can as the other algorithms are 10x slower +# the Float32 cutoff can be simplified because it overlaps with the large argument expansion so there isn't a need for more complicated statements +airybi_power_series_cutoff(x::T, y::T) where T <: Float64 = (abs(y) > -1.4*(x + 5.5) && abs(y) < -2.2*(x - 4)) || (x > zero(T) && abs(y) < 3) +airybi_power_series_cutoff(x::T, y::T) where T <: Float32 = abs(complex(x, y)) < 9 + function airybi_power_series(x::ComplexOrReal{T}) where T S = eltype(x) MaxIter = 3000 @@ -320,6 +240,102 @@ function airybiprime_power_series(x::ComplexOrReal{T}) where T return (ai1*3^(T(1)/6) + ai2*3^(-T(7)/6)) end +##### +##### Large argument expansion for airy functions +##### +airy_large_argument_cutoff(z::ComplexOrReal{Float64}) = abs(z) > 8 +airy_large_argument_cutoff(z::ComplexOrReal{Float32}) = abs(z) > 4 + +function airyai_large_argument(x::Real) + x < zero(x) && return real(airyai_large_argument(complex(x))) + return airy_large_arg_a(abs(x)) +end +function airyai_large_argument(z::Complex{T}) where T + x, y = real(z), imag(z) + a = airy_large_arg_a(z) + if x < zero(T) && abs(y) < 5.5 + b = airy_large_arg_b(z) + y >= zero(T) ? (return a + im*b) : (return a - im*b) + end + return a +end + +function airyaiprime_large_argument(x::Real) + x < zero(x) && return real(airyaiprime_large_argument(complex(x))) + return airy_large_arg_c(abs(x)) +end +function airyaiprime_large_argument(z::Complex{T}) where T + x, y = real(z), imag(z) + c = airy_large_arg_c(z) + if x < zero(T) && abs(y) < 5.5 + d = airy_large_arg_d(z) + y >= zero(T) ? (return c + im*d) : (return c - im*d) + end + return c +end + +function airybi_large_argument(x::Real) + if x < zero(x) + return 2*real(airy_large_arg_b(complex(x))) + else + return 2*(airy_large_arg_b(x)) + end +end + +function airybi_large_argument(z::Complex{T}) where T + x, y = real(z), imag(z) + b = airy_large_arg_b(z) + abs(y) <= 1.7*(x - 6) && return 2*b + + check_conj = false + if y < zero(T) + z = conj(z) + y = abs(y) + check_conj = true + end + + a = airy_large_arg_a(z) + if x < zero(T) && y < 5 + out = b + im*a + check_conj && (out = conj(out)) + return out + else + out = 2*b + im*a + check_conj && (out = conj(out)) + return out + end +end +function airybiprime_large_argument(x::Real) + if x < zero(x) + return 2*real(airy_large_arg_d(complex(x))) + else + return 2*(airy_large_arg_d(x)) + end +end +function airybiprime_large_argument(z::Complex{T}) where T + x, y = real(z), imag(z) + d = airy_large_arg_d(z) + abs(y) <= 1.7*(x - 6) && return 2*d + + check_conj = false + if y < zero(T) + z = conj(z) + y = abs(y) + check_conj = true + end + + c = airy_large_arg_c(z) + if x < zero(T) && y < 5 + out = d + im*c + check_conj && (out = conj(out)) + return out + else + out = 2*d + im*c + check_conj && (out = conj(out)) + return out + end +end + function airy_large_arg_a(x::ComplexOrReal{T}) where T S = eltype(x) MaxIter = 3000 diff --git a/src/besseli.jl b/src/besseli.jl index 87b1cb2..272c1c5 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -306,9 +306,9 @@ end besseli_large_argument_scaled(v, x::T) where T = T(_besseli_large_argument(v, x) / sqrt(2 * (π * x))) -function _besseli_large_argument(v, x::T) where T +function _besseli_large_argument(v, x::ComplexOrReal{T}) where T MaxIter = 5000 - S = promote_type(T, Float64) + S = eltype(x) v, x = S(v), S(x) fv2 = 4 * v^2 @@ -339,9 +339,10 @@ besseli_large_argument_cutoff(nu, x::Float32) = x > 18.0f0 && x > nu^2 / 19.5f0 Computes ``I_{nu}(x)`` using the power series for any value of nu. """ -function besseli_power_series(v, x::T) where T +function besseli_power_series(v, x::ComplexOrReal{T}) where T MaxIter = 3000 - out = zero(T) + S = eltype(x) + out = zero(S) xs = (x/2)^v a = xs / gamma(v + one(T)) t2 = (x/2)^2 diff --git a/src/besselk.jl b/src/besselk.jl index 0bfbe6b..8e87baa 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -339,10 +339,11 @@ end # a modified version of the I_{nu} power series to compute both I_{nu} and I_{nu-1} # use this along with the continued fractions for besselk -function besseli_power_series_inu_inum1(v, x::T) where T +function besseli_power_series_inu_inum1(v, x::ComplexOrReal{T}) where T MaxIter = 3000 - out = zero(T) - out2 = zero(T) + S = eltype(x) + out = zero(S) + out2 = zero(S) x2 = x / 2 xs = x2^v gmx = xs / gamma(v) @@ -361,11 +362,12 @@ end # computes K_{nu+1}/K_{nu} using continued fractions and the modified Lentz method # generally slow to converge for small x -function besselk_ratio_knu_knup1(v, x::T) where T +function besselk_ratio_knu_knup1(v, x::ComplexOrReal{T}) where T MaxIter = 1000 - (hn, Dn, Cn) = (1e-50, zero(v), 1e-50) + S = eltype(x) + (hn, Dn, Cn) = (S(1e-150), zero(S), S(1e-150)) - jf = one(T) + jf = one(S) vv = v * v for _ in 1:MaxIter an = (vv - ((2*jf - 1)^2) * T(0.25)) @@ -398,9 +400,12 @@ Computes ``K_{nu}(x)`` using the power series when nu is not an integer. In general, this is most accurate for small arguments and when nu > x. No checks are performed on nu so this is not accurate when nu is an integer. """ -function besselk_power_series(v, x::T) where T +besselk_power_series(v, x::Float32) = Float32(besselk_power_series(v, Float64(x))) +besselk_power_series(v, x::ComplexF32) = ComplexF32(besselk_power_series(v, ComplexF64(x))) + +function besselk_power_series(v, x::ComplexOrReal{T}) where T MaxIter = 1000 - S = promote_type(T, Float64) + S = eltype(x) v, x = S(v), S(x) z = x / 2 @@ -431,7 +436,7 @@ function besselk_power_series(v, x::T) where T xd2_pow *= zz fact_k *= k + one(S) end - return T(out) + return out end besselk_power_series_cutoff(nu, x::Float64) = x < 2.0 || nu > 1.6x - 1.0 besselk_power_series_cutoff(nu, x::Float32) = x < 10.0f0 || nu > 1.65f0*x - 8.0f0 @@ -452,9 +457,9 @@ end besselk_large_argument_scaled(v, x::T) where T = T(_besselk_large_argument(v, x) * sqrt(pi / 2x)) -function _besselk_large_argument(v, x::T) where T +function _besselk_large_argument(v, x::ComplexOrReal{T}) where T MaxIter = 5000 - S = promote_type(T, Float64) + S = eltype(x) v, x = S(v), S(x) fv2 = 4 * v^2 diff --git a/src/recurrence.jl b/src/recurrence.jl index c959537..4321cf8 100644 --- a/src/recurrence.jl +++ b/src/recurrence.jl @@ -41,8 +41,6 @@ end return jnup1, jnu end -#= -# currently not used # backward recurrence relation for besselk and besseli # outputs both (bessel(x, nu_end), bessel(x, nu_end-1) # x = 0.1; k0 = besseli(10,x); k1 = besseli(11,x); @@ -55,4 +53,3 @@ end end return jnup1, jnu end -=# From 77dfa1f19b89cac60eba444edff1e7e89cd9666c Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Sat, 24 Sep 2022 17:56:02 -0400 Subject: [PATCH 04/14] add complex tests fix conjugates --- src/Bessels.jl | 2 + src/airy.jl | 105 ++++++++++++++++++++++++++++++------------ src/math_constants.jl | 5 ++ test/airy_test.jl | 27 +++++++---- 4 files changed, 101 insertions(+), 38 deletions(-) diff --git a/src/Bessels.jl b/src/Bessels.jl index 736c5a5..63e7536 100644 --- a/src/Bessels.jl +++ b/src/Bessels.jl @@ -34,6 +34,8 @@ export airyaiprime export airybi export airybiprime +const ComplexOrReal{T} = Union{T,Complex{T}} + include("besseli.jl") include("besselj.jl") include("besselk.jl") diff --git a/src/airy.jl b/src/airy.jl index ac816f9..aa46e9c 100644 --- a/src/airy.jl +++ b/src/airy.jl @@ -9,24 +9,21 @@ # For negative arguments these definitions are prone to some cancellation leading to higher errors. # In the future, these could be replaced with more custom routines as they depend on a single variable. # - - - -GAMMA_TWO_THIRDS(::Type{Float64}) = 1.3541179394264005 -GAMMA_ONE_THIRD(::Type{Float64}) = 2.6789385347077475 -GAMMA_ONE_THIRD(::Type{T}) where T <: AbstractFloat = T(big"2.678938534707747633655692940974677644128689377957301100950428327590417610167733") -GAMMA_TWO_THIRDS(::Type{T}) where T <: AbstractFloat = T(big"1.354117939426400416945288028154513785519327266056793698394022467963782965401746") - - -const ComplexOrReal{T} = Union{T,Complex{T}} -# returns both (airyai, airyaiprime) using the power series definition for small arguments +""" + airyai(x) +Airy function of the first kind ``\\operatorname{Ai}(x)``. +""" +airyai(z::Number) = _airyai(float(z)) -airyai(z::Number) = _airyai(float(x)) +_airyai(x::Float16) = Float16(_airyai(Float32(x))) +_airyai(z::ComplexF16) = ComplexF16(_airyai(ComplexF32(z))) -_airyai(x::T) where T <: ComplexOrReal{Float16} = T(_airyai(Float32(x))) - -function airyai2(z::ComplexOrReal{T}) where T +function _airyai(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} + if ~isfinite(z) + isnan(z) && return z + isinf(z) && return throw(DomainError(z, "airyai(z) not defined at infinity")) + end x, y = real(z), imag(z) zabs = abs(z) @@ -43,17 +40,31 @@ function airyai2(z::ComplexOrReal{T}) where T # use computation for real numbers then convert to input type for stability # for imag(z) != 0 use rotation identity (http://dlmf.nist.gov/9.2.E14) if iszero(y) - xx = 2 * x^(T(3)/2) / 3 + xabs = abs(x) + xx = 2 * xabs^(T(3)/2) / 3 Jv, Yv = besseljy_positive_args(one(T)/3, xx) Jmv = (Jv - sqrt(T(3)) * Yv) / 2 - return convert(eltype(z), sqrt(abs(x)) * (Jmv + Jv) / 3) + return convert(eltype(z), sqrt(xabs) * (Jmv + Jv) / 3) else - return exp(pi*im/3) * airyai2(-z*exp(pi*im/3)) + exp(-pi*im/3) * airyai2(-z*exp(-pi*im/3)) + return exp(pi*im/3) * _airyai(-z*exp(pi*im/3)) + exp(-pi*im/3) * _airyai(-z*exp(-pi*im/3)) end end end -function airyaiprime2(z::ComplexOrReal{T}) where T +""" + airyaiprime(x) +Derivative of the Airy function of the first kind ``\\operatorname{Ai}'(x)``. +""" +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} + if ~isfinite(z) + isnan(z) && return z + isinf(z) && return throw(DomainError(z, "airyai(z) not defined at infinity")) + end x, y = real(z), imag(z) zabs = abs(z) @@ -70,17 +81,31 @@ function airyaiprime2(z::ComplexOrReal{T}) where T # use computation for real numbers then convert to input type for stability # for imag(z) != 0 use rotation identity (http://dlmf.nist.gov/9.2.E14) if iszero(y) - xx = 2 * x^(T(3)/2) / 3 + xabs = abs(x) + xx = 2 * xabs^(T(3)/2) / 3 Jv, Yv = besseljy_positive_args(T(2)/3, xx) Jmv = -(Jv + sqrt(T(3))*Yv) / 2 - return convert(eltype(z), abs(x) * (Jv - Jmv) / 3 + return convert(eltype(z), xabs * (Jv - Jmv) / 3) else - return -(exp(2pi*im/3)*airyaiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*airyaiprime(-z*exp(-pi*im/3))) + return -(exp(2pi*im/3)*_airyaiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*_airyaiprime(-z*exp(-pi*im/3))) end end end -function airybi2(z::ComplexOrReal{T}) where T +""" + airybi(x) +Airy function of the second kind ``\\operatorname{Bi}(x)``. +""" +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} + if ~isfinite(z) + isnan(z) && return z + isinf(z) && return throw(DomainError(z, "airyai(z) not defined at infinity")) + end x, y = real(z), imag(z) zabs = abs(z) @@ -100,16 +125,28 @@ function airybi2(z::ComplexOrReal{T}) where T return sqrt(z/3) * (inu + inu2) else if iszero(y) - xx = 2 * x^(T(3)/2) / 3 + xabs = abs(x) + xx = 2 * xabs^(T(3)/2) / 3 Jv, Yv = besseljy_positive_args(one(T)/3, xx) Jmv = (Jv - sqrt(T(3)) * Yv) / 2 - return convert(eltype(z), sqrt(abs(x)/3) * (Jmv - Jv)) + return convert(eltype(z), sqrt(xabs/3) * (Jmv - Jv)) else - return exp(pi*im/3) * airybi2(-z*exp(pi*im/3)) + exp(-pi*im/3) * airybi2(-z*exp(-pi*im/3)) + return exp(pi*im/3) * _airybi(-z*exp(pi*im/3)) + exp(-pi*im/3) * _airybi(-z*exp(-pi*im/3)) end + end end -function airybiprime2(z::ComplexOrReal{T}) where T +""" + airybiprime(x) +Derivative of the Airy function of the second kind ``\\operatorname{Bi}'(x)``. +""" + +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} x, y = real(z), imag(z) zabs = abs(z) @@ -129,13 +166,15 @@ function airybiprime2(z::ComplexOrReal{T}) where T return z / sqrt(3) * (inu + inu2) else if iszero(y) - xx = 2 * x^(T(3)/2) / 3 + xabs = abs(x) + xx = 2 * xabs^(T(3)/2) / 3 Jv, Yv = besseljy_positive_args(T(2)/3, xx) Jmv = -(Jv + sqrt(T(3))*Yv) / 2 - return convert(eltype(z), abs(x) * (Jv + Jmv) / sqrt(T(3))) + return convert(eltype(z), xabs * (Jv + Jmv) / sqrt(T(3))) else - return -(exp(2pi*im/3)*airybiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*airybiprime(-z*exp(-pi*im/3))) + return -(exp(2pi*im/3)*_airybiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*_airybiprime(-z*exp(-pi*im/3))) end + end end ##### @@ -148,6 +187,7 @@ airyai_power_series_cutoff(x::T, y::T) where T <: Float32 = x < 4.5f0 && abs(y) function airyai_power_series(x::ComplexOrReal{T}) where T S = eltype(x) + iszero(x) && return S(0.3550280538878172) MaxIter = 3000 ai1 = zero(S) ai2 = zero(S) @@ -172,6 +212,7 @@ end function airyaiprime_power_series(x::ComplexOrReal{T}) where T S = eltype(x) + iszero(x) && return S(-0.2588194037928068) MaxIter = 3000 ai1 = zero(S) ai2 = zero(S) @@ -203,6 +244,7 @@ airybi_power_series_cutoff(x::T, y::T) where T <: Float32 = abs(complex(x, y)) < function airybi_power_series(x::ComplexOrReal{T}) where T S = eltype(x) + iszero(x) && return S(0.6149266274460007) MaxIter = 3000 ai1 = zero(S) ai2 = zero(S) @@ -222,6 +264,7 @@ function airybi_power_series(x::ComplexOrReal{T}) where T end function airybiprime_power_series(x::ComplexOrReal{T}) where T S = eltype(x) + iszero(x) && return S(0.4482883573538264) MaxIter = 3000 ai1 = zero(S) ai2 = zero(S) @@ -290,6 +333,7 @@ function airybi_large_argument(z::Complex{T}) where T check_conj = false if y < zero(T) z = conj(z) + b = conj(b) y = abs(y) check_conj = true end @@ -320,6 +364,7 @@ function airybiprime_large_argument(z::Complex{T}) where T check_conj = false if y < zero(T) z = conj(z) + d = conj(d) y = abs(y) check_conj = true end diff --git a/src/math_constants.jl b/src/math_constants.jl index f98e221..0e30f0c 100644 --- a/src/math_constants.jl +++ b/src/math_constants.jl @@ -22,3 +22,8 @@ const TWOOPI(::Type{Float32}) = 0.636619772367581343075535f0 const THPIO4(::Type{Float32}) = 2.35619449019234492885f0 const SQ2O2(::Type{Float32}) = 0.7071067811865476f0 const SQPIO2(::Type{Float32}) = 1.25331413731550025f0 + +const GAMMA_TWO_THIRDS(::Type{Float64}) = 1.3541179394264005 +const GAMMA_ONE_THIRD(::Type{Float64}) = 2.6789385347077475 +const GAMMA_TWO_THIRDS(::Type{Float32}) = 1.3541179394264005f0 +const GAMMA_ONE_THIRD(::Type{Float32}) = 2.6789385347077475f0 diff --git a/test/airy_test.jl b/test/airy_test.jl index aa09707..e78d751 100644 --- a/test/airy_test.jl +++ b/test/airy_test.jl @@ -3,24 +3,35 @@ # 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=1e-12) - @test isapprox(airyai(-x), SpecialFunctions.airyai(-x), rtol=1e-9) + @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=1e-12) - @test isapprox(airyaiprime(-x), SpecialFunctions.airyaiprime(-x), rtol=1e-9) + @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=1e-12) - @test isapprox(airybi(-x), SpecialFunctions.airybi(-x), rtol=5e-8) + @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=1e-12) - @test isapprox(airybiprime(-x), SpecialFunctions.airybiprime(-x), rtol=1e-9) + @test isapprox(airybiprime(x), SpecialFunctions.airybiprime(x), rtol=2e-13) + @test isapprox(airybiprime(-x), SpecialFunctions.airybiprime(-x), rtol=5e-12) +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 + z = x*exp(im*a) + @test isapprox(airyai(z), SpecialFunctions.airyai(z), rtol=5e-10) + @test isapprox(airyaiprime(z), SpecialFunctions.airyaiprime(z), rtol=5e-10) + @test isapprox(airybi(z), SpecialFunctions.airybi(z), rtol=1e-11) + @test isapprox(airybiprime(z), SpecialFunctions.airybiprime(z), rtol=1e-11) end # test Inf +# no longer set up +#= @test iszero(airyai(Inf)) @test iszero(airyaiprime(Inf)) @test isinf(airybi(Inf)) @test isinf(airybiprime(Inf)) +=# # test Float16 types @test airyai(Float16(1.2)) isa Float16 From 9d7227e8b2fd9d329c438e95680907f461653cbb Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Sat, 24 Sep 2022 18:23:54 -0400 Subject: [PATCH 05/14] fix type promotions --- src/besseli.jl | 4 ++-- src/besselk.jl | 6 +++++- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index 272c1c5..b4cae34 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -306,9 +306,9 @@ end besseli_large_argument_scaled(v, x::T) where T = T(_besseli_large_argument(v, x) / sqrt(2 * (π * x))) -function _besseli_large_argument(v, x::ComplexOrReal{T}) where T +function _besseli_large_argument(v, x::T) where T MaxIter = 5000 - S = eltype(x) + S = promote_type(T, Float64) v, x = S(v), S(x) fv2 = 4 * v^2 diff --git a/src/besselk.jl b/src/besselk.jl index 8e87baa..eca9deb 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -362,10 +362,12 @@ end # computes K_{nu+1}/K_{nu} using continued fractions and the modified Lentz method # generally slow to converge for small x +besselk_ratio_knu_knup1(v, x::Float32) = Float32(besselk_ratio_knu_knup1(v, Float64(x))) +besselk_ratio_knu_knup1(v, x::ComplexF32) = ComplexF32(besselk_ratio_knu_knup1(v, ComplexF64(x))) function besselk_ratio_knu_knup1(v, x::ComplexOrReal{T}) where T MaxIter = 1000 S = eltype(x) - (hn, Dn, Cn) = (S(1e-150), zero(S), S(1e-150)) + (hn, Dn, Cn) = (S(1e-50), zero(S), S(1e-50)) jf = one(S) vv = v * v @@ -457,6 +459,8 @@ end besselk_large_argument_scaled(v, x::T) where T = T(_besselk_large_argument(v, x) * sqrt(pi / 2x)) +_besselk_large_argument(v, x::Float32) = Float32(_besselk_large_argument(v, Float64(x))) +_besselk_large_argument(v, x::ComplexF32) = ComplexF32(_besselk_large_argument(v, ComplexF64(x))) function _besselk_large_argument(v, x::ComplexOrReal{T}) where T MaxIter = 5000 S = eltype(x) From a2ff7b857113056e6a6e56d4df0d86d2801e7927 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Sat, 24 Sep 2022 21:40:24 -0400 Subject: [PATCH 06/14] remove old functions --- src/besseli.jl | 39 --------------------- src/gamma.jl | 93 -------------------------------------------------- 2 files changed, 132 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index b4cae34..3ac648f 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -395,42 +395,3 @@ function steed(n, x::T) where T return h end =# -function _besseli(nu, x::Complex{T}) where T <: Union{Float32, Float64} - isinteger(nu) && return _besseli(Int(nu), x) - ~isfinite(x) && return x - abs_nu = abs(nu) - abs_x = abs(x) - - if nu >= 0 - if x >= 0 - return besseli_positive_args(abs_nu, abs_x) - else - return cispi(abs_nu) * besseli_positive_args(abs_nu, abs_x) - end - else - return throw(DomainError(nu, "nu must be positive if x is complex")) - #= - if x >= 0 - return besseli_positive_args(abs_nu, abs_x) + 2 / π * sinpi(abs_nu) * besselk_positive_args(abs_nu, abs_x) - else - Iv = besseli_positive_args(abs_nu, abs_x) - Kv = besselk_positive_args(abs_nu, abs_x) - return cispi(abs_nu) * Iv + 2 / π * sinpi(abs_nu) * (cispi(-abs_nu) * Kv - im * π * Iv) - end - =# - end -end -function _besseli(nu::Integer, x::Complex{T}) where T <: Union{Float32, Float64} - ~isfinite(x) && return x - abs_nu = abs(nu) - abs_x = abs(x) - sg = iseven(abs_nu) ? 1 : -1 - - if x >= 0 - return besseli_positive_args(abs_nu, abs_x) - else - return sg * besseli_positive_args(abs_nu, abs_x) - end -end - -Base.eps(::Type{Complex{T}}) where {T<:AbstractFloat} = eps(T) \ No newline at end of file diff --git a/src/gamma.jl b/src/gamma.jl index a0668bc..28ed201 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -67,96 +67,3 @@ function small_gamma(x) return z * p / q end -_gamma(z::Complex) = exp(loggamma(z)) - -loggamma(z::Complex) = _loggamma(float(z)) - -_loggamma(x::Complex{Float32}) = ComplexF32(gamma(ComplexF64(x))) - -# copied from https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/src/gamma.jl - -# Compute the logΓ(z) function using a combination of the asymptotic series, -# the Taylor series around z=1 and z=2, the reflection formula, and the shift formula. -# Many details of these techniques are discussed in D. E. G. Hare, -# "Computing the principal branch of log-Gamma," J. Algorithms 25, pp. 221-236 (1997), -# and similar techniques are used (in a somewhat different way) by the -# SciPy loggamma function. The key identities are also described -# at http://functions.wolfram.com/GammaBetaErf/LogGamma/ -function _loggamma(z::Complex{Float64}) - x, y = reim(z) - yabs = abs(y) - if !isfinite(x) || !isfinite(y) # Inf or NaN - if isinf(x) && isfinite(y) - return Complex(x, x > 0 ? (y == 0 ? y : copysign(Inf, y)) : copysign(Inf, -y)) - elseif isfinite(x) && isinf(y) - return Complex(-Inf, y) - else - return Complex(NaN, NaN) - end - elseif x > 7 || yabs > 7 # use the Stirling asymptotic series for sufficiently large x or |y| - return loggamma_asymptotic(z) - elseif x < 0.1 # use reflection formula to transform to x > 0 - if x == 0 && y == 0 # return Inf with the correct imaginary part for z == 0 - return Complex(Inf, signbit(x) ? copysign(Float64(π), -y) : -y) - end - # the 2pi * floor(...) stuff is to choose the correct branch cut for log(sinpi(z)) - return Complex(Float64(logπ), copysign(Float64(twoπ), y) * floor(0.5*x+0.25)) - - log(sinpi(z)) - loggamma(1-z) - elseif abs(x - 1) + yabs < 0.1 - # taylor series at z=1 - # ... coefficients are [-eulergamma; [(-1)^k * zeta(k)/k for k in 2:15]] - w = Complex(x - 1, y) - return w * @evalpoly(w, -5.7721566490153286060651188e-01,8.2246703342411321823620794e-01, - -4.0068563438653142846657956e-01,2.705808084277845478790009e-01, - -2.0738555102867398526627303e-01,1.6955717699740818995241986e-01, - -1.4404989676884611811997107e-01,1.2550966952474304242233559e-01, - -1.1133426586956469049087244e-01,1.000994575127818085337147e-01, - -9.0954017145829042232609344e-02,8.3353840546109004024886499e-02, - -7.6932516411352191472827157e-02,7.1432946295361336059232779e-02, - -6.6668705882420468032903454e-02) - elseif abs(x - 2) + yabs < 0.1 - # taylor series at z=2 - # ... coefficients are [1-eulergamma; [(-1)^k * (zeta(k)-1)/k for k in 2:12]] - w = Complex(x - 2, y) - return w * @evalpoly(w, 4.2278433509846713939348812e-01,3.2246703342411321823620794e-01, - -6.7352301053198095133246196e-02,2.0580808427784547879000897e-02, - -7.3855510286739852662729527e-03,2.8905103307415232857531201e-03, - -1.1927539117032609771139825e-03,5.0966952474304242233558822e-04, - -2.2315475845357937976132853e-04,9.9457512781808533714662972e-05, - -4.4926236738133141700224489e-05,2.0507212775670691553131246e-05) - end - # use recurrence relation loggamma(z) = loggamma(z+1) - log(z) to shift to x > 7 for asymptotic series - shiftprod = Complex(x,yabs) - x += 1 - sb = false # == signbit(imag(shiftprod)) == signbit(yabs) - # To use log(product of shifts) rather than sum(logs of shifts), - # we need to keep track of the number of + to - sign flips in - # imag(shiftprod), as described in Hare (1997), proposition 2.2. - signflips = 0 - while x <= 7 - shiftprod *= Complex(x,yabs) - sb′ = signbit(imag(shiftprod)) - signflips += sb′ & (sb′ != sb) - sb = sb′ - x += 1 - end - shift = log(shiftprod) - if signbit(y) # if y is negative, conjugate the shift - shift = Complex(real(shift), signflips*-6.2831853071795864769252842 - imag(shift)) - else - shift = Complex(real(shift), imag(shift) + signflips*6.2831853071795864769252842) - end - return loggamma_asymptotic(Complex(x,y)) - shift -end - -# asymptotic series for log(gamma(z)), valid for sufficiently large real(z) or |imag(z)| -function loggamma_asymptotic(z::Complex{Float64}) - zinv = inv(z) - t = zinv*zinv - # coefficients are bernoulli[2:n+1] .// (2*(1:n).*(2*(1:n) - 1)) - return (z-0.5)*log(z) - z + 9.1893853320467274178032927e-01 + # <-- log(2pi)/2 - zinv*@evalpoly(t, 8.3333333333333333333333368e-02,-2.7777777777777777777777776e-03, - 7.9365079365079365079365075e-04,-5.9523809523809523809523806e-04, - 8.4175084175084175084175104e-04,-1.9175269175269175269175262e-03, - 6.4102564102564102564102561e-03,-2.9550653594771241830065352e-02) -end \ No newline at end of file From 6e72ef621d2faaba8fa1e4e831dcbbdaa8118427 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Sat, 24 Sep 2022 21:59:30 -0400 Subject: [PATCH 07/14] remove gamma calls --- src/airy.jl | 12 ++++++++---- src/gamma.jl | 1 - src/math_constants.jl | 5 +++++ 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/airy.jl b/src/airy.jl index aa46e9c..8ac2d4c 100644 --- a/src/airy.jl +++ b/src/airy.jl @@ -387,7 +387,7 @@ function airy_large_arg_a(x::ComplexOrReal{T}) where T xsqr = sqrt(x) out = zero(S) - t = gamma(one(T) / 6) * gamma(T(5) / 6) / 4 + t = GAMMA_ONE_SIXTH(T) * GAMMA_FIVE_SIXTHS(T) / 4 a = 4*xsqr*x for i in 0:MaxIter out += t @@ -403,7 +403,7 @@ function airy_large_arg_b(x::ComplexOrReal{T}) where T xsqr = sqrt(x) out = zero(S) - t = gamma(one(T) / 6) * gamma(T(5) / 6) / 4 + t = GAMMA_ONE_SIXTH(T) * GAMMA_FIVE_SIXTHS(T) / 4 a = 4*xsqr*x for i in 0:MaxIter out += t @@ -419,7 +419,9 @@ function airy_large_arg_c(x::ComplexOrReal{T}) where T xsqr = sqrt(x) out = zero(S) - t = gamma(-one(T) / 6) * gamma(T(7) / 6) / 4 + # use identities of gamma to relate to defined constants + # t = gamma(-one(T) / 6) * gamma(T(7) / 6) / 4 + t = -GAMMA_FIVE_SIXTHS(T) * GAMMA_ONE_SIXTH(T) / 4 a = 4*xsqr*x for i in 0:MaxIter out += t @@ -435,7 +437,9 @@ function airy_large_arg_d(x::ComplexOrReal{T}) where T xsqr = sqrt(x) out = zero(S) - t = gamma(-one(T) / 6) * gamma(T(7) / 6) / 4 + # use identities of gamma to relate to defined constants + # t = gamma(-one(T) / 6) * gamma(T(7) / 6) / 4 + t = -GAMMA_FIVE_SIXTHS(T) * GAMMA_ONE_SIXTH(T) / 4 a = 4*xsqr*x for i in 0:MaxIter out += t diff --git a/src/gamma.jl b/src/gamma.jl index 28ed201..25aa086 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -66,4 +66,3 @@ function small_gamma(x) q = evalpoly(x, Q) return z * p / q end - diff --git a/src/math_constants.jl b/src/math_constants.jl index 0e30f0c..30d06ea 100644 --- a/src/math_constants.jl +++ b/src/math_constants.jl @@ -25,5 +25,10 @@ const SQPIO2(::Type{Float32}) = 1.25331413731550025f0 const GAMMA_TWO_THIRDS(::Type{Float64}) = 1.3541179394264005 const GAMMA_ONE_THIRD(::Type{Float64}) = 2.6789385347077475 +const GAMMA_ONE_SIXTH(::Type{Float64}) = 5.566316001780235 +const GAMMA_FIVE_SIXTHS(::Type{Float64}) = 1.128787029908126 + const GAMMA_TWO_THIRDS(::Type{Float32}) = 1.3541179394264005f0 const GAMMA_ONE_THIRD(::Type{Float32}) = 2.6789385347077475f0 +const GAMMA_ONE_SIXTH(::Type{Float32}) = 5.566316001780235f0 +const GAMMA_FIVE_SIXTHS(::Type{Float32}) = 1.128787029908126f0 From b838a200adf7ca05fda367d43a50e392435c9df5 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 26 Sep 2022 16:39:47 -0400 Subject: [PATCH 08/14] add definitions at infinity --- src/airy.jl | 31 ++++++++++++++++++++++++++++--- test/airy_test.jl | 14 +++++++++++--- 2 files changed, 39 insertions(+), 6 deletions(-) diff --git a/src/airy.jl b/src/airy.jl index 8ac2d4c..c789516 100644 --- a/src/airy.jl +++ b/src/airy.jl @@ -22,7 +22,11 @@ _airyai(z::ComplexF16) = ComplexF16(_airyai(ComplexF32(z))) function _airyai(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} if ~isfinite(z) isnan(z) && return z - isinf(z) && return throw(DomainError(z, "airyai(z) not defined at infinity")) + if abs(angle(z)) < 2π/3 + return exp(-z) + else + return 1 / z + end end x, y = real(z), imag(z) zabs = abs(z) @@ -63,7 +67,11 @@ _airyaiprime(z::ComplexF16) = ComplexF16(_airyaiprime(ComplexF32(z))) function _airyaiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} if ~isfinite(z) isnan(z) && return z - isinf(z) && return throw(DomainError(z, "airyai(z) not defined at infinity")) + if abs(angle(z)) < 2π/3 + return -exp(-z) + else + return 1 / z + end end x, y = real(z), imag(z) zabs = abs(z) @@ -104,7 +112,11 @@ _airybi(z::ComplexF16) = ComplexF16(_airybi(ComplexF32(z))) function _airybi(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} if ~isfinite(z) isnan(z) && return z - isinf(z) && return throw(DomainError(z, "airyai(z) not defined at infinity")) + if abs(angle(z)) < 2π/3 + return exp(z) + else + return 1 / z + end end x, y = real(z), imag(z) zabs = abs(z) @@ -147,6 +159,14 @@ _airybiprime(x::Float16) = Float16(_airybiprime(Float32(x))) _airybiprime(z::ComplexF16) = ComplexF16(_airybiprime(ComplexF32(z))) function _airybiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} + if ~isfinite(z) + isnan(z) && return z + if abs(angle(z)) < 2π/3 + return exp(z) + else + return -1 / z + end + end x, y = real(z), imag(z) zabs = abs(z) @@ -262,6 +282,11 @@ function airybi_power_series(x::ComplexOrReal{T}) where T end return (ai1*3^(-T(1)/6) + ai2*3^(-T(5)/6)) end + +##### +##### Power series for airybiprime(x) +##### + function airybiprime_power_series(x::ComplexOrReal{T}) where T S = eltype(x) iszero(x) && return S(0.4482883573538264) diff --git a/test/airy_test.jl b/test/airy_test.jl index e78d751..03ac49a 100644 --- a/test/airy_test.jl +++ b/test/airy_test.jl @@ -25,13 +25,21 @@ 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 end # test Inf -# no longer set up -#= + @test iszero(airyai(Inf)) @test iszero(airyaiprime(Inf)) @test isinf(airybi(Inf)) @test isinf(airybiprime(Inf)) -=# + +@test airyai(Inf + 0.0im) === exp(-(Inf + 0.0im)) +@test airyaiprime(Inf + 0.0im) === -exp(-(Inf + 0.0im)) +@test airyai(-Inf + 0.0im) === 1 / (-Inf + 0.0im) +@test airyaiprime(-Inf + 0.0im) === 1 / (-Inf + 0.0im) +@test airybi(Inf + Inf*im) === exp((Inf + Inf*im)) +@test airybi(-Inf + 10.0*im) === 1 / (-Inf + 10.0*im) +@test airybiprime(Inf + 0.0*im) === exp((Inf + 0.0*im)) +@test airybiprime(-Inf + 0.0*im) === -1 / (-Inf + 0.0*im) + # test Float16 types @test airyai(Float16(1.2)) isa Float16 From b0390d599c1572738d2ac0786c2ac82af94f190e Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 26 Sep 2022 16:56:16 -0400 Subject: [PATCH 09/14] update docs and add reference --- src/airy.jl | 30 +++++++++++++++++++++++------- src/besseli.jl | 2 +- 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/src/airy.jl b/src/airy.jl index c789516..df49d85 100644 --- a/src/airy.jl +++ b/src/airy.jl @@ -1,14 +1,25 @@ # Airy functions # -# airyai(x), airybi(nu, x) -# airyaiprime(x), airybiprime(nu, x) +# airyai(z), airybi(nu, z) +# airyaiprime(z), airybiprime(nu, z) # -# A numerical routine to compute the airy functions and their derivatives. -# These routines use their relations to other special functions using https://dlmf.nist.gov/9.6. -# Specifically see (NIST 9.6.E1 - 9.6.E9) for computation from the defined bessel functions. -# For negative arguments these definitions are prone to some cancellation leading to higher errors. -# In the future, these could be replaced with more custom routines as they depend on a single variable. +# 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 +# for small arguments and a large argument expansion for (x > ~10). The primary difference between [1] +# and what is used here is that the regions where the power series and large argument expansions +# do not provide good results they are filled by relation to other special functions (besselk and besseli) +# using https://dlmf.nist.gov/9.6 (NIST 9.6.E1 - 9.6.E9). In this case the power series of besseli is used and then besselk +# is calculated using the continued fraction approach. This method is described in more detail in src/besselk.jl. +# However, care must be taken when computing besseli because when the imaginary component is much larger than the real part +# cancellation will occur. This can be overcome by shifting the order of besseli to be much larger and then using the power series +# and downward recurrence to get besseli(1/3, x). Another difficult region is when -10 Date: Tue, 27 Sep 2022 11:38:34 -0400 Subject: [PATCH 10/14] add Float32 tests --- src/airy.jl | 63 +++++++++++++++++++++++++++++++---------------- test/airy_test.jl | 19 ++++++++++++++ 2 files changed, 61 insertions(+), 21 deletions(-) diff --git a/src/airy.jl b/src/airy.jl index df49d85..c7149ab 100644 --- a/src/airy.jl +++ b/src/airy.jl @@ -216,7 +216,7 @@ end airyai_power_series_cutoff(x::T, y::T) where T <: Float64 = x < 2 && abs(y) > -1.4*(x + 5.5) airyai_power_series_cutoff(x::T, y::T) where T <: Float32 = x < 4.5f0 && abs(y) > -1.4f0*(x + 9.5f0) -function airyai_power_series(x::ComplexOrReal{T}) where T +function airyai_power_series(x::ComplexOrReal{T}; tol=eps(T)) where T S = eltype(x) iszero(x) && return S(0.3550280538878172) MaxIter = 3000 @@ -230,18 +230,20 @@ function airyai_power_series(x::ComplexOrReal{T}) where T for i in 0:MaxIter ai1 += t ai2 += t2 - abs(t) < eps(T) * abs(ai1) && break + abs(t) < tol * abs(ai1) && break t *= x3 * inv(9*(i + one(T))*(i + T(2)/3)) t2 *= x3 * inv(9*(i + one(T))*(i + T(4)/3)) end return (ai1*3^(-T(2)/3) - ai2*3^(-T(4)/3)) end +airyai_power_series(x::Float32) = Float32(airyai_power_series(Float64(x), tol=eps(Float32))) +airyai_power_series(x::ComplexF32) = ComplexF32(airyai_power_series(ComplexF64(x), tol=eps(Float32))) ##### ##### Power series for airyaiprime(x) ##### -function airyaiprime_power_series(x::ComplexOrReal{T}) where T +function airyaiprime_power_series(x::ComplexOrReal{T}; tol=eps(T)) where T S = eltype(x) iszero(x) && return S(-0.2588194037928068) MaxIter = 3000 @@ -255,12 +257,14 @@ function airyaiprime_power_series(x::ComplexOrReal{T}) where T for i in 0:MaxIter ai1 += t ai2 += t2 - abs(t) < eps(T) * abs(ai1) && break + abs(t) < tol * abs(ai1) && break t *= x3 * inv(9*(i + one(T))*(i + T(1)/3)) t2 *= x3 * inv(9*(i + one(T))*(i + T(5)/3)) end return -ai1*3^(-T(1)/3) + ai2*3^(-T(5)/3) end +airyaiprime_power_series(x::Float32) = Float32(airyaiprime_power_series(Float64(x), tol=eps(Float32))) +airyaiprime_power_series(x::ComplexF32) = ComplexF32(airyaiprime_power_series(ComplexF64(x), tol=eps(Float32))) ##### ##### Power series for airybi(x) @@ -268,12 +272,12 @@ end # cutoffs for power series valid for both airybi and airybiprime # has a more complicated validity as it works well close to positive real line and for small negative arguments also works for angle(z) ~ 2pi/3 -# the statements are somewhat complicated but we absolutely want to hit this branch when we can as the other algorithms are 10x slower +# the statements are somewhat complicated but we want to hit this branch when we can as the other algorithms are 10x slower # the Float32 cutoff can be simplified because it overlaps with the large argument expansion so there isn't a need for more complicated statements airybi_power_series_cutoff(x::T, y::T) where T <: Float64 = (abs(y) > -1.4*(x + 5.5) && abs(y) < -2.2*(x - 4)) || (x > zero(T) && abs(y) < 3) airybi_power_series_cutoff(x::T, y::T) where T <: Float32 = abs(complex(x, y)) < 9 -function airybi_power_series(x::ComplexOrReal{T}) where T +function airybi_power_series(x::ComplexOrReal{T}; tol=eps(T)) where T S = eltype(x) iszero(x) && return S(0.6149266274460007) MaxIter = 3000 @@ -287,18 +291,20 @@ function airybi_power_series(x::ComplexOrReal{T}) where T for i in 0:MaxIter ai1 += t ai2 += t2 - abs(t) < eps(T) * abs(ai1) && break + abs(t) < tol * abs(ai1) && break t *= x3 * inv(9*(i + one(T))*(i + T(2)/3)) t2 *= x3 * inv(9*(i + one(T))*(i + T(4)/3)) end return (ai1*3^(-T(1)/6) + ai2*3^(-T(5)/6)) end +airybi_power_series(x::Float32) = Float32(airybi_power_series(Float64(x), tol=eps(Float32))) +airybi_power_series(x::ComplexF32) = ComplexF32(airybi_power_series(ComplexF64(x), tol=eps(Float32))) ##### ##### Power series for airybiprime(x) ##### -function airybiprime_power_series(x::ComplexOrReal{T}) where T +function airybiprime_power_series(x::ComplexOrReal{T}; tol=eps(T)) where T S = eltype(x) iszero(x) && return S(0.4482883573538264) MaxIter = 3000 @@ -312,12 +318,14 @@ function airybiprime_power_series(x::ComplexOrReal{T}) where T for i in 0:MaxIter ai1 += t ai2 += t2 - abs(t) < eps(T) * abs(ai1) && break + abs(t) < tol * abs(ai1) && break t *= x3 * inv(9*(i + one(T))*(i + T(1)/3)) t2 *= x3 * inv(9*(i + one(T))*(i + T(5)/3)) end return (ai1*3^(T(1)/6) + ai2*3^(-T(7)/6)) end +airybiprime_power_series(x::Float32) = Float32(airybiprime_power_series(Float64(x), tol=eps(Float32))) +airybiprime_power_series(x::ComplexF32) = ComplexF32(airybiprime_power_series(ComplexF64(x), tol=eps(Float32))) ##### ##### Large argument expansion for airy functions @@ -422,7 +430,7 @@ function airybiprime_large_argument(z::Complex{T}) where T end # see equations 24 and relations using eq 25 and 26 in [1] -function airy_large_arg_a(x::ComplexOrReal{T}) where T +function airy_large_arg_a(x::ComplexOrReal{T}; tol=eps(T)*40) where T S = eltype(x) MaxIter = 3000 xsqr = sqrt(x) @@ -432,13 +440,13 @@ function airy_large_arg_a(x::ComplexOrReal{T}) where T a = 4*xsqr*x for i in 0:MaxIter out += t - abs(t) < eps(T)*50 * abs(out) && break + abs(t) < tol*abs(out) && break t *= -3*(i + one(T)/6) * (i + T(5)/6) / (a*(i + one(T))) end - return out * exp(-a / 6) / (pi^(3/2) * sqrt(xsqr)) + return out * exp(-a / 6) / (π^(T(3)/2) * sqrt(xsqr)) end -function airy_large_arg_b(x::ComplexOrReal{T}) where T +function airy_large_arg_b(x::ComplexOrReal{T}; tol=eps(T)*40) where T S = eltype(x) MaxIter = 3000 xsqr = sqrt(x) @@ -448,13 +456,13 @@ function airy_large_arg_b(x::ComplexOrReal{T}) where T a = 4*xsqr*x for i in 0:MaxIter out += t - abs(t) < eps(T)*50 * abs(out) && break + abs(t) < tol*abs(out) && break t *= 3*(i + one(T)/6) * (i + T(5)/6) / (a*(i + one(T))) end - return out * exp(a / 6) / (pi^(3/2) * sqrt(xsqr)) + return out * exp(a / 6) / (π^(T(3)/2) * sqrt(xsqr)) end -function airy_large_arg_c(x::ComplexOrReal{T}) where T +function airy_large_arg_c(x::ComplexOrReal{T}; tol=eps(T)*40) where T S = eltype(x) MaxIter = 3000 xsqr = sqrt(x) @@ -466,13 +474,13 @@ function airy_large_arg_c(x::ComplexOrReal{T}) where T a = 4*xsqr*x for i in 0:MaxIter out += t - abs(t) < eps(T)*50* abs(out) && break + abs(t) < tol*abs(out) && break t *= -3*(i - one(T)/6) * (i + T(7)/6) / (a*(i + one(T))) end - return out * exp(-a / 6) * sqrt(xsqr) / pi^(3/2) + return out * exp(-a / 6) * sqrt(xsqr) / π^(T(3)/2) end -function airy_large_arg_d(x::ComplexOrReal{T}) where T +function airy_large_arg_d(x::ComplexOrReal{T}; tol=eps(T)*40) where T S = eltype(x) MaxIter = 3000 xsqr = sqrt(x) @@ -484,12 +492,25 @@ function airy_large_arg_d(x::ComplexOrReal{T}) where T a = 4*xsqr*x for i in 0:MaxIter out += t - abs(t) < eps(T)*50 * abs(out) && break + abs(t) < tol*abs(out) && break t *= 3*(i - one(T)/6) * (i + T(7)/6) / (a*(i + one(T))) end - return -out * exp(a / 6) * sqrt(xsqr) / pi^(3/2) + return -out * exp(a / 6) * sqrt(xsqr) / π^(T(3)/2) end +# negative arguments of airy functions oscillate around zero so as x -> -Inf it is more prone to cancellation +# to give best accuracy it is best to promote to Float64 numbers until the Float32 tolerance +airy_large_arg_a(x::Float32) = (airy_large_arg_a(Float64(x), tol=eps(Float32))) +airy_large_arg_a(x::ComplexF32) = (airy_large_arg_a(ComplexF64(x), tol=eps(Float32))) + +airy_large_arg_b(x::Float32) = Float32(airy_large_arg_b(Float64(x), tol=eps(Float32))) +airy_large_arg_b(x::ComplexF32) = ComplexF32(airy_large_arg_b(ComplexF64(x), tol=eps(Float32))) + +airy_large_arg_c(x::Float32) = Float32(airy_large_arg_c(Float64(x), tol=eps(Float32))) +airy_large_arg_c(x::ComplexF32) = ComplexF32(airy_large_arg_c(ComplexF64(x), tol=eps(Float32))) + +airy_large_arg_d(x::Float32) = Float32(airy_large_arg_d(Float64(x), tol=eps(Float32))) +airy_large_arg_d(x::ComplexF32) = ComplexF32(airy_large_arg_d(ComplexF64(x), tol=eps(Float32))) # calculates besselk from the power series of besseli using the continued fraction and wronskian # this shift the order higher first to avoid cancellation in the power series of besseli along the imaginary axis # for real arguments this is not needed because besseli can be computed stably over the entire real axis diff --git a/test/airy_test.jl b/test/airy_test.jl index 03ac49a..e63f0b4 100644 --- a/test/airy_test.jl +++ b/test/airy_test.jl @@ -16,6 +16,21 @@ for x in [0.0; 1e-17:0.1:100.0] @test isapprox(airybiprime(-x), SpecialFunctions.airybiprime(-x), rtol=5e-12) 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) +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 z = x*exp(im*a) @test isapprox(airyai(z), SpecialFunctions.airyai(z), rtol=5e-10) @@ -43,6 +58,10 @@ 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 From c195ec6479224ecd0679db4d4c8c6483f42c6ac0 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Tue, 27 Sep 2022 11:57:37 -0400 Subject: [PATCH 11/14] add news and vers. bump --- NEWS.md | 7 ++++++- Project.toml | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index bd7f17f..575cc8d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,7 +8,12 @@ For pre-1.0.0 releases we will increment the minor version when we export any ne For bug fixes, performance enhancements, or fixes to unexported functions we will increment the patch version. **Note**: The exported API can be considered very stable and likely will not change without serious consideration. -# Unreleased +# Unrealeased + +# Version 0.2.2 + +### Added + - Support for all airy functions and derivatives in entire complex plane ([PR #51](https://github.com/JuliaMath/Bessels.jl/pull/51)). These are specialized routines for airy functions instead of relying on connection to besselk. These are a couple digits more accurate and faster than previous versions. # Version 0.2.1 diff --git a/Project.toml b/Project.toml index bf23bbf..264beb9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Bessels" uuid = "0e736298-9ec6-45e8-9647-e4fc86a2fe38" authors = ["Michael Helton and contributors"] -version = "0.2.1" +version = "0.2.2" [compat] julia = "1.6" From 3e2e295ae1409a33f4a4a0f69a10c0310a831a45 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Tue, 27 Sep 2022 12:12:39 -0400 Subject: [PATCH 12/14] fix typo --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 575cc8d..ee07c90 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,7 +8,7 @@ For pre-1.0.0 releases we will increment the minor version when we export any ne For bug fixes, performance enhancements, or fixes to unexported functions we will increment the patch version. **Note**: The exported API can be considered very stable and likely will not change without serious consideration. -# Unrealeased +# Unreleased # Version 0.2.2 From 4ad55a59f8835315f999238864f8b5511765b0a2 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 10 Oct 2022 13:24:52 -0400 Subject: [PATCH 13/14] add cispi and fix 3/2 --- src/airy.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/airy.jl b/src/airy.jl index c7149ab..39ee672 100644 --- a/src/airy.jl +++ b/src/airy.jl @@ -47,7 +47,7 @@ function _airyai(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} if x > zero(T) # use relation to besselk (http://dlmf.nist.gov/9.6.E1) - zz = 2 * z^(T(3)/2) / 3 + zz = 2 * z * sqrt(z) / 3 return sqrt(z / 3) * besselk_continued_fraction_shift(one(T)/3, zz) / π else # z is close to the negative real axis @@ -56,12 +56,12 @@ function _airyai(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} # for imag(z) != 0 use rotation identity (http://dlmf.nist.gov/9.2.E14) if iszero(y) xabs = abs(x) - xx = 2 * xabs^(T(3)/2) / 3 + xx = 2 * xabs * sqrt(xabs) / 3 Jv, Yv = besseljy_positive_args(one(T)/3, xx) Jmv = (Jv - sqrt(T(3)) * Yv) / 2 return convert(eltype(z), sqrt(xabs) * (Jmv + Jv) / 3) else - return exp(pi*im/3) * _airyai(-z*exp(pi*im/3)) + exp(-pi*im/3) * _airyai(-z*exp(-pi*im/3)) + return cispi(one(T)/3) * _airyai(-z*cispi(one(T)/3)) + cispi(-one(T)/3) * _airyai(-z*cispi(-one(T)/3)) end end end @@ -92,7 +92,7 @@ function _airyaiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} if x > zero(T) # use relation to besselk (http://dlmf.nist.gov/9.6.E2) - zz = 2 * z^(T(3)/2) / 3 + zz = 2 * z * sqrt(z) / 3 return -z * besselk_continued_fraction_shift(T(2)/3, zz) / (π * sqrt(T(3))) else # z is close to the negative real axis @@ -101,12 +101,12 @@ function _airyaiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} # for imag(z) != 0 use rotation identity (http://dlmf.nist.gov/9.2.E14) if iszero(y) xabs = abs(x) - xx = 2 * xabs^(T(3)/2) / 3 + xx = 2 * xabs * sqrt(xabs) / 3 Jv, Yv = besseljy_positive_args(T(2)/3, xx) Jmv = -(Jv + sqrt(T(3))*Yv) / 2 return convert(eltype(z), xabs * (Jv - Jmv) / 3) else - return -(exp(2pi*im/3)*_airyaiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*_airyaiprime(-z*exp(-pi*im/3))) + return -(cispi(T(2)/3) * _airyaiprime(-z * cispi(one(T)/3)) + cispi(-T(2)/3) * _airyaiprime(-z * cispi(-one(T)/3))) end end end @@ -137,7 +137,7 @@ function _airybi(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} airybi_power_series_cutoff(x, y) && return airybi_power_series(z) if x > zero(T) - zz = 2 * z^(T(3)/2) / 3 + zz = 2 * z * sqrt(z) / 3 shift = 20 order = one(T)/3 inu, inum1 = besseli_power_series_inu_inum1(order + shift, zz) @@ -149,12 +149,12 @@ function _airybi(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} else if iszero(y) xabs = abs(x) - xx = 2 * xabs^(T(3)/2) / 3 + xx = 2 * xabs * sqrt(xabs) / 3 Jv, Yv = besseljy_positive_args(one(T)/3, xx) Jmv = (Jv - sqrt(T(3)) * Yv) / 2 return convert(eltype(z), sqrt(xabs/3) * (Jmv - Jv)) else - return exp(pi*im/3) * _airybi(-z*exp(pi*im/3)) + exp(-pi*im/3) * _airybi(-z*exp(-pi*im/3)) + return cispi(one(T)/3) * _airybi(-z * cispi(one(T)/3)) + cispi(-one(T)/3) * _airybi(-z*cispi(-one(T)/3)) end end end @@ -186,7 +186,7 @@ function _airybiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} airybi_power_series_cutoff(x, y) && return airybiprime_power_series(z) if x > zero(T) - zz = 2 * z^(T(3)/2) / 3 + zz = 2 * z * sqrt(z) / 3 shift = 20 order = T(2)/3 inu, inum1 = besseli_power_series_inu_inum1(order + shift, zz) @@ -198,12 +198,12 @@ function _airybiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} else if iszero(y) xabs = abs(x) - xx = 2 * xabs^(T(3)/2) / 3 + xx = 2 * xabs * sqrt(xabs) / 3 Jv, Yv = besseljy_positive_args(T(2)/3, xx) Jmv = -(Jv + sqrt(T(3))*Yv) / 2 return convert(eltype(z), xabs * (Jv + Jmv) / sqrt(T(3))) else - return -(exp(2pi*im/3)*_airybiprime(-z*exp(pi*im/3)) + exp(-2pi*im/3)*_airybiprime(-z*exp(-pi*im/3))) + return -(cispi(T(2)/3) * _airybiprime(-z*cispi(one(T)/3)) + cispi(-T(2)/3) * _airybiprime(-z*cispi(-one(T)/3))) end end end @@ -443,7 +443,7 @@ function airy_large_arg_a(x::ComplexOrReal{T}; tol=eps(T)*40) where T abs(t) < tol*abs(out) && break t *= -3*(i + one(T)/6) * (i + T(5)/6) / (a*(i + one(T))) end - return out * exp(-a / 6) / (π^(T(3)/2) * sqrt(xsqr)) + return out * exp(-a / 6) / (sqrt(T(π)^3) * sqrt(xsqr)) end function airy_large_arg_b(x::ComplexOrReal{T}; tol=eps(T)*40) where T @@ -459,7 +459,7 @@ function airy_large_arg_b(x::ComplexOrReal{T}; tol=eps(T)*40) where T abs(t) < tol*abs(out) && break t *= 3*(i + one(T)/6) * (i + T(5)/6) / (a*(i + one(T))) end - return out * exp(a / 6) / (π^(T(3)/2) * sqrt(xsqr)) + return out * exp(a / 6) / (sqrt(T(π)^3) * sqrt(xsqr)) end function airy_large_arg_c(x::ComplexOrReal{T}; tol=eps(T)*40) where T @@ -477,7 +477,7 @@ function airy_large_arg_c(x::ComplexOrReal{T}; tol=eps(T)*40) where T abs(t) < tol*abs(out) && break t *= -3*(i - one(T)/6) * (i + T(7)/6) / (a*(i + one(T))) end - return out * exp(-a / 6) * sqrt(xsqr) / π^(T(3)/2) + return out * exp(-a / 6) * sqrt(xsqr) / sqrt(T(π)^3) end function airy_large_arg_d(x::ComplexOrReal{T}; tol=eps(T)*40) where T @@ -495,7 +495,7 @@ function airy_large_arg_d(x::ComplexOrReal{T}; tol=eps(T)*40) where T abs(t) < tol*abs(out) && break t *= 3*(i - one(T)/6) * (i + T(7)/6) / (a*(i + one(T))) end - return -out * exp(a / 6) * sqrt(xsqr) / π^(T(3)/2) + return -out * exp(a / 6) * sqrt(xsqr) / sqrt(T(π)^3) end # negative arguments of airy functions oscillate around zero so as x -> -Inf it is more prone to cancellation From 9cd17e1be598fc8a430acde8fd2b9fdf25b5f4ee Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 10 Oct 2022 23:54:10 -0400 Subject: [PATCH 14/14] fix formating --- src/airy.jl | 24 +++++------------------- 1 file changed, 5 insertions(+), 19 deletions(-) diff --git a/src/airy.jl b/src/airy.jl index 39ee672..1c8f920 100644 --- a/src/airy.jl +++ b/src/airy.jl @@ -32,23 +32,20 @@ _airyai(z::ComplexF16) = ComplexF16(_airyai(ComplexF32(z))) function _airyai(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} if ~isfinite(z) - isnan(z) && return z - if abs(angle(z)) < 2π/3 + if abs(angle(z)) < 2*T(π)/3 return exp(-z) else return 1 / z end end x, y = real(z), imag(z) - zabs = abs(z) - airy_large_argument_cutoff(z) && return airyai_large_argument(z) airyai_power_series_cutoff(x, y) && return airyai_power_series(z) if x > zero(T) # use relation to besselk (http://dlmf.nist.gov/9.6.E1) zz = 2 * z * sqrt(z) / 3 - return sqrt(z / 3) * besselk_continued_fraction_shift(one(T)/3, zz) / π + return sqrt(z / 3) * besselk_continued_fraction_shift(one(T)/3, zz) / T(π) else # z is close to the negative real axis # for imag(z) == 0 use reflection to compute in terms of bessel functions of first kind (http://dlmf.nist.gov/9.6.E5) @@ -77,23 +74,20 @@ _airyaiprime(z::ComplexF16) = ComplexF16(_airyaiprime(ComplexF32(z))) function _airyaiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} if ~isfinite(z) - isnan(z) && return z - if abs(angle(z)) < 2π/3 + if abs(angle(z)) < 2*T(π)/3 return -exp(-z) else return 1 / z end end x, y = real(z), imag(z) - zabs = abs(z) - airy_large_argument_cutoff(z) && return airyaiprime_large_argument(z) airyai_power_series_cutoff(x, y) && return airyaiprime_power_series(z) if x > zero(T) # use relation to besselk (http://dlmf.nist.gov/9.6.E2) zz = 2 * z * sqrt(z) / 3 - return -z * besselk_continued_fraction_shift(T(2)/3, zz) / (π * sqrt(T(3))) + return -z * besselk_continued_fraction_shift(T(2)/3, zz) / (T(π) * sqrt(T(3))) else # z is close to the negative real axis # for imag(z) == 0 use reflection to compute in terms of bessel functions of first kind (http://dlmf.nist.gov/9.6.E5) @@ -122,7 +116,6 @@ _airybi(z::ComplexF16) = ComplexF16(_airybi(ComplexF32(z))) function _airybi(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} if ~isfinite(z) - isnan(z) && return z if abs(angle(z)) < 2π/3 return exp(z) else @@ -130,10 +123,7 @@ function _airybi(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} end end x, y = real(z), imag(z) - zabs = abs(z) - airy_large_argument_cutoff(z) && return airybi_large_argument(z) - airybi_power_series_cutoff(x, y) && return airybi_power_series(z) if x > zero(T) @@ -171,18 +161,14 @@ _airybiprime(z::ComplexF16) = ComplexF16(_airybiprime(ComplexF32(z))) function _airybiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} if ~isfinite(z) - isnan(z) && return z - if abs(angle(z)) < 2π/3 + if abs(angle(z)) < 2*T(π)/3 return exp(z) else return -1 / z end end x, y = real(z), imag(z) - zabs = abs(z) - airy_large_argument_cutoff(z) && return airybiprime_large_argument(z) - airybi_power_series_cutoff(x, y) && return airybiprime_power_series(z) if x > zero(T)