optimize @evalpoly for SIMD?

[stevengj]

Right now, @evalpoly(x, coefs...) uses Horner's rule for real x, and a more complicated algorithm for complex x, and is performance-critical for evaluation of special functions.

However, it seems likely to be faster to try and exploit SIMD instructions for these polynomial evaluations, e.g. with Estrin's algorithm or some other technique.

(The nice thing about code generation / macros is that we can investigate fancy methods for this sort of thing, protected by an isa(x,Float64) check that will get optimized out at compile-time.)

[yuyichao]

Will the new algorithm still be optimum when used in a vectorized loop?

[ViralBShah]

This may be relevant for the new libm.jl work. Cc @musm

[stevengj]

@yuyichao, probably? Even if the compiler inlines a special-function call like erfinv, I'm skeptical that compilers are powerful enough to vectorize the polynomial evaluation. e.g. I don't see any evidence of this in:

julia> function foo(a)
          s = zero(eltype(a))
          @simd for i in eachindex(a)
              s += @inbounds @evalpoly(a[i], 2.3, 4.2, 1.6, 8.9, 0.8)
          end
          return s
       end
foo (generic function with 1 method)

julia> @code_native foo(rand(10))

if I'm reading the asm output correctly.

[yuyichao]

I'm not worrying about inlining special functions that use @evalpoly, I'm worrying about the case where @evalpoly is used directly in the loop.

The code you posted has two problems.

1. It's accumulating nothing.

   Until RFC: Make inbounds macros expression-like more often #15558 is merged you need to put the @inbounds on the loop or before the assignment.

2. It does seem that the loop body is somehow too complicated for llvm to vectorize in this case (might worth investigating) However, if you are doing purely elementwise operations, a very valid use case IMHO and is even more relavant now with the broadcast syntax, i.e.

   julia> function foo2(a)
              s = zero(eltype(a))
              @simd for i in eachindex(a)
                  @inbounds a[i] = @evalpoly(a[i], 2.3, 4.2, 1.6, 8.9, 0.8)
              end
              return s
          end

   Then LLVM vectorizes it without any problem (3.8.1 at least)

P.S. I strongly recommend checking code_llvm even for people who know assembly well. It's much easier to notice the problem with @inbounds in the llvm ir.

[stevengj]

Ah, thanks. Even so, it might be nice to have a SIMD variant of @evalpoly for use in special functions.

[yuyichao]

for use in special functions.

Or cases not meant to/possible to be vectorized in a loop in general, I agree. It would be nice if the new algorithm won't regret in the vectorize loop case but I guess the answer is that it might?

[musm]

Estrin's method ... I got bored in class~~~but the following is already more efficient than @horner~~~

Only works for odd polynomials for now. I don't think simd is being attached however, could anyone confirm?

#only for odd polynomials
#c in descending order
using BenchmarkTools

x = 2.0
c = collect(Float64, 8:-1:1)
function estrin(x,c)
    n = length(c)
    m = n÷2
    ap = c
    a = similar(c)
    X = x
    while n > 1
        @simd for i = 1:m
            @inbounds a[m+1-i] = muladd(ap[n+1-(2i-1)], X, ap[n+1-(2i-1)-1])
        end
        ap[:] = a[:]
        X = X*X
        n = m
        m = m÷2
    end
    return ap[1]
end

f(x,c) = (@evalpoly x c[1] c[2] c[3] c[4] c[5] c[6] c[7] c[8])

@benchmark estrin(x,c)
@benchmark f(x,c)

There are a lot of other opportunities for parallelism for simple polynomial multiplication, though maybe simpler is better; for e.g. all the X*X for the various powers can also be computed in parallel.

[yuyichao]

What are you trying to achieve with globals?

[musm]

My bad, fixed, thanks!

[yuyichao]

There's also unnecessary copy/allocation. ap should be updated in place.

[yuyichao]

Or even better, IIUC, you can just swap the two buffers.

[yuyichao]

I don't think you can get vectorization with non-unit stride (there's trick to kind of get it working though).

I just realize this is the evalpoly issue (I thought this is one of the special/libm function ones......) in this case any allocation is already a huge loose. @evalpoly generates an expression at parse time and the replacement should do that too.

[musm]

would

    @inbounds ap[1:m] = a[1:m]

swap the buffers? or am I misunderstanding what you are saying ?