From 8077d116e22f12342bc392c54eedc470d9d5f2ef Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Wed, 17 Feb 2016 16:07:15 +0000 Subject: [PATCH] Don't silently convert `BigFloat`s in Bessel and Airy functions (fixes #11512). --- base/mpfr.jl | 1 + base/special/bessel.jl | 213 ++++++++++++++++++++++------------------- test/math.jl | 40 +++++++- 3 files changed, 150 insertions(+), 104 deletions(-) diff --git a/base/mpfr.jl b/base/mpfr.jl index 0dcc262a46c39..fefd138e34629 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -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() diff --git a/base/special/bessel.jl b/base/special/bessel.jl index f4248b70f93ae..778ef468e8b33 100644 --- a/base/special/bessel.jl +++ b/base/special/bessel.jl @@ -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) @@ -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)) @@ -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) @@ -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)) @@ -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) @@ -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 @@ -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) @@ -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 diff --git a/test/math.jl b/test/math.jl index c1bdc14a58d63..95aa40cbec123 100644 --- a/test/math.jl +++ b/test/math.jl @@ -314,7 +314,7 @@ end @test_throws Base.Math.AmosException airybi(200) @test_throws ArgumentError airy(5,one(Complex128)) z = 1.8 + 1.0im -for elty in [Complex64,Complex128, Complex{BigFloat}] +for elty in [Complex64,Complex128] @test_approx_eq airy(convert(elty,1.8)) 0.0470362168668458052247 z = convert(elty,z) @test_approx_eq airyx(z) airyx(0,z) @@ -324,6 +324,7 @@ for elty in [Complex64,Complex128, Complex{BigFloat}] @test_approx_eq airyx(3, z) airy(3, z) * exp(-abs(real(2/3 * z * sqrt(z)))) @test_throws ArgumentError airyx(5,z) end +@test_throws MethodError airy(complex(big(1.0))) # bessely0, bessely1, besselj0, besselj1 @test_approx_eq besselj0(Float32(2.0)) besselj0(Float64(2.0)) @@ -339,6 +340,11 @@ end @test_approx_eq bessely0(2.0 + im) bessely(0, 2.0 + im) @test_approx_eq bessely1(2.0 + im) bessely(1, 2.0 + im) +@test_throws MethodError besselj(1.2,big(1.0)) +@test_throws MethodError besselj(1,complex(big(1.0))) +@test_throws MethodError besseljx(1,big(1.0)) +@test_throws MethodError besseljx(1,complex(big(1.0))) + # besselh true_h133 = 0.30906272225525164362 - 0.53854161610503161800im @test_approx_eq besselh(3,1,3) true_h133 @@ -347,6 +353,10 @@ true_h133 = 0.30906272225525164362 - 0.53854161610503161800im @test_approx_eq besselh(-3,2,3) -conj(true_h133) @test_throws Base.Math.AmosException besselh(1,0) +@test_throws MethodError besselh(1,big(1.0)) +@test_throws MethodError besselh(1,complex(big(1.0))) +@test_throws MethodError besselhx(1,big(1.0)) +@test_throws MethodError besselhx(1,complex(big(1.0))) # besseli true_i33 = 0.95975362949600785698 @@ -357,6 +367,12 @@ true_i33 = 0.95975362949600785698 @test_throws Base.Math.AmosException besseli(1,1000) @test_throws DomainError besseli(0.4,-1.0) +@test_throws MethodError besseli(1,big(1.0)) +@test_throws MethodError besseli(1,complex(big(1.0))) +@test_throws MethodError besselix(1,big(1.0)) +@test_throws MethodError besselix(1,complex(big(1.0))) + + # besselj @test besselj(0,0) == 1 for i = 1:5 @@ -385,9 +401,8 @@ j43 = besselj(4,3.) @test_approx_eq besselj(3.2, 1.3+0.6im) 0.01135309305831220201 + 0.03927719044393515275im @test_approx_eq besselj(1, 3im) 3.953370217402609396im @test_approx_eq besselj(1.0,3im) besselj(1,3im) -@test besselj(big(1.0),3im) ≈ besselj(1,3im) -@test besselj(big(0.1), complex(-0.4)) ≈ 0.820421842809028916 + 0.266571215948350899im @test_throws Base.Math.AmosException besselj(20,1000im) +@test_throws MethodError besselj(big(1.0),3im) # besselk true_k33 = 0.12217037575718356792 @@ -401,6 +416,12 @@ true_k3m3 = -0.1221703757571835679 - 3.0151549516807985776im # issue #6564 @test besselk(1.0,0.0) == Inf +@test_throws MethodError besselk(1,big(1.0)) +@test_throws MethodError besselk(1,complex(big(1.0))) +@test_throws MethodError besselkx(1,big(1.0)) +@test_throws MethodError besselkx(1,complex(big(1.0))) + + # bessely y33 = bessely(3,3.) @test bessely(3,3) == y33 @@ -415,12 +436,23 @@ y33 = bessely(3,3.) @test_throws DomainError bessely(0.4,Float32(-1.0)) @test_throws DomainError bessely(1,Float32(-1.0)) +@test_throws MethodError bessely(1.2,big(1.0)) +@test_throws MethodError bessely(1,complex(big(1.0))) +@test_throws MethodError besselyx(1,big(1.0)) +@test_throws MethodError besselyx(1,complex(big(1.0))) + + #besselhx -for elty in [Complex64,Complex128, Complex{BigFloat}] +for elty in [Complex64,Complex128] z = convert(elty, 1.0 + 1.9im) @test_approx_eq besselhx(1.0, 1, z) convert(elty,-0.5949634147786144 - 0.18451272807835967im) end +@test_throws MethodError besselh(1,1,big(1.0)) +@test_throws MethodError besselh(1,1,complex(big(1.0))) +@test_throws MethodError besselhx(1,1,big(1.0)) +@test_throws MethodError besselhx(1,1,complex(big(1.0))) + # issue #6653 for f in (besselj,bessely,besseli,besselk,hankelh1,hankelh2) @test_approx_eq f(0,1) f(0,Complex128(1))