Hello,

I wrote a function which computes the exact value of a multivariate polynomial over a simplex. I’m open to any comment on my code.

I also implemented this method in Python. I’m sure it is correct, I multi-checked.

Integration on simplices is important because any convex polyhedron can be partitionned into simplices, thanks to the Delaunay tessellation. Hence integration on simplices allows integration on any convex polyhedron.

By the way, does there exist a Julia package for integration of an arbitrary function over a simplex?

```
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)
```