Skip to content

Commit

Permalink
Merge pull request #15119 from JuliaLang/sb/bessel
Browse files Browse the repository at this point in the history 
Don't silently convert `BigFloat`s in Bessel and Airy functions
  • Loading branch information
@simonbyrne
simonbyrne committed Feb 26, 2016
2 parents 009cede + 8077d11 commit 5e44ea2
Show file tree
Hide file tree
Showing 3 changed files with 150 additions and 104 deletions.
1 change: 1 addition & 0 deletions base/mpfr.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -465,6 +465,7 @@ function airyai(x::BigFloat)
ccall((:mpfr_ai, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
return z
end
airy(x::BigFloat) = airyai(x)

function ldexp(x::BigFloat, n::Clong)
z = BigFloat()
Expand Down
213 changes: 113 additions & 100 deletions base/special/bessel.jl
View file
Original file line number Diff line number Diff line change
@@ -1,32 +1,11 @@
# This file is a part of Julia. License is MIT: http://julialang.org/license

for jy in ("j","y"), nu in (0,1)
jynu = Expr(:quote, symbol(jy,nu))
jynuf = Expr(:quote, symbol(jy,nu,"f"))
bjynu = symbol("bessel",jy,nu)
if jy == "y"
@eval begin
$bjynu(x::Float64) = nan_dom_err(ccall(($jynu,libm), Float64, (Float64,), x), x)
$bjynu(x::Float32) = nan_dom_err(ccall(($jynuf,libm), Float32, (Float32,), x), x)
end
else
@eval begin
$bjynu(x::Float64) = ccall(($jynu,libm), Float64, (Float64,), x)
$bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x)
end
end
@eval begin
$bjynu(x::Real) = $bjynu(float(x))
$bjynu(x::Complex) = $(symbol("bessel",jy))($nu,x)
@vectorize_1arg Number $bjynu
end
end


type AmosException <: Exception
info::Int32
end

## Airy functions

let
const ai::Array{Float64,1} = Array(Float64,2)
const ae::Array{Int32,1} = Array(Int32,2)
Expand Down Expand Up @@ -61,7 +40,7 @@ let
end
end

function airy(k::Int, z::Complex128)
function airy(k::Integer, z::Complex128)
id = Int32(k==1 || k==3)
if k == 0 || k == 1
return _airy(z, id, Int32(1))
Expand All @@ -72,8 +51,6 @@ function airy(k::Int, z::Complex128)
end
end

airy(z) = airy(0,z)
@vectorize_1arg Number airy
airyprime(z) = airy(1,z)
@vectorize_1arg Number airyprime
airyai(z) = airy(0,z)
Expand All @@ -85,13 +62,7 @@ airybi(z) = airy(2,z)
airybiprime(z) = airy(3,z)
@vectorize_1arg Number airybiprime

airy(k::Number, x::AbstractFloat) = oftype(x, real(airy(k, complex(x))))
airy(k::Number, x::Real) = airy(k, float(x))
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)
function airyx(k::Integer, z::Complex128)
id = Int32(k==1 || k==3)
if k == 0 || k == 1
return _airy(z, id, Int32(2))
Expand All @@ -102,14 +73,48 @@ function airyx(k::Int, z::Complex128)
end
end

airyx(z) = airyx(0,z)
@vectorize_1arg Number airyx
for afn in (:airy,:airyx)
@eval begin
$afn(k::Integer, z::Complex) = $afn(k, float(z))
$afn{T<:AbstractFloat}(k::Integer, z::Complex{T}) = throw(MethodError($afn,(k,z)))
$afn(k::Integer, z::Complex64) = Complex64($afn(k, Complex128(z)))

$afn(k::Integer, x::Real) = $afn(k, float(x))
$afn(k::Integer, x::AbstractFloat) = real($afn(k, complex(x)))

$afn(z) = $afn(0,z)
@vectorize_1arg Number $afn
@vectorize_2arg Number $afn
end
end

## Bessel functions

# besselj0, besselj1, bessely0, bessely1
for jy in ("j","y"), nu in (0,1)
jynu = Expr(:quote, symbol(jy,nu))
jynuf = Expr(:quote, symbol(jy,nu,"f"))
bjynu = symbol("bessel",jy,nu)
if jy == "y"
@eval begin
$bjynu(x::Float64) = nan_dom_err(ccall(($jynu,libm), Float64, (Float64,), x), x)
$bjynu(x::Float32) = nan_dom_err(ccall(($jynuf,libm), Float32, (Float32,), x), x)
end
else
@eval begin
$bjynu(x::Float64) = ccall(($jynu,libm), Float64, (Float64,), x)
$bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x)
end
end
@eval begin
$bjynu(x::Real) = $bjynu(float(x))
$bjynu(x::Complex) = $(symbol("bessel",jy))($nu,x)
@vectorize_1arg Number $bjynu
end
end



airyx(k::Number, x::AbstractFloat) = 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)
Expand Down Expand Up @@ -227,13 +232,9 @@ function besselj(nu::Float64, z::Complex128)
end
end

besselj(nu::Integer, x::AbstractFloat) = typemin(Int32) <= nu <= typemax(Int32) ?
oftype(x, ccall((:jn, libm), Float64, (Cint, Float64), nu, x)) :
besselj(Float64(nu), x)
besselj(nu::Cint, x::Float64) = ccall((:jn, libm), Float64, (Cint, Float64), nu, x)
besselj(nu::Cint, x::Float32) = ccall((:jn, libm), Float32, (Cint, Float32), nu, x)

besselj(nu::Integer, x::Float32) = typemin(Int32) <= nu <= typemax(Int32) ?
ccall((:jnf, libm), Float32, (Cint, Float32), nu, x) :
besselj(Float64(nu), x)

function besseljx(nu::Float64, z::Complex128)
if nu < 0
Expand All @@ -247,6 +248,19 @@ 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::Cint, x::Float64)
if x < 0
throw(DomainError())
end
ccall((:yn, libm), Float64, (Cint, Float64), nu, x)
end
function bessely(nu::Cint, x::Float32)
if x < 0
throw(DomainError())
end
ccall((:ynf, libm), Float32, (Cint, Float32), nu, x)
end

function bessely(nu::Float64, z::Complex128)
if nu < 0
return _bessely(-nu,z,Int32(1))*cospi(nu) - _besselj(-nu,z,Int32(1))*sinpi(nu)
Expand All @@ -263,115 +277,114 @@ function besselyx(nu::Float64, z::Complex128)
end
end

besselh(nu, z) = besselh(nu, 1, z)
besselh(nu::Real, k::Integer, z::Complex64) = Complex64(besselh(Float64(nu), k, Complex128(z)))
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

hankelh1(nu, z) = besselh(nu, 1, z)
@vectorize_2arg Number hankelh1

hankelh2(nu, z) = besselh(nu, 2, z)
@vectorize_2arg Number hankelh2

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))

hankelh1x(nu, z) = besselhx(nu, 1, z)
@vectorize_2arg Number hankelh1x

hankelh2x(nu, z) = besselhx(nu, 2, z)
@vectorize_2arg Number hankelh2x

function besseli(nu::Real, x::AbstractFloat)
if x < 0 && !isinteger(nu)
throw(DomainError())
end
oftype(x, real(besseli(Float64(nu), Complex128(x))))
real(besseli(float(nu), complex(x)))
end

function besselix(nu::Real, x::AbstractFloat)
if x < 0 && !isinteger(nu)
throw(DomainError())
end
oftype(x, real(besselix(Float64(nu), Complex128(x))))
real(besselix(float(nu), complex(x)))
end

function besselj(nu::AbstractFloat, x::AbstractFloat)
function besselj(nu::Real, x::AbstractFloat)
if isinteger(nu)
if typemin(Int32) <= nu <= typemax(Int32)
return besselj(Int(nu), x)
if typemin(Cint) <= nu <= typemax(Cint)
return besselj(Cint(nu), x)
end
elseif x < 0
throw(DomainError())
end
oftype(x, real(besselj(Float64(nu), Complex128(x))))
real(besselj(float(nu), complex(x)))
end

function besseljx(nu::Real, x::AbstractFloat)
if x < 0 && !isinteger(nu)
throw(DomainError())
end
oftype(x, real(besseljx(Float64(nu), Complex128(x))))
real(besseljx(float(nu), complex(x)))
end

function besselk(nu::Real, x::AbstractFloat)
if x < 0
throw(DomainError())
end
if x == 0
elseif x == 0
return oftype(x, Inf)
end
oftype(x, real(besselk(Float64(nu), Complex128(x))))
real(besselk(float(nu), complex(x)))
end

function besselkx(nu::Real, x::AbstractFloat)
if x < 0
throw(DomainError())
end
if x == 0
elseif x == 0
return oftype(x, Inf)
end
oftype(x, real(besselkx(Float64(nu), Complex128(x))))
real(besselkx(float(nu), complex(x)))
end

function bessely(nu::Real, x::AbstractFloat)
if x < 0
throw(DomainError())
elseif isinteger(nu) && typemin(Cint) <= nu <= typemax(Cint)
return bessely(Cint(nu), x)
end
if isinteger(nu) && typemin(Int32) <= nu <= typemax(Int32)
return bessely(Int(nu), x)
end
oftype(x, real(bessely(Float64(nu), Complex128(x))))
end
function bessely(nu::Integer, x::AbstractFloat)
if x < 0
throw(DomainError())
end
return oftype(x, ccall((:yn, libm), Float64, (Cint, Float64), nu, x))
end
function bessely(nu::Integer, x::Float32)
if x < 0
throw(DomainError())
end
return ccall((:ynf, libm), Float32, (Cint, Float32), nu, x)
real(bessely(float(nu), complex(x)))
end

function besselyx(nu::Real, x::AbstractFloat)
if x < 0
throw(DomainError())
end
oftype(x, real(besselyx(Float64(nu), Complex128(x))))
real(besselyx(float(nu), complex(x)))
end

for f in ("i", "ix", "j", "jx", "k", "kx", "y", "yx")
bfn = symbol("bessel", f)
@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))
$bfn(nu::Real, x::Real) = $bfn(nu, float(x))
function $bfn(nu::Real, z::Complex)
Tf = promote_type(float(typeof(nu)),float(typeof(real(z))))
$bfn(Tf(nu), Complex{Tf}(z))
end
$bfn{T<:AbstractFloat}(k::T, z::Complex{T}) = throw(MethodError($bfn,(k,z)))
$bfn(nu::Float32, x::Complex64) = Complex64($bfn(Float64(nu), Complex128(x)))
@vectorize_2arg Number $bfn
end
end


for bfn in (:besselh, :besselhx)
@eval begin
$bfn(nu, z) = $bfn(nu, 1, z)
$bfn(nu::Real, k::Integer, x::Real) = $bfn(nu, k, float(x))
$bfn(nu::Real, k::Integer, x::AbstractFloat) = $bfn(float(nu), k, complex(x))

function $bfn(nu::Real, k::Integer, z::Complex)
Tf = promote_type(float(typeof(nu)),float(typeof(real(z))))
$bfn(Tf(nu), k, Complex{Tf}(z))
end

$bfn{T<:AbstractFloat}(nu::T, k::Integer, z::Complex{T}) = throw(MethodError($bfn,(nu,k,z)))
$bfn(nu::Float32, k::Integer, x::Complex64) = Complex64($bfn(Float64(nu), k, Complex128(x)))
@vectorize_2arg Number $bfn
end
end


hankelh1(nu, z) = besselh(nu, 1, z)
@vectorize_2arg Number hankelh1

hankelh2(nu, z) = besselh(nu, 2, z)
@vectorize_2arg Number hankelh2

hankelh1x(nu, z) = besselhx(nu, 1, z)
@vectorize_2arg Number hankelh1x

hankelh2x(nu, z) = besselhx(nu, 2, z)
@vectorize_2arg Number hankelh2x
Loading

0 comments on commit 5e44ea2

Please sign in to comment.