Question

I have a routine that can return an array of any FloatingPoint type (or really any numeric type but it's intended for floating points). It computes the convergents of a continued fraction. My problem is that I want the compiler to optimize it for each return type, e.g. Float64, BigFloat, etc. rather than using generic code.

function floatconvergents(a::Vector{Int}, n_quotients::Int, t::Type)
    r = Array(t, n_quotients)
    u = Array(t, n_quotients)
    v = Array(t, n_quotients)
    r[1],u[1],v[1] = a[1],0,1
    for k=2:n_quotients
        u[k] = 1 / (a[k] + u[k-1])
        r[k] = (a[k]*r[k-1] + v[k-1]) * u[k]
        v[k] = r[k-1] * u[k]
    end
    return r
end

The only solution that seems feasible to me is to pass the result array as an argument and make it a parametric function

function floatconvergents!{T<:FloatingPoint}(a::Vector{Int}, n_quotients::Int, r::Vector{T})

The only problem I see with this is that there's no reason for supplying the result array as an argument except to specify its type. It just doesn't seem like the "right" thing to do.

Was it helpful?

Solution

This will dispatch on the type, thus it will specialize the code:

function floatconvergents{T<:FloatingPoint}(a::Vector{Int}, n_quotients::Int, ::Type{T})
    r = Array(T, n_quotients)
    u = Array(T, n_quotients)
    v = Array(T, n_quotients)
    r[1],u[1],v[1] = a[1],0,1
    for k=2:n_quotients
        u[k] = 1 / (a[k] + u[k-1])
        r[k] = (a[k]*r[k-1] + v[k-1]) * u[k]
        v[k] = r[k-1] * u[k]
    end
    return r
end

It can be called, for example, by floatconvergents(some_vector, some_int, Float64).

Proof:

julia> f{T<:FloatingPoint}(::Type{T}) = T===Float64 ? 1 : 0
julia> code_llvm(f, (Type{Float64},))
define i64 @julia_f15644(%jl_value_t*) {
top:
  ret i64 1, !dbg !2042
}
julia> code_llvm(f, (Type{BigFloat},))
define i64 @julia_f15645(%jl_value_t*) {
top:
  ret i64 0, !dbg !2045
}
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top