Skip to content

Commit

Permalink
Bug fixing only
Browse files Browse the repository at this point in the history
  • Loading branch information
@simonp0420
simonp0420 committed Jan 13, 2014
1 parent 8ecfde5 commit eae2b84
Show file tree
Hide file tree
Showing 3 changed files with 71 additions and 25 deletions.
51 changes: 42 additions & 9 deletions base/math.jl
View file
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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

0 comments on commit eae2b84

Please sign in to comment.