# I have a much faster version of evalpoly. Why is it faster?

**URL:** <https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899>\
**Category:** Performance\
**Created:** [April 23, 2022, 3:45pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899 "2022-04-23T15:45:34Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![Luapulu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luapulu/32/17669_2.png) [@Luapulu](https://discourse.julialang.org/u/Luapulu)\
**Post date:** [April 23, 2022, 3:45pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/1 "2022-04-23T15:45:34Z")

</div>

I was playing around with different ways to implement the `evalpoly` function.

**TLDR**

- My best version is about twice as fast as the version currently in base and works with custom array indices.
- In this case while loops seem to be faster than for loops.
- the base function `evalpoly` doesn’t use `@inbounds`, eventhough it improves performance with little risk of invalid indices. (The coefficient vector must guarantee that `firstindex`, `lastindex` and all integer indices in between are valid)
- chunking Horner’s method into chunks of two may improve performance further.

Here is a plot summarising performance.

 ![evalpoly-perf](https://global.discourse-cdn.com/julialang/original/3X/2/5/253a0e8f3fa87c21c52869240c4900d095ca4c15.png)

`evalpoly` refers to the version in base Julia. Here are the implementations of the others.

> **while loop**
>
> ```julia
> function ep_while(x, cs::AbstractVector)
> i = lastindex(cs)
> 
> @inbounds out = cs[i]
> i -= 1
> 
> while i >= firstindex(cs)
> @inbounds out = muladd(out, x, cs[i])
> i -= 1
> end
> 
> return out
> end
> 
> ```

> **four at a time**
>
> ```julia
> function ep_four(x, cs::AbstractVector)
> i = lastindex(cs)
> 
> @inbounds out = cs[i]
> i -= 1
> 
> fi = firstindex(cs)
> @inbounds while i > fi + 3
> out = muladd(out, x, cs[i])
> i -= 1
> out = muladd(out, x, cs[i])
> i -= 1
> out = muladd(out, x, cs[i])
> i -= 1
> out = muladd(out, x, cs[i])
> i -= 1
> end
> 
> while i >= fi
> @inbounds out = muladd(out, x, cs[i])
> i -= 1
> end
> 
> return out
> end
> 
> ```

> **Two at a time**
>
> ```julia
> function ep_two(x, cs::AbstractVector)
> i = lastindex(cs)
> 
> out = @inbounds cs[i]
> i -= 1
> 
> fi = firstindex(cs)
> @inbounds while i > fi
> out = muladd(out, x, cs[i])
> out = muladd(out, x, cs[i-1])
> 
> i -= 2
> end
> 
> return i == fi ? muladd(out, x, @inbounds(cs[fi])) : out
> end
> 
> ```

> **for loop**
>
> ```julia
> function ep_for(x, cs::AbstractVector)
> i = lastindex(cs)
> 	
> @inbounds out = cs[i]
> i -= 1
> 	
> for j in i:-1:firstindex(cs)
> out = @inbounds muladd(out, x, cs[j])
> end
> 
> return out
> end
> 
> ```

> **instruction level parallelism**
>
> ```julia
> function ep_ilp(x, cs::AbstractVector)
> fi = firstindex(cs)
> i = lastindex(cs)
> 	
> x2 = x*x
> x3 = x2*x
> 	
> @inbounds out = cs[i]
> i -= 1
> 
> @inbounds while i >= fi + 3
> # the idea here is to remove data dependencies
> # to make use of instruction level parallelism
> out = out * x3 + cs[i] * x2 + cs[i-1] * x + cs[i-2]
> i -= 3
> end
> 
> while i >= fi
> @inbounds out = muladd(out, x, cs[i])
> i -= 1
> end
> 
> return out
> end
> 
> ```

The plot uses the following benchmark

> **plot benchmark**
>
> ```julia
> setup(n) = (randn(Float64), [randn(Float64) for _ in 1:n])
> 
> bdata_base = map(ns) do k
> b = @benchmark evalpoly(x, c) setup=((x, c) = setup($k)) seconds=0.1
> minimum(b).time
> end
> 
> ```

Finally, I have a benchmark which tests performance for random numbers of coefficients. This may be closer to real world code if the number of coefficients is not constant. This benchmark yields:

```julia
evalpoly : 1.0
while loop : 0.52
four at a time : 0.65
two at a time : 0.55
for loop : 0.91
instruction level parallelism : 0.59

```

where the numbers represent the minimum time divided by the base Julia implementation’s time. Here is the benchmark code:

> **variable coefficient number benchmark**
>
> ```julia
> nsamples = 5000
> n = 20
> 
> function setup_bench(n, nsamples)
> xs = randn(Float64, nsamples)
> js = rand(1:n-1, nsamples)
> cs = [randn(Float64, j) for j in 2:n]
> 
> return xs, cs, js
> end
> 
> function bench!(f::F, (xs, cs, js)) where F
> s = 0.0
> for i in 1:length(js)
> @inbounds s += f(xs[i], cs[js[i]])
> end
> 
> return s
> end
> 
> @benchmark bench!(evalpoly, s) setup=(s = setup_bench(n, nsamples)) evals=1 samples=100_000
> 
> ```

**Questions**

- Why is a while loop so much faster?
- Might it be worth chunking into groups of two?
- Can others replicate my findings?
- Shall I make a pull request for the simple while loop implementation?

I’m on the following system

> **Version Info**
>
> ```julia
> Julia Version 1.7.2
> Commit bf53498635 (2022-02-06 15:21 UTC)
> Platform Info:
> OS: macOS (x86_64-apple-darwin19.5.0)
> CPU: Intel(R) Core(TM) i5-8257U CPU @ 1.40GHz
> WORD_SIZE: 64
> LIBM: libopenlibm
> LLVM: libLLVM-12.0.1 (ORCJIT, skylake)
> Environment:
> JULIA_NUM_THREADS = 8
> 
> ```

---

<div class="post-metadata">

**Author:** ![Luapulu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luapulu/32/17669_2.png) [@Luapulu](https://discourse.julialang.org/u/Luapulu)\
**Post date:** [April 23, 2022, 4:18pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/2 "2022-04-23T16:18:50Z")

</div>

Tagging @Mason (and @jeff.bezanson) since they wrote (and commited) the base implementation.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [April 23, 2022, 4:31pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/3 "2022-04-23T16:31:19Z")

</div>

A lot of the answer is that `evalpoly` has mostly been optimized for `Tuple` coefficients since then the whole thing gets unrolled. An `Array` is pretty much guaranteed to be at least 2x slower because it has to actually get the elements out.

---

<div class="post-metadata">

**Author:** ![Luapulu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luapulu/32/17669_2.png) [@Luapulu](https://discourse.julialang.org/u/Luapulu)\
**Post date:** [April 23, 2022, 4:35pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/4 "2022-04-23T16:35:13Z")

</div>

This is confirmed by the fact that the not `@generated` branch and the `evalpoly(x, cs::AbstractVector)` both go to the same `_evalpoly` implementation.

It’s still curious that a while loop is so much faster though.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 23, 2022, 4:47pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/5 "2022-04-23T16:47:20Z")

</div>

It’s not surprising to me that you can get a factor of 2 compared to the simple `for` loop that we have now, and it might be worth a PR.

The Base implementation was mainly optimized for the case where the coefficients are a tuple, in which case the whole computation can be unrolled, because the main usage of `evalpoly` (formerly `@evalpoly`) is currently for things like SpecialFunctions.jl where the polynomial is known statically. (An even more optimized implementation for that case is [SIMDPoly.jl](https://github.com/augustt198/SIMDPoly.jl), or [StaticPolynomials.jl](https://github.com/JuliaAlgebra/StaticPolynomials.jl) for a multivariate generalization.)

For arbitrary (dynamic length) polynomials, you’re usually better off not representing them in coefficient form (i.e. in a monomial basis), but instead using something like Chebyshev bases (e.g. as in [ApproxFun.jl](https://github.com/JuliaApproximation/ApproxFun.jl) or [FastChebInterp.jl](https://github.com/stevengj/FastChebInterp.jl)), in which case the analogue of Horner’s rule is a [Clenshaw recurrence](https://en.wikipedia.org/wiki/Clenshaw_algorithm). The famous exception, of course, is evaluating polynomials on the complex unit circle, in which case you have some type of discrete Fourier transform.

---

<div class="post-metadata">

**Author:** ![Luapulu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luapulu/32/17669_2.png) [@Luapulu](https://discourse.julialang.org/u/Luapulu)\
**Post date:** [April 23, 2022, 8:06pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/6 "2022-04-23T20:06:05Z")

</div>

I made a [PR](https://github.com/JuliaLang/julia/pull/45070#issue-1213442604).

---

<div class="post-metadata">

**Author:** ![gvdeynde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gvdeynde/32/50975_2.png) [@gvdeynde](https://discourse.julialang.org/u/gvdeynde)\
**Post date:** [May 9, 2022, 2:13pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/7 "2022-05-09T14:13:34Z")

</div>

I’m completely new to Julia (coming from Python) and I was wondering why there is no vectorized implementation for the argument `x` (i.e. evaluating the polynomial in multiple points). Google’ing led me to this (recent) question/post.

Is it intentional that a user should loop over `x` themselves (so first create a result vector `similar` to `x` and then explicitly loop? Or am I so newbie that I’m missing out on some concepts of Julia… Might be very well the case. Thanks for any clarification.

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [May 9, 2022, 2:19pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/8 "2022-05-09T14:19:12Z")

</div>

Are you thinking about broadcasting? Short version, if you have a scalar function `f(x)`, the “vectorized” version is `f.(xArray)`.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [May 9, 2022, 2:27pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/9 "2022-05-09T14:27:16Z")

</div>

> [@gvdeynde](#):
>
> why there is no vectorized implementation

you may be aware of this already, but “vectorized” in Python most of the times simply means “loop in C/C++”, since Julia is not-slow, a native Julia function can be easily applied to multiple inputs, via explicit loop, or “broadcasting” syntax `.`

---

<div class="post-metadata">

**Author:** ![gvdeynde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gvdeynde/32/50975_2.png) [@gvdeynde](https://discourse.julialang.org/u/gvdeynde)\
**Post date:** [May 9, 2022, 2:28pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/10 "2022-05-09T14:28:53Z")

</div>

Yep. But it doesn’t work for evalpoly… Or at least, I don’t know how to tell the code that the `x` argument should be broadcasted and not the `p`.

```julia
julia> p = rand(100); x = rand(1000); evalpoly.(x,p)
ERROR: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 1000 and 100")Stacktrace:
 [1] _bcs1
   @ .\broadcast.jl:516 [inlined]
 [2] _bcs
   @ .\broadcast.jl:510 [inlined]
 [3] broadcast_shape
   @ .\broadcast.jl:504 [inlined]
 [4] combine_axes
   @ .\broadcast.jl:499 [inlined]
 [5] instantiate
   @ .\broadcast.jl:281 [inlined]
 [6] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(evalpoly), Tuple{Vector{Float64}, Vector{Float64}}})
   @ Base.Broadcast .\broadcast.jl:860
 [7] top-level scope
   @ REPL[4]:1

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 9, 2022, 2:29pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/11 "2022-05-09T14:29:29Z")

</div>

> [@gvdeynde](#):
>
> Yep. But it doesn’t work for evalpoly… Or at least, I don’t know how to tell the code that the `x` argument should be broadcasted and not the `p` .

`evalpoly.(x,Ref(p))` or ` evalpoly.(x,(p,))`

See also [How to broadcast over only certain function arguments?](https://discourse.julialang.org/t/how-to-broadcast-over-only-certain-function-arguments/19274)

---

<div class="post-metadata">

**Author:** ![gvdeynde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gvdeynde/32/50975_2.png) [@gvdeynde](https://discourse.julialang.org/u/gvdeynde)\
**Post date:** [May 9, 2022, 2:32pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/12 "2022-05-09T14:32:11Z")

</div>

Thanks! I did not pick that up from reading the manual… Sorry.

Could you briefly explain what both options do exactly? Or point me to the right page in the manual?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 9, 2022, 2:33pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/13 "2022-05-09T14:33:54Z")

</div>

> [@gvdeynde](#):
>
> Or point me to the right page in the manual?

See [Broadcasting](https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting) in the manual:

> Sometimes, you want a container (like an array) that would normally participate in broadcast to be “protected” from broadcast’s behavior of iterating over all of its elements. By placing it inside another container (like a single element [`Tuple`](https://docs.julialang.org/en/v1/base/base/#Core.Tuple)) broadcast will treat it as a single value.

(I noticed that it doesn’t mention `Ref`, so I just submitted [a PR](https://github.com/JuliaLang/julia/pull/45244).)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [May 9, 2022, 3:27pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/14 "2022-05-09T15:27:04Z")

</div>

> [@stevengj](#):
>
> I noticed that it doesn’t mention `Ref`

Is there a problem with using also square brackets as a shielding container?

```julia
evalpoly.(x, [p]) 

```

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [May 9, 2022, 3:36pm UTC](https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899/15 "2022-05-09T15:36:16Z")

</div>

The compiler can currently eliminate the allocation of a `Ref`, but not an `Array`.

This should change one day, in which case creating an `Array` will be fine from a performance-perspective.

However, a `Vector` is also a 1-dimensional container, while `Ref`s are scalars.  
The length of a vector is not known at compile time, and that information also currently is also not const-propped from the constructor, either.