How to speed up this high-dimensional definite integral?

This is what I get in my 7-years old Intel I3 in my macOS 10.13.2 and MMA ver 11.2:

Table[AbsoluteTiming[vars = Array[x, n]; 
Integrate[(Apply[Plus, vars])^-p, vars \[Element] 
ImplicitRegion[And @@ (# > 0 & /@ vars) && Apply[Plus, vars] < 1,Evaluate@vars]]], {n, 7}]

$\left( \begin{array}{cc} 0.488669 & \text{ConditionalExpression}\left[\frac{1}{1-p},\Re(p)<1\right] \\ 0.978526 & \text{ConditionalExpression}\left[\frac{\frac{1}{2-p}-1}{p-1},\Re(p)<2\right] \\ 2.07194 & \text{ConditionalExpression}\left[\frac{\frac{1}{3-p}-\frac{p}{2}}{p^2-3 p+2},\Re(p)<3\right] \\ 4.77045 & \text{ConditionalExpression}\left[\frac{p^2-2 p-\frac{6}{4-p}+3}{-6 p^3+36 p^2-66 p+36},\Re(p)<4\right] \\ 20.0045 & \text{ConditionalExpression}\left[-\frac{1}{24 (p-5)},\Re(p)<5\right] \\ 33.4052 & \text{ConditionalExpression}\left[-\frac{1}{120 (p-6)},\Re(p)<6\right] \\ 84.0427 & \text{ConditionalExpression}\left[\frac{\frac{720}{7-p}-p ((p-7) p ((p-7) p+28)+252)}{720 (p-6) (p-5) (p-4) (p-3) (p-2) (p-1)},\Re(p)<7\right] \\ \end{array} \right)$


I was able to speed things up incredibly by observing that the integrand only depends on the sum $ x_1 + x_2 + \ldots + x_n $, and we can take advantage of this by introducing a gratuitous integration over a new variable and throwing in a Dirac $ \delta $ "function", which satisfies the following relation:

$$ \int_a^b dx\,f(x)\,\delta(x-c) = \begin{cases} f(c) & a < c < b \\ 0 & c < a \text{ or } b < c \end{cases} $$

With a little coaxing, you can reproduce this in Mathematica, using DiracDelta:

Integrate[f[x]*DiracDelta[x - c], {x, a, b},
   Assumptions -> {a < b, #}] & /@ {
   a < c < b,
   c < a || b < a
  }
(* {f[c], 0} *)

Now, use a specific application of the defining relation. For $ 0 < x_1 + \ldots + x_n < 1 $, we have

$$ \frac{1}{\left(x_1 + \ldots + x_n\right)^{\,p}} = \int_0^1 dy \frac{\delta\left(y - (x_1 + \ldots + x_n\right)}{y^p} $$

We can insert this into the original expression, and swap around the order of the integration.

$$ \idotsint\limits_{V} \frac{dx_1\cdots\,dx_n}{(x_1 + \ldots + x_n)^p } = \int_0^1 \frac{dy}{y^p} \int_0^\infty dx_1 \dots \int_0^\infty dx_n \,\,\, \delta\left(y - (x_1 + \ldots + x_n\right) $$

Now we can do the integral over the $x_i$s separately

deltaFunctionIntegral[n_Integer?Positive, y_] :=
 With[{vars = Array[\[FormalX], n]},
  Apply[
   Integrate[DiracDelta[y - Total[vars]], ##, 
     Assumptions -> Positive[y]] &,
   Map[{#, 0, Infinity} &, vars]]];

Mathematica will blaze through it:

Table[deltaFunctionIntegral[n, y] // AbsoluteTiming, {n, 1, 10}]
(* 
   {{0.024917, 1}, 
    {0.033459, y}, 
    {0.042226, y^2/2}, 
    {0.045033, y^3/6}, 
    {0.050282, y^4/24}, 
    {0.05671, y^5/120}, 
    {0.063595, y^6/720}, 
    {0.073157, y^7/5040}, 
    {0.080864, y^8/40320}, 
    {0.086789, y^9/362880}}
*)

Now we can do the integral with this factor and things remain speedy:

Table[
  With[{factor = deltaFunctionIntegral[n, y]},
    Integrate[factor/y^p, {y, 0, 1}]] // AbsoluteTiming
  {n, 1, 10}]
(* 
   {{0.288039, ConditionalExpression[1/(1 - p), Re[p] < 1]}, 
    {0.300244, ConditionalExpression[1/(2 - p), Re[p] < 2]}, 
    {0.314402, ConditionalExpression[1/(2 (3 - p)), Re[p] < 3]}, 
    {0.32506, ConditionalExpression[1/(6 (4 - p)), Re[p] < 4]},
    {0.318678, ConditionalExpression[1/(24 (5 - p)), Re[p] < 5]}, 
    {0.339789, ConditionalExpression[1/(120 (6 - p)), Re[p] < 6]}, 
    {0.346237, ConditionalExpression[1/(720 (7 - p)), Re[p] < 7]}, 
    {0.35325, ConditionalExpression[1/(5040 (8 - p)), Re[p] < 8]}, 
    {0.365873, ConditionalExpression[1/(40320 (9 - p)), Re[p] < 9]}, 
    {0.367103, ConditionalExpression[1/(362880 (10 - p)), Re[p] < 10]}}
*)

I know how to show that deltaFunctionIntegral[n, y] gives $ y^{n-1}/(n-1)! $ outside of Mathematica for general $ n $, but I don't know how to coax the result out of Mathematica yet.


(This answer is just a simplification of @Pillsy's answer)

You could make the change of variables $x_i \to y x_i$ to convert the integral to $\underset{x\in \Delta^{n-1}}{\int } y^{n-1} dx$ (here $\Delta^{n-1}$ is the usual nomenclature for Simplex[n-1]) which evaluates very quickly. For example:

Integrate[y^99, x ∈ Simplex[99]] //AbsoluteTiming

{0.02938, y^99/ 933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000}

By comparison, deltaFunctionIntegral takes about 1000 times longer:

deltaFunctionIntegral[100, y] //AbsoluteTiming

{31.0896, y^99/ 933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000}