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

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)

R implementation:

library(spray)

integratePolynomialonSimplex <- function(P, S) {
  n <- ncol(S)
  v <- S[n+1L, ]
  B <- t(S[1L:n, ]) - v
  gens <- lapply(1L:n, function(i) lone(i, n))
  newvars <- vector("list", n)
  for(i in 1L:n) {
    newvar <- v[i]
    Bi <- B[i, ]
    for(j in 1L:n) {
      newvar <- newvar + Bi[j] * gens[[j]]
    }
    newvars[[i]] <- newvar
  }
  Q <- 0
  exponents <- P[["index"]]
  coeffs    <- P[["value"]] 
  for(i in 1L:nrow(exponents)) {
    coef <- coeffs[i]
    powers <- exponents[i, ]
    term <- 1
    for(j in 1L:n) {
      term <- term * newvars[[j]]^powers[j] 
    }
    Q <- Q + coef * term
  }
  s <- 0
  exponents <- Q[["index"]]
  coeffs    <- Q[["value"]] 
  for(i in 1L:nrow(exponents)) {
    coef <- coeffs[i]
    powers <- exponents[i, ]
    j <- sum(powers)
    if(j == 0L) {
      s <- s + coef
      next
    }
    coef <- coef * prod(factorial(powers))
    s <- s + coef / prod((n+1L):(n+j))
  }
  abs(det(B)) / factorial(n) * s
}

############################# --- EXAMPLE --- ##################################
library(spray)

# variables
x <- lone(1, 3)
y <- lone(2, 3)
z <- lone(3, 3)
# polynomial
P <- x^4 + y + 2*x*y^2 - 3*z

# simplex (tetrahedron) vertices
v1 <- c(1, 1, 1)
v2 <- c(2, 2, 3)
v3 <- c(3, 4, 5)
v4 <- c(3, 2, 1)
# simplex
S <- rbind(v1, v2, v3, v4)

# integral
integratePolynomialonSimplex(P, S)
```

from math import factorial
from sympy import Poly
import numpy as np

def term(Q, monom):
    coef = Q.coeff_monomial(monom)
    powers = list(monom)
    j = sum(powers)
    if j == 0:
        return coef
    coef = coef * np.prod(list(map(factorial, powers)))
    n = len(monom)
    return coef / np.prod(list(range(n+1, n+j+1)))

def integral(P, S):
    gens = P.gens
    n = len(gens)
    dico = {}
    v = S[n,:]
    columns = []
    for i in range(n):
        columns.append(S[i,:] - v)    
    B = np.column_stack(tuple(columns))
    for i in range(n):
        newvar = v4[i]
        for j in range(n):
            newvar = newvar + B[i,j]*Poly(gens[j], gens, domain="RR")
        dico[gens[i]] = newvar.as_expr()
    Q = P.subs(dico, simultaneous=True).as_expr().as_poly(gens)
    monoms = Q.monoms()
    s = 0.0
    for monom in monoms:
        s = s + term(Q, monom)
    return np.abs(np.linalg.det(B)) / factorial(n) * s

###############################################################################
# tetrahedron vertices
v1 = np.array([1.0, 1.0, 1.0]) 
v2 = np.array([2.0, 2.0, 3.0]) 
v3 = np.array([3.0, 4.0, 5.0]) 
v4 = np.array([3.0, 2.0, 1.0])
S = np.array([v1, v2, v3, v4])

# polynomial to integrate
from sympy.abc import x, y, z
P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR")

# integral
val = integral(P, S)
print(val)
```

from math import factorial
from sympy import Poly
import numpy as np

def term(Q, monom):
    coef = Q.coeff_monomial(monom)
    powers = list(monom)
    j = sum(powers)
    if j == 0:
        return coef
    coef = coef * np.prod(list(map(factorial, powers)))
    n = len(monom)
    return coef / np.prod(list(range(n+1, n+j+1)))

def integral(P, S):
    gens = P.gens
    n = len(gens)
    dico = {}
    v = S[n,:]
    columns = []
    for i in range(n):
        columns.append(S[i,:] - v)    
    B = np.column_stack(tuple(columns))
    for i in range(n):
        newvar = v4[i]
        for j in range(n):
            newvar = newvar + B[i,j]*Poly(gens[j], gens, domain="RR")
        dico[gens[i]] = newvar.as_expr()
    Q = P.subs(dico, simultaneous=True).as_expr().as_poly(gens)
    monoms = Q.monoms()
    s = 0.0
    for monom in monoms:
        s = s + term(Q, monom)
    return np.abs(np.linalg.det(B)) / factorial(n) * s

###############################################################################
# tetrahedron vertices
v1 = np.array([1.0, 1.0, 1.0]) 
v2 = np.array([2.0, 2.0, 3.0]) 
v3 = np.array([3.0, 4.0, 5.0]) 
v4 = np.array([3.0, 2.0, 1.0])
S = np.array([v1, v2, v3, v4])

# polynomial to integrate
from sympy.abc import x, y, z
P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR")

# integral
val = integral(P, S)
print(val)

Julia implementation:

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

from math import factorial
from sympy import Poly
from sympy.abc import x, y, z
import numpy as np

def term(Q, monom):
    coef = Q.coeff_monomial(monom)
    powers = list(monom)
    j = sum(powers)
    if j == 0:
        return coef
    coef = coef * np.prod(list(map(factorial, powers)))
    return coef / np.prod(list(range(4, 4+j)))

def integral(P, v1, v2, v3, v4):
    X = Poly(x, x, y, z, domain="RR")
    Y = Poly(y, x, y, z, domain="RR")
    Z = Poly(z, x, y, z, domain="RR")
    B = np.column_stack((v1-v4, v2-v4, v3-v4))
    newx = B[0,0]*X + B[0,1]*Y + B[0,2]*Z + v4[0]
    newy = B[1,0]*X + B[1,1]*Y + B[1,2]*Z + v4[1]
    newz = B[2,0]*X + B[2,1]*Y + B[2,2]*Z + v4[2]
    let = {x: newx.as_expr(), y: newy.as_expr(), z: newz.as_expr()}
    Q = P.subs(let, simultaneous=True).as_expr().as_poly(x, y, z)
    monoms = Q.monoms()
    s = 0.0
    for monom in monoms:
        s = s + term(Q, monom)
    return np.abs(np.linalg.det(B)) / 6 * s

###############################################################################
# tetrahedron vertices
v1 = np.array([1.0, 1.0, 1.0]) 
v2 = np.array([2.0, 2.0, 3.0]) 
v3 = np.array([3.0, 4.0, 5.0]) 
v4 = np.array([3.0, 2.0, 1.0])

# polynomial to integrate    
P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR")

# integral
val = integral(P, v1, v2, v3, v4)
print(val)
``` 

from math import factorial
from sympy import Poly
from sympy.abc import x, y, z
import numpy as np

def term(Q, monom):
    coef = Q.coeff_monomial(monom)
    powers = list(monom)
    j = sum(powers)
    if j == 0:
        return coef
    coef = coef * np.prod(list(map(factorial, powers)))
    return coef / np.prod(list(range(4, 4+j)))

def integral(P, v1, v2, v3, v4):
    X = Poly(x, x, y, z, domain="RR")
    Y = Poly(y, x, y, z, domain="RR")
    Z = Poly(z, x, y, z, domain="RR")
    B = np.column_stack((v1-v4, v2-v4, v3-v4))
    newx = B[0,0]*X + B[0,1]*Y + B[0,2]*Z + v4[0]
    newy = B[1,0]*X + B[1,1]*Y + B[1,2]*Z + v4[1]
    newz = B[2,0]*X + B[2,1]*Y + B[2,2]*Z + v4[2]
    let = {x: newx.as_expr(), y: newy.as_expr(), z: newz.as_expr()}
    Q = P.subs(let, simultaneous=True).as_expr().as_poly(x, y, z)
    monoms = Q.monoms()
    s = 0.0
    for monom in monoms:
        s = s + term(Q, monom)
    return np.abs(np.linalg.det(B)) / 6 * s

###############################################################################
# tetrahedron vertices
v1 = np.array([1.0, 1.0, 1.0]) 
v2 = np.array([2.0, 2.0, 3.0]) 
v3 = np.array([3.0, 4.0, 5.0]) 
v4 = np.array([3.0, 2.0, 1.0])

# polynomial to integrate    
P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR")

# integral
val = integral(P, v1, v2, v3, v4)
print(val)

EDIT

For an arbitrary $n$:

from math import factorial
from sympy import Poly
import numpy as np

def term(Q, monom):
    coef = Q.coeff_monomial(monom)
    powers = list(monom)
    j = sum(powers)
    if j == 0:
        return coef
    coef = coef * np.prod(list(map(factorial, powers)))
    n = len(monom)
    return coef / np.prod(list(range(n+1, n+j+1)))

def integral(P, S):
    gens = P.gens
    n = len(gens)
    dico = {}
    v = S[n,:]
    columns = []
    for i in range(n):
        columns.append(S[i,:] - v)    
    B = np.column_stack(tuple(columns))
    for i in range(n):
        newvar = v4[i]
        for j in range(n):
            newvar = newvar + B[i,j]*Poly(gens[j], gens, domain="RR")
        dico[gens[i]] = newvar.as_expr()
    Q = P.subs(dico, simultaneous=True).as_expr().as_poly(gens)
    monoms = Q.monoms()
    s = 0.0
    for monom in monoms:
        s = s + term(Q, monom)
    return np.abs(np.linalg.det(B)) / factorial(n) * s

###############################################################################
# tetrahedron vertices
v1 = np.array([1.0, 1.0, 1.0]) 
v2 = np.array([2.0, 2.0, 3.0]) 
v3 = np.array([3.0, 4.0, 5.0]) 
v4 = np.array([3.0, 2.0, 1.0])
S = np.array([v1, v2, v3, v4])

# polynomial to integrate
from sympy.abc import x, y, z
P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR")

# integral
val = integral(P, S)
print(val)
```