New evaluation routine

benchmark/evaluation.jl

@@ -0,0 +1,23 @@
+using BenchmarkTools
+import DynamicPolynomials: @polyvar
+using FixedPolynomials
+
+@polyvar x y
+f1 = x + y + 1
+f2 = (x + 1)^4 + y
+f3 = 50x^3+83x^2*y+24x*y^2+y^3+392x^2+414x*y+50y^2-28x+59y-100
+f4 = x+x^2+y+x*y^3+x^4*y^2+3x^3*y+10*x*y+10x^2*y+10x*y^2+15x^2*y^2+10x^3*y^2
+
+function make(f)
+    p = Polynomial{Float64}(f)
+    w = rand(2)
+    result = @benchmark evaluate($p, $w)
+    println(STDOUT, f)
+    show(STDOUT, MIME"text/plain"(), result)
+    println("\n")
+end
+
+make(f1)
+make(f2)
+make(f3)
+make(f4)

src/poly.jl

@@ -33,9 +33,32 @@ struct Polynomial{T<:Number}
     variables::Vector{Symbol}
     homogenized::Bool
 
+    # this are only internal data structures to avoid allocations
+    # Let M the exponents matrix. Then we can associate a matrix M' in which each
+    # row of m was sorted separetly in ascending order.
+    # `_phi` is the  map from M' to M s.t. M[i,j] == M'[i, _phi[i,j]]
+    # For each row of M' we can only store the difference to the previous value.
+    # This is `_sorteddiff`.
+    # To evaluate a polynomial we then use '_phi' and '_sorteddiff' to compute the matrix
+    # `x.^M` and store the result in `_values`.
+    _phi::Matrix{Int}
+    _sorteddiff::Matrix{Int}
+    _values::Matrix{T}
+
     function Polynomial{T}(exponents::Matrix{Int}, coefficients::Vector{T}, variables::Vector{Symbol}, homogenized::Bool) where {T<:Number}
         sorted_cols = sort!([1:size(exponents,2);], lt=((i, j) -> lt_total_degree(exponents[:,i], exponents[:,j])), rev=true)
-        new(exponents[:, sorted_cols], coefficients[sorted_cols], variables, homogenized)
+
+        exps = exponents[:, sorted_cols]
+        coeffs = coefficients[sorted_cols]
+
+        psi = vcat([sortperm(exps[k,:])' for k=1:size(exps, 1)]...)
+        phi = [findfirst(x -> x == j, psi[i,:]) for i=1:size(exps, 1), j=1:size(exps, 2)]
+        sorted = vcat([sort(exps[i,:])' for i in 1:size(exps,1)]...)
+        sorteddiff = sorted
+        sorteddiff[:,2:end] -= sorteddiff[:,1:end-1]
+        values = zeros(T, size(exps))
+
+        new(exps, coeffs, variables, homogenized, phi, sorteddiff, values)
     end
 end
 
@@ -168,15 +191,27 @@ Evaluates `p` at `x`, i.e. ``p(x)``.
 `Polynomial` is also callable, therefore you can also evaluate it via `p(x)`.
 """
 function evaluate(p::Polynomial{S}, x::AbstractVector{T}) where {S<:Number, T<:Number}
-    cfs = coefficients(p)
-    exps = exponents(p)
-    nvars, nterms = size(exps)
-    res = zero(promote_type(S,T))
-    for j = 1:nterms
-        @inbounds term = p.coefficients[j]
-        for i = 1:nvars
-            k = exps[i, j]
-            @inbounds term *= x[i]^k
+    m, n = size(p.exponents)
+    values = p._values
+    sorteddiff = p._sorteddiff
+    phi = p._phi
+
+    for i = 1:m
+        @inbounds values[i,1] = x[i]^sorteddiff[i,1]
+    end
+
+    for i=1:m, k=2:n
+        @inbounds l = sorteddiff[i,k]
+        @inbounds v = values[i,k - 1]
+        @inbounds values[i,k] = l == 0 ? v : (l == 1 ? v * x[i] : v * x[i]^sorteddiff[i,k])
+    end
+
+    res = 0.0
+    c = p.coefficients
+    for k = 1:n
+        @inbounds term = c[k]
+        for i = 1:m
+            @inbounds term *= p._values[i, phi[i,k]]
         end
         res += term
     end