diff --git a/base/math.jl b/base/math.jl index 6f4ea70caaf70..67187a5997edf 100644 --- a/base/math.jl +++ b/base/math.jl @@ -432,6 +432,11 @@ for jy in ("j","y"), nu in (0,1) end end + +type AmosException <: Exception + info::Int32 +end + let const ai::Array{Float64,1} = Array(Float64,2) const ae::Array{Int32,1} = Array(Int32,2) @@ -446,16 +451,24 @@ function airy(k::Int, z::Complex128) &id, &1, pointer(ai,1), pointer(ai,2), pointer(ae,1), pointer(ae,2)) - return complex(ai[1],ai[2]) + if ae[2] == 0 || ae[2] == 3 + return complex(ai[1],ai[2]) + else + throw(AmosException(ae[2])) + end elseif k == 2 || k == 3 ccall((:zbiry_,openspecfun), Void, (Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), + Ptr{Float64}, Ptr{Float64}, Ptr{Int32}), &real(z), &imag(z), &id, &1, pointer(ai,1), pointer(ai,2), - pointer(ae,1), pointer(ae,2)) - return complex(ai[1],ai[2]) + pointer(ae,1)) + if ae[1] == 0 || ae[1] == 3 # ignore underflow + return complex(ai[1],ai[2]) + else + throw(AmosException(ae[2])) + end else error("invalid argument") end @@ -492,7 +505,11 @@ function _besselh(nu::Float64, k::Integer, z::Complex128) &real(z), &imag(z), &nu, &1, &k, &1, pointer(cy,1), pointer(cy,2), pointer(ae,1), pointer(ae,2)) - return complex(cy[1],cy[2]) + if ae[2] == 0 || ae[2] == 3 + return complex(cy[1],cy[2]) + else + throw(AmosException(ae[2])) + end end function _besseli(nu::Float64, z::Complex128) @@ -502,7 +519,11 @@ function _besseli(nu::Float64, z::Complex128) &real(z), &imag(z), &nu, &1, &1, pointer(cy,1), pointer(cy,2), pointer(ae,1), pointer(ae,2)) - return complex(cy[1],cy[2]) + if ae[2] == 0 || ae[2] == 3 + return complex(cy[1],cy[2]) + else + throw(AmosException(ae[2])) + end end function _besselj(nu::Float64, z::Complex128) @@ -512,7 +533,11 @@ function _besselj(nu::Float64, z::Complex128) &real(z), &imag(z), &nu, &1, &1, pointer(cy,1), pointer(cy,2), pointer(ae,1), pointer(ae,2)) - return complex(cy[1],cy[2]) + if ae[2] == 0 || ae[2] == 3 + return complex(cy[1],cy[2]) + else + throw(AmosException(ae[2])) + end end function _besselk(nu::Float64, z::Complex128) @@ -522,7 +547,11 @@ function _besselk(nu::Float64, z::Complex128) &real(z), &imag(z), &nu, &1, &1, pointer(cy,1), pointer(cy,2), pointer(ae,1), pointer(ae,2)) - return complex(cy[1],cy[2]) + if ae[2] == 0 || ae[2] == 3 + return complex(cy[1],cy[2]) + else + throw(AmosException(ae[2])) + end end function _bessely(nu::Float64, z::Complex128) @@ -534,7 +563,11 @@ function _bessely(nu::Float64, z::Complex128) pointer(cy,1), pointer(cy,2), pointer(ae,1), pointer(wrk,1), pointer(wrk,2), pointer(ae,2)) - return complex(cy[1],cy[2]) + if ae[2] == 0 || ae[2] == 3 + return complex(cy[1],cy[2]) + else + throw(AmosException(ae[2])) + end end function besselh(nu::Float64, k::Integer, z::Complex128) diff --git a/doc/manual/mathematical-operations.rst b/doc/manual/mathematical-operations.rst index cb905d48c0c47..669dc8c3821ee 100644 --- a/doc/manual/mathematical-operations.rst +++ b/doc/manual/mathematical-operations.rst @@ -447,20 +447,24 @@ Function Description ``lbeta(x,y)`` accurate ``log(beta(x,y))`` for large ``x`` or ``y`` ``eta(x)`` the `Dirichlet eta function `_ at ``x`` ``zeta(x)`` the `Riemann zeta function `_ at ``x`` -``airy(x)``, ``airyai(x)`` the `Airy Ai function `_ at ``x`` -``airyprime(x)``, ``airyaiprime(x)`` the derivative of the Airy Ai function at ``x`` -``airybi(x)`` the `Airy Bi function `_ at ``x`` -``airybiprime(x)`` the derivative of the Airy Bi function at ``x`` -``airy(k,x)`` the ``k``-th derivative of the Airy Ai function at ``x`` -``besselj(nu,z)`` the `Bessel function `_ of the first kind of order ``nu`` at ``z`` -``besselj0(z)`` ``besselj(0,z)`` -``besselj1(z)`` ``besselj(1,z)`` -``bessely(nu,z)`` the `Bessel function `_ of the second kind of order ``nu`` at ``z`` -``bessely0(z)`` ``bessely(0,z)`` -``bessely1(z)`` ``bessely(1,z)`` -``besselh(nu,k,z)`` the `Bessel function `_ of the third kind (a.k.a. Hankel function) of order ``nu`` at ``z``; ``k`` must be either ``1`` or ``2`` -``hankelh1(nu,z)`` ``besselh(nu, 1, z)`` -``hankelh2(nu,z)`` ``besselh(nu, 2, z)`` -``besseli(nu,z)`` the modified `Bessel function `_ of the first kind of order ``nu`` at ``z`` -``besselk(nu,z)`` the modified `Bessel function `_ of the second kind of order ``nu`` at ``z`` +|airylist| the `Airy Ai function `_ at ``z`` +|airyprimelist| the derivative of the Airy Ai function at ``z`` +``airybi(z)``, ``airy(2,z)`` the `Airy Bi function `_ at ``z`` +``airybiprime(z)``, ``airy(3,z)`` the derivative of the Airy Bi function at ``z`` +``besselj(nu,z)`` the `Bessel function `_ of the first kind of order ``nu`` at ``z`` +``besselj0(z)`` ``besselj(0,z)`` +``besselj1(z)`` ``besselj(1,z)`` +``bessely(nu,z)`` the `Bessel function `_ of the second kind of order ``nu`` at ``z`` +``bessely0(z)`` ``bessely(0,z)`` +``bessely1(z)`` ``bessely(1,z)`` +``besselh(nu,k,z)`` the `Bessel function `_ of the third kind (a.k.a. Hankel function) of order ``nu`` at ``z``; ``k`` must be either ``1`` or ``2`` +``hankelh1(nu,z)`` ``besselh(nu, 1, z)`` +``hankelh2(nu,z)`` ``besselh(nu, 2, z)`` +``besseli(nu,z)`` the modified `Bessel function `_ of the first kind of order ``nu`` at ``z`` +``besselk(nu,z)`` the modified `Bessel function `_ of the second kind of order ``nu`` at ``z`` ====================================== ============================================================================== + +.. |airylist| replace:: ``airy(z)``, ``airyai(z)``, ``airy(0,z)`` +.. |airyprimelist| replace:: ``airyprime(z)``, ``airyaiprime(z)``, ``airy(1,z)`` + + diff --git a/test/math.jl b/test/math.jl index a64a67e0f3fa2..3982711a10c4d 100644 --- a/test/math.jl +++ b/test/math.jl @@ -49,6 +49,8 @@ end @test_approx_eq airyprime(1.8) -0.0685247801186109345638 @test_approx_eq airybi(1.8) 2.595869356743906290060 @test_approx_eq airybiprime(1.8) 2.98554005084659907283 +@test_throws airy(200im) +@test_throws airybi(200) # besselh true_h133 = 0.30906272225525164362 - 0.53854161610503161800im @@ -56,6 +58,8 @@ true_h133 = 0.30906272225525164362 - 0.53854161610503161800im @test_approx_eq besselh(-3,1,3) -true_h133 @test_approx_eq besselh(3,2,3) conj(true_h133) @test_approx_eq besselh(-3,2,3) -conj(true_h133) +@test_throws besselh(1,0) + # besseli true_i33 = 0.95975362949600785698 @@ -63,6 +67,7 @@ true_i33 = 0.95975362949600785698 @test_approx_eq besseli(-3,3) true_i33 @test_approx_eq besseli(3,-3) -true_i33 @test_approx_eq besseli(-3,-3) -true_i33 +@test_throws besseli(1,1000) # besselj @test besselj(0,0) == 1 @@ -89,6 +94,7 @@ j43 = besselj(4,3.) @test_approx_eq besselj(0.1, complex(-0.4)) 0.820421842809028916 + 0.266571215948350899im @test_approx_eq besselj(3.2, 1.3+0.6im) 0.01135309305831220201 + 0.03927719044393515275im @test_approx_eq besselj(1, 3im) 3.953370217402609396im +@test_throws besselj(20,1000im) # besselk true_k33 = 0.12217037575718356792 @@ -98,6 +104,8 @@ true_k3m3 = -0.1221703757571835679 - 3.0151549516807985776im @test_throws besselk(3,-3) @test_approx_eq besselk(3,complex(-3)) true_k3m3 @test_approx_eq besselk(-3,complex(-3)) true_k3m3 +@test_throws besselk(200,0.01) + # bessely y33 = bessely(3,3.) @@ -106,6 +114,7 @@ y33 = bessely(3,3.) @test_approx_eq y33 -0.53854161610503161800 @test_throws bessely(3,-3) @test_approx_eq bessely(3,complex(-3)) 0.53854161610503161800 - 0.61812544451050328724im +@test_throws bessely(200.5,0.1) # beta, lbeta @test_approx_eq beta(3/2,7/2) 5π/128