Exact integral of a polynomial over a simplex, with Julia

integrate_polynomial_on_simplex.jl

using TypedPolynomials
using LinearAlgebra

function integratePolynomialOnSimplex(P, S)
    gens = variables(P)
    n = length(gens)
    v = S[end]    
    B = Array{Float64}(undef, n, 0)
    for i in 1:n
        B = hcat(B, S[i] - v)
    end
    Q = P(gens => v + B * vec(gens))
    s = 0.0
    for t in terms(Q)
        coef = TypedPolynomials.coefficient(t)
        powers = TypedPolynomials.exponents(t)
        j = sum(powers)
        if j == 0
            s = s + coef
            continue
        end
        coef = coef * prod(factorial.(powers))
        s = s + coef / prod((n+1):(n+j))
    end
    return abs(LinearAlgebra.det(B)) / factorial(n) * s
end


#### --- EXAMPLE --- ####

# define the polynomial to be integrated
@polyvar x y z
P = x^4 + y + 2*x*y^2 - 3*z

#= be careful
if your polynomial does not involve one of the variables, 
e.g. P(x,y,z) = x⁴ + 2xy², it must be defined as a 
polynomial in x, y, and z; you can do:
@polyvar x y z
P = x^4 + 2*x*y^2 + 0.0*z
then check that variables(P) returns (x, y, z)
 =#

# simplex vertices
v1 = [1.0, 1.0, 1.0] 
v2 = [2.0, 2.0, 3.0] 
v3 = [3.0, 4.0, 5.0] 
v4 = [3.0, 2.0, 1.0]
# simplex
S = [v1, v2, v3, v4]

# integral
integratePolynomialOnSimplex(P, S)