Test type-stability of functions

[giordano]

Use @inferred to test type-stability of functions.

• besselj and bessely 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.

[giordano]

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

[musm]

Excellent catch @giordano

[musm]

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

[giordano]

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

[musm]

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

[giordano]

So we're saying the same thing 🙂

[musm]

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.

[musm]

Anyways I think the discussion got sidetracked, since these are besides the point of this PR.

[cossio]

bump

[ViralBShah]

@giordano Should we get this merged? If so, would it be possible for you to rebase?

[andreasnoack]

Superseded by #549, which carries this idea forward — thanks @giordano, and sorry it took nine years.

Rather than rebasing, #549 adds a dedicated test/type_stability.jl, because test/runtests.jl has since been split into per-topic files, so the diff here no longer applies to anything. Two things changed relative to this PR:

• it asserts @inferred(f(x)) isa T instead of a bare @inferred. Of the 43 methods that currently get this wrong, only 14 fail a bare @inferred; the other 29 are inferred perfectly and just return a Float64, so the weaker check would have missed them.
• the Bessel cases are marked @test_broken instead of blocking the PR, so the sweep can land now and the list can be worked down one entry at a time.

The instabilities listed in the description here are resolved for digamma, trigamma, invdigamma, polygamma, eta and zeta — all of them infer Float32 from Float32 today. besselj and bessely are still broken, exactly as diagnosed in the discussion above: the missing else branch you pointed at is why besselj(2, 1.0f0) is Union{Float32,Float64}, and it is also why besselj(2, Float16(1)) recurses until the stack overflows (#547).

---

Comment written by Claude Code on behalf of @andreasnoack.