Computing polynomials: outperforming C++

When you want to convince people that Julia is the best language for some task, some are easily swayed by the beautiful syntax, rank polymorphism, and flirtations with functional programming. In all honesty, those people need little convincing. Others want to see hard numbers!

Here's a very small demo that shows how Julia can outperform any traditional language: computing polynomials.

Suppose we have a function described by

where we have the coefficients stored in some dense vector.

struct Polynomial{T}
    c::Vector{T}
end

We can evaluate this function as follows:

function (f::Polynomial{T})(x::T) where T<:Number
    result = f.c[1]
    xpow = x
    for c in f.c[2:end-1]
        result += xpow * c
        xpow *= x
    end
    result + xpow * f.c[end]
end

This is very much how you would implement such a function in most other languages. In Julia however, we can do something much nicer: generate a function.

function expand(f::Polynomial{T}) where T<:Number
    :(function (x::$T)
        r = $(f.c[1])
        xp = x
        $((:(r += xp*$c; xp*=x) for c in f.c[2:end-1])...)
        r + xp * $(f.c[end])
    end)
end

Now, when we expand(f), say we have f = Polynomial(1.0, 0.5, 0.333, 0.25) the generated function will look as follows:

:(function (x::Float64,)
      r = 1.0
      xp = x
      begin
          r += xp * 0.5
          xp *= x
      end
      begin
          r += xp * 0.333
          xp *= x
      end
      r + xp * 0.25
  end)

To run this, we say f_unroll = eval(expand(f)). Not only does this run much faster than our previous loopy implementation, it even outperforms C++ (GCC at -O3) by more than a factor of two (depending on your machine). The reason is that we have no loop in our generated code, and more importantly, we are not accessing any memory. We have done something that a C++ programmer wouldn't dream of doing: write a code generator, compile the specialised code, and fly to the moon!

A completely literate and Documentered implementation of this demo can be found here.

Edit: Thanks to Nicos Pitsianis for point out to me that there is a more efficient method called Horner's method. This method applies to both Julia and C++ implementations, reducing the amount of multiplications needed. Using this method reduces the difference between Julia and C++ a little bit, but increasing the order of the polynomial grows it back again.