How to return self recursive function?

Let’s benchmark before we guess:

Code (above + some fixes and minor changes)
using BenchmarkTools

struct RecSeq{N,T,U}
    a    :: U
    init :: NTuple{N,T}
end

# This defines the behavior of `RecSeq` as a callable object
# (previously, it was the closure body)
function (f::RecSeq{N})(n) where {N}
    n < N && return f.init[n+1]
    u = ntuple(Val(N)) do k # the function object is naturally named
        f(n-k)              # which makes it easy for it to recursively
    end                     # "call" itself
    f.a(u...)
end

fib_struct = RecSeq((uₙ₋₁, uₙ₋₂) -> uₙ₋₁ + uₙ₋₂, (0,1))

function recseq_self(a::T, init) where{T<:Function}
    let N=length(init), init=init
        return n-> begin
            n < N && return init[n+1]
            return a(var"#self#", n)
        end
    end
end

fib_self = recseq_self((f,n) -> f(n-1)+f(n-2), (0,1))

function recseq_functional(a::T, init) where{T<:Function}
    let N=length(init), init=init
        recfun = (self, n)-> begin
            n < N && return init[n+1]
            a(n->self(self,n), n)
        end
        return n->recfun(recfun, n) # essentially Y Combinator applied to recfun
    end
end

fib_functional = recseq_self((f,n) -> f(n-1)+f(n-2), (0,1))

@show fib_functional.(0:10) == fib_struct.(0:10)
@show fib_functional.(0:10) == fib_self.(0:10)
const testN = 20
println("Struct-based")
@btime $(fib_struct)($testN)
println("Self-based")
@btime $(fib_self)($testN)
println("Y-combinator-based")
@btime $(fib_functional)($testN)

Output (tests with N=20):

Struct-based
  71.939 μs (0 allocations: 0 bytes)
Self-based
  42.185 μs (0 allocations: 0 bytes)
Y-combinator-based
  42.395 μs (0 allocations: 0 bytes)

So actually the struct-based approach appears to be slowest, which I did not expect.
Also I would have thought that the “Y-Combinator” variant is slower than the “self”-based solution but this is only a very minor distance. This tells us that function calls are not a major contributor to the total runtime, I think.

I want to add that this is a very simple case where Julia was able to properly infer all types. Once that fails programs get much slower. As a rule of thumb, recursive programs are harder to infer than non-recursive ones. So for complex code, a more direct style might lead to faster code due to imperfect inference. (recent thread on inference problems in a recursive context: Failure to optimize due to type recursion)

1 Like