Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

RFC: Add error checking to Amos package-based Airy and Bessel methods. Correct # arguments in ccall to :zbiry. #4967

Merged
merged 1 commit into from Jan 21, 2014
Merged

RFC: Add error checking to Amos package-based Airy and Bessel methods. Correct # arguments in ccall to :zbiry. #4967

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
51 changes: 42 additions & 9 deletions base/math.jl
View file
Open in desktop
Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down
36 changes: 20 additions & 16 deletions doc/manual/mathematical-operations.rst
View file
Open in desktop
Expand Up @@ -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 <http://en.wikipedia.org/wiki/Dirichlet_eta_function>`_ at ``x``
``zeta(x)`` the `Riemann zeta function <http://en.wikipedia.org/wiki/Riemann_zeta_function>`_ at ``x``
``airy(x)``, ``airyai(x)`` the `Airy Ai function <http://en.wikipedia.org/wiki/Airy_function>`_ at ``x``
``airyprime(x)``, ``airyaiprime(x)`` the derivative of the Airy Ai function at ``x``
``airybi(x)`` the `Airy Bi function <http://en.wikipedia.org/wiki/Airy_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 <http://en.wikipedia.org/wiki/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 <http://en.wikipedia.org/wiki/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 <http://en.wikipedia.org/wiki/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 <http://en.wikipedia.org/wiki/Bessel_function>`_ of the first kind of order ``nu`` at ``z``
``besselk(nu,z)`` the modified `Bessel function <http://en.wikipedia.org/wiki/Bessel_function>`_ of the second kind of order ``nu`` at ``z``
|airylist| the `Airy Ai function <http://en.wikipedia.org/wiki/Airy_function>`_ at ``z``
|airyprimelist| the derivative of the Airy Ai function at ``z``
``airybi(z)``, ``airy(2,z)`` the `Airy Bi function <http://en.wikipedia.org/wiki/Airy_function>`_ at ``z``
``airybiprime(z)``, ``airy(3,z)`` the derivative of the Airy Bi function at ``z``
``besselj(nu,z)`` the `Bessel function <http://en.wikipedia.org/wiki/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 <http://en.wikipedia.org/wiki/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 <http://en.wikipedia.org/wiki/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 <http://en.wikipedia.org/wiki/Bessel_function>`_ of the first kind of order ``nu`` at ``z``
``besselk(nu,z)`` the modified `Bessel function <http://en.wikipedia.org/wiki/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)``


9 changes: 9 additions & 0 deletions test/math.jl
View file
Open in desktop
Expand Up @@ -49,20 +49,25 @@ 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
@test_approx_eq besselh(3,1,3) true_h133
@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
@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_approx_eq besseli(-3,-3) -true_i33
@test_throws besseli(1,1000)

# besselj
@test besselj(0,0) == 1
Expand All @@ -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
Expand All @@ -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.)
Expand All @@ -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
Expand Down