Skip to content

Commit

Permalink
Fix #6665
Browse files Browse the repository at this point in the history 
Provides airyx in addition to bessel[jyikh]x
  • Loading branch information
@jdrugo
jdrugo committed Apr 29, 2014
1 parent fd25446 commit d52ee34
Show file tree
Hide file tree
Showing 5 changed files with 197 additions and 47 deletions.
6 changes: 6 additions & 0 deletions base/exports.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -448,13 +448,19 @@ export
airybi,
airybiprime,
airyprime,
airyx,
besselh,
besselhx,
besseli,
besselix,
besselj,
besseljx,
besselj0,
besselj1,
besselk,
besselkx,
bessely,
besselyx,
bessely0,
bessely1,
beta,
Expand Down
179 changes: 134 additions & 45 deletions base/math.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,10 @@ export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan,
ceil, floor, trunc, round, significand,
lgamma, hypot, gamma, lfact, max, min, minmax, ldexp, frexp,
clamp, modf, ^, mod2pi,
airy, airyai, airyprime, airyaiprime, airybi, airybiprime,
airy, airyai, airyprime, airyaiprime, airybi, airybiprime, airyx,
besselj0, besselj1, besselj, bessely0, bessely1, bessely,
hankelh1, hankelh2, besseli, besselk, besselh,
besselhx, besselix, besseljx, besselkx, besselyx,
beta, lbeta, eta, zeta, polygamma, invdigamma, digamma, trigamma,
erfinv, erfcinv

Expand Down Expand Up @@ -451,40 +452,47 @@ end
let
const ai::Array{Float64,1} = Array(Float64,2)
const ae::Array{Int32,1} = Array(Int32,2)
global airy
function airy(k::Int, z::Complex128)
id = int32(k==1 || k==3)
if k == 0 || k == 1
global _airy, _biry
function _airy(z::Complex128, id::Int32, kode::Int32)
ccall((:zairy_,openspecfun), Void,
(Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
&real(z), &imag(z),
&id, &1,
&id, &kode,
pointer(ai,1), pointer(ai,2),
pointer(ae,1), pointer(ae,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
end
end
function _biry(z::Complex128, id::Int32, kode::Int32)
ccall((:zbiry_,openspecfun), Void,
(Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}),
&real(z), &imag(z),
&id, &1,
&id, &kode,
pointer(ai,1), pointer(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
end
end

function airy(k::Int, z::Complex128)
id = int32(k==1 || k==3)
if k == 0 || k == 1
return _airy(z, id, int32(1))
elseif k == 2 || k == 3
return _biry(z, id, int32(1))
else
error("invalid argument")
end
end
end

airy(z) = airy(0,z)
@vectorize_1arg Number airy
Expand All @@ -505,15 +513,35 @@ airy(k::Number, z::Complex64) = complex64(airy(k, complex128(z)))
airy(k::Number, z::Complex) = airy(convert(Int,k), complex128(z))
@vectorize_2arg Number airy

function airyx(k::Int, z::Complex128)
id = int32(k==1 || k==3)
if k == 0 || k == 1
return _airy(z, id, int32(2))
elseif k == 2 || k == 3
return _biry(z, id, int32(2))
else
error("invalid argument")
end
end

airyx(z) = airyx(0,z)
@vectorize_1arg Number airyx

airyx(k::Number, x::FloatingPoint) = oftype(x, real(airyx(k, complex(x))))
airyx(k::Number, x::Real) = airyx(k, float(x))
airyx(k::Number, z::Complex64) = complex64(airyx(k, complex128(z)))
airyx(k::Number, z::Complex) = airyx(convert(Int,k), complex128(z))
@vectorize_2arg Number airyx

const cy = Array(Float64,2)
const ae = Array(Int32,2)
const wrk = Array(Float64,2)

function _besselh(nu::Float64, k::Integer, z::Complex128)
function _besselh(nu::Float64, k::Int32, z::Complex128, kode::Int32)
ccall((:zbesh_,openspecfun), Void,
(Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
&real(z), &imag(z), &nu, &1, &k, &1,
&real(z), &imag(z), &nu, &kode, &k, &1,
pointer(cy,1), pointer(cy,2),
pointer(ae,1), pointer(ae,2))
if ae[2] == 0 || ae[2] == 3
Expand All @@ -523,11 +551,11 @@ function _besselh(nu::Float64, k::Integer, z::Complex128)
end
end

function _besseli(nu::Float64, z::Complex128)
function _besseli(nu::Float64, z::Complex128, kode::Int32)
ccall((:zbesi_,openspecfun), Void,
(Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
&real(z), &imag(z), &nu, &1, &1,
&real(z), &imag(z), &nu, &kode, &1,
pointer(cy,1), pointer(cy,2),
pointer(ae,1), pointer(ae,2))
if ae[2] == 0 || ae[2] == 3
Expand All @@ -537,11 +565,11 @@ function _besseli(nu::Float64, z::Complex128)
end
end

function _besselj(nu::Float64, z::Complex128)
function _besselj(nu::Float64, z::Complex128, kode::Int32)
ccall((:zbesj_,openspecfun), Void,
(Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
&real(z), &imag(z), &nu, &1, &1,
&real(z), &imag(z), &nu, &kode, &1,
pointer(cy,1), pointer(cy,2),
pointer(ae,1), pointer(ae,2))
if ae[2] == 0 || ae[2] == 3
Expand All @@ -551,11 +579,11 @@ function _besselj(nu::Float64, z::Complex128)
end
end

function _besselk(nu::Float64, z::Complex128)
function _besselk(nu::Float64, z::Complex128, kode::Int32)
ccall((:zbesk_,openspecfun), Void,
(Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
&real(z), &imag(z), &nu, &1, &1,
&real(z), &imag(z), &nu, &kode, &1,
pointer(cy,1), pointer(cy,2),
pointer(ae,1), pointer(ae,2))
if ae[2] == 0 || ae[2] == 3
Expand All @@ -565,12 +593,12 @@ function _besselk(nu::Float64, z::Complex128)
end
end

function _bessely(nu::Float64, z::Complex128)
function _bessely(nu::Float64, z::Complex128, kode::Int32)
ccall((:zbesy_,openspecfun), Void,
(Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32},
Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32},
Ptr{Float64}, Ptr{Float64}, Ptr{Int32}),
&real(z), &imag(z), &nu, &1, &1,
&real(z), &imag(z), &nu, &kode, &1,
pointer(cy,1), pointer(cy,2),
pointer(ae,1), pointer(wrk,1),
pointer(wrk,2), pointer(ae,2))
Expand All @@ -584,24 +612,40 @@ end
function besselh(nu::Float64, k::Integer, z::Complex128)
if nu < 0
s = (k == 1) ? 1 : -1
return _besselh(-nu, k, z) * complex(cospi(nu),-s*sinpi(nu))
return _besselh(-nu,int32(k),z,int32(1)) * complex(cospi(nu),-s*sinpi(nu))
end
return _besselh(nu, k, z)
return _besselh(nu,int32(k),z,int32(1))
end

function besselhx(nu::Float64, k::Integer, z::Complex128)
if nu < 0
s = (k == 1) ? 1 : -1
return _besselh(-nu,int32(k),z,int32(2)) * complex(cospi(nu),-s*sinpi(nu))
end
return _besselh(nu,int32(k),z,int32(2))
end

function besseli(nu::Float64, z::Complex128)
if nu < 0
return _besseli(-nu,z) - 2_besselk(-nu,z)*sinpi(nu)/pi
return _besseli(-nu,z,int32(1)) - 2_besselk(-nu,z,int32(1))*sinpi(nu)/pi
else
return _besseli(nu,z,int32(1))
end
end

function besselix(nu::Float64, z::Complex128)
if nu < 0
return _besseli(-nu,z,int32(2)) - 2_besselk(-nu,z,int32(2))*exp(-abs(real(z))-z)*sinpi(nu)/pi
else
return _besseli(nu, z)
return _besseli(nu,z,int32(2))
end
end

function besselj(nu::Float64, z::Complex128)
if nu < 0
return _besselj(-nu,z)cos(pi*nu) + _bessely(-nu,z)*sinpi(nu)
return _besselj(-nu,z,int32(1))*cospi(nu) + _bessely(-nu,z,int32(1))*sinpi(nu)
else
return _besselj(nu, z)
return _besselj(nu,z,int32(1))
end
end

Expand All @@ -613,13 +657,31 @@ besselj(nu::Integer, x::Float32) = typemin(Int32) <= nu <= typemax(Int32) ?
ccall((:jnf, libm), Float32, (Cint, Float32), nu, x) :
besselj(float64(nu), x)

besselk(nu::Float64, z::Complex128) = _besselk(abs(nu), z)
function besseljx(nu::Float64, z::Complex128)
if nu < 0
return _besselj(-nu,z,int32(2))*cospi(nu) + _bessely(-nu,z,int32(2))*sinpi(nu)
else
return _besselj(nu,z,int32(2))
end
end

besselk(nu::Float64, z::Complex128) = _besselk(abs(nu), z, int32(1))

besselkx(nu::Float64, z::Complex128) = _besselk(abs(nu), z, int32(2))

function bessely(nu::Float64, z::Complex128)
if nu < 0
return _bessely(-nu,z)*cospi(nu) - _besselj(-nu,z)*sinpi(nu)
return _bessely(-nu,z,int32(1))*cospi(nu) - _besselj(-nu,z,int32(1))*sinpi(nu)
else
return _bessely(nu, z)
return _bessely(nu,z,int32(1))
end
end

function besselyx(nu::Float64, z::Complex128)
if nu < 0
return _bessely(-nu,z,int32(2))*cospi(nu) - _besselj(-nu,z,int32(2))*sinpi(nu)
else
return _bessely(nu,z,int32(2))
end
end

Expand All @@ -629,16 +691,25 @@ besselh(nu::Real, k::Integer, z::Complex) = besselh(float64(nu), k, complex128(z
besselh(nu::Real, k::Integer, x::Real) = besselh(float64(nu), k, complex128(x))
@vectorize_2arg Number besselh

besseli(nu::Real, z::Complex64) = complex64(besseli(float64(nu), complex128(z)))
besseli(nu::Real, z::Complex) = besseli(float64(nu), complex128(z))
besseli(nu::Real, x::Integer) = besseli(nu, float64(x))
besselhx(nu, z) = besselhx(nu, 1, z)
besselhx(nu::Real, k::Integer, z::Complex64) = complex64(besselhx(float64(nu), k, complex128(z)))
besselhx(nu::Real, k::Integer, z::Complex) = besselhx(float64(nu), k, complex128(z))
besselhx(nu::Real, k::Integer, x::Real) = besselhx(float64(nu), k, complex128(x))
@vectorize_2arg Number besselhx

function besseli(nu::Real, x::FloatingPoint)
if x < 0 && !isinteger(nu)
throw(DomainError())
end
oftype(x, real(besseli(float64(nu), complex128(x))))
end
@vectorize_2arg Number besseli

function besselix(nu::Real, x::FloatingPoint)
if x < 0 && !isinteger(nu)
throw(DomainError())
end
oftype(x, real(besselix(float64(nu), complex128(x))))
end

function besselj(nu::FloatingPoint, x::FloatingPoint)
if isinteger(nu)
Expand All @@ -651,25 +722,27 @@ function besselj(nu::FloatingPoint, x::FloatingPoint)
oftype(x, real(besselj(float64(nu), complex128(x))))
end

besselj(nu::Real, z::Complex64) = complex64(besselj(float64(nu), complex128(z)))
besselj(nu::Real, z::Complex) = besselj(float64(nu), complex128(z))
besselj(nu::Real, x::Integer) = besselj(nu, float(x))
@vectorize_2arg Number besselj
function besseljx(nu::Real, x::FloatingPoint)
if x < 0 && !isinteger(nu)
throw(DomainError())
end
oftype(x, real(besseljx(float64(nu), complex128(x))))
end

besselk(nu::Real, z::Complex64) = complex64(besselk(float64(nu), complex128(z)))
besselk(nu::Real, z::Complex) = besselk(float64(nu), complex128(z))
besselk(nu::Real, x::Integer) = besselk(nu, float64(x))
function besselk(nu::Real, x::FloatingPoint)
if x < 0
throw(DomainError())
end
oftype(x, real(besselk(float64(nu), complex128(x))))
end
@vectorize_2arg Number besselk

bessely(nu::Real, z::Complex64) = complex64(bessely(float64(nu), complex128(z)))
bessely(nu::Real, z::Complex) = bessely(float64(nu), complex128(z))
bessely(nu::Real, x::Integer) = bessely(nu, float64(x))
function besselkx(nu::Real, x::FloatingPoint)
if x < 0
throw(DomainError())
end
oftype(x, real(besselkx(float64(nu), complex128(x))))
end

function bessely(nu::Real, x::FloatingPoint)
if x < 0
throw(DomainError())
Expand All @@ -691,7 +764,23 @@ function bessely(nu::Integer, x::Float32)
end
return ccall((:ynf, libm), Float32, (Cint, Float32), nu, x)
end
@vectorize_2arg Number bessely

function besselyx(nu::Real, x::FloatingPoint)
if x < 0
throw(DomainError())
end
oftype(x, real(besselyx(float64(nu), complex128(x))))
end

for f in ("i", "j", "k", "y"), g in ("", "x")
bfn = symbol(string("bessel", f, g))
@eval begin
$bfn(nu::Real, z::Complex64) = complex64($bfn(float64(nu), complex128(z)))
$bfn(nu::Real, z::Complex) = $bfn(float64(nu), complex128(z))
$bfn(nu::Real, x::Integer) = $bfn(nu, float64(x))
@vectorize_2arg Number $bfn
end
end

hankelh1(nu, z) = besselh(nu, 1, z)
@vectorize_2arg Number hankelh1
Expand Down
11 changes: 9 additions & 2 deletions doc/manual/mathematical-operations.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -452,17 +452,24 @@ Function Description
|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``
``airyx(0,z)``, ``airyx(1,z)`` ``airy(0,z)``, ``airy(1,z)`` scaled by ``exp(2/3 * z * sqrt(z))``
``airyx(2,z)``, ``airyx(3,z)`` ``airy(2,z)``, ``airy(3,z)`` scaled by ``exp(-abs(real(2/3 * z * sqrt(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)``
``besseljx(nu,z)`` ``besselj(nu,z) * exp(-abs(imag(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)``
``besselyx(nu,z)`` ``bessely(nu,z) * exp(-abs(imag(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)``
``hankelh1(nu,z)`` ``besselh(nu,1,z)``
``hankelh2(nu,z)`` ``besselh(nu,2,z)``
``besselhx(nu,k,z)`` ``besselh(nu,k,z) * exp(s * z * im)`` with ``s=-1`` for ``k==1``, ``s=1`` for ``k==2``
``besseli(nu,z)`` the modified `Bessel function <http://en.wikipedia.org/wiki/Bessel_function>`_ of the first kind of order ``nu`` at ``z``
``besselix(nu,z)`` ``besseli(nu,z) * exp(-abs(real(z)))``
``besselk(nu,z)`` the modified `Bessel function <http://en.wikipedia.org/wiki/Bessel_function>`_ of the second kind of order ``nu`` at ``z``
``besselkx(nu,z)`` ``besselk(nu,z) * exp(z)``
====================================== ==============================================================================

.. |airylist| replace:: ``airy(z)``, ``airyai(z)``, ``airy(0,z)``
Expand Down
Loading

0 comments on commit d52ee34

Please sign in to comment.