SpecialFunctions.jl
SpecialFunctions.jl copied to clipboard
Test type-stability of functions
Use @inferred
to test type-stability of functions.
-
besselj
andbessely
are unstable in Julia 0.6 and 0.7. -
digamma
,trigamma
,invdigamma
,polygamma
,eta
,zeta
are unstable in Julia 0.7 only.
Type-instability of digamma
was already reported in issue #42.
I believe the instabilities in besselj
and bessely
are caused by the conditionals
if typemin(Cint) <= nu <= typemax(Cint)
for which there is no else
clause, so the return type is probably undetermined (actually, it should be an error).
The tests pass in Julia 0.6 with the following patch:
diff --git a/src/bessel.jl b/src/bessel.jl
index 958ac37..4f52082 100644
--- a/src/bessel.jl
+++ b/src/bessel.jl
@@ -388,6 +388,8 @@ function besselj(nu::Real, x::AbstractFloat)
if isinteger(nu)
if typemin(Cint) <= nu <= typemax(Cint)
return besselj(Cint(nu), x)
+ else
+ error()
end
elseif x < 0
throw(DomainError(x, "`x` must be nonnegative."))
@@ -443,8 +445,12 @@ Bessel function of the second kind of order `nu`, ``Y_\\nu(x)``.
function bessely(nu::Real, x::AbstractFloat)
if x < 0
throw(DomainError(x, "`x` must be nonnegative."))
- elseif isinteger(nu) && typemin(Cint) <= nu <= typemax(Cint)
- return bessely(Cint(nu), x)
+ elseif isinteger(nu)
+ if typemin(Cint) <= nu <= typemax(Cint)
+ return bessely(Cint(nu), x)
+ else
+ error()
+ end
end
real(bessely(float(nu), complex(x)))
end
Well, in place of error()
there should be a more meaningful error (AmosException
?)
Excellent catch @giordano
Hmm actually I think the problem is a little different. The problem is not with the conditionals but rather
bessely(nu::Real, z::Complex) in SpecialFunctions at C:\Users\Mus\.julia\v0.7\SpecialFunctions\src\bessel.jl:470
always return a Float64
, which throws off type stability of bessely(nu::Real, x::AbstractFloat)
and this is due to the promotion Tf = promote_type(float(typeof(nu)),float(typeof(real(z))))
Uhm, not sure I understand what you mean:
julia> besselj(-3f0, complex(3f0))
-0.30906272f0 - 3.7849267f-17im
julia> bessely(-3f0, complex(3f0))
0.5385416f0 + 0.0f0im
are correctly Complex{Float32}
.
Maybe the issue is how those methods are called? float(nu)
looks a bit suspicious: when the conditional is false but nu::Int
then float(nu)
is always Float64
, whatever the type of x
Consider
julia> x = 2f0
2.0f0
julia> nu = 2.0
2.0
julia> bessely(nu,x)
-0.6174081f0
This calls bessely(Cint(nu), x)
which has the correct return type (within the isinteger(nu) && typemin(Cint) <= nu <= typemax(Cint)
) branch. However
julia> real(bessely(float(nu), complex(x)))
-0.6174081041906828
which is called on line 449 and this causes the type instability
So we're saying the same thing 🙂
Basically, but I don't think the patch you propose actually fixes the type instability; may need to stick oftype(x, real(bessely(float(nu), complex(x))))
on line 449, which fixes it at least for Float32
inputs, but maybe there is a cleaner fix.
Anyways I think the discussion got sidetracked, since these are besides the point of this PR.
bump
@giordano Should we get this merged? If so, would it be possible for you to rebase?