# Vectorization of multifold summation to speedup

I searched this website but didn't find any suitable answer describing how one can speed up summation in Mathematica using vectorization techniques and other techniques.

I often have to numerically sum over a multi-fold series of the hypergeometric type in my research work. One toy example is

lim = 150;
Sum[
Gamma[1 + n1 + n2 + n3]/(n1! n2! n3!) (0.1)^n1 (0.1)^n2 (0.1)^
n3, {n1, 0, lim}, {n2, 0, lim}, {n3, 0, lim}] // AbsoluteTiming


which takes about 42 sec on my laptop. The only way I know to speed-up is by using ParallelSum instead of Sum, which takes 9 sec, thanks to my 8 core processor. I want to know if there are any tricks or techniques to speed-up?

• The best way is probably to exploit the fact that the summands decay quite rapidly. If I set lim = 20 the summation is 100 times as fast and has still a precision of 15 digits... Dec 3, 2020 at 17:33

This may be useful as an idea, but @Henrik's observation suggests a superior approach on the case at hand. The following takes a single-core run of the OP's code from around 40s down to a little over 1s on 4 cores. Using Map instead of ParallelMap takes about three times as long. The code uses the highly efficient Outer[List,...] to construct the array of inputs and LogGamma instead of Gamma to keep the numerics in the machine-number realm.

lim = 150;
log10 = Log[0.1];
ParallelMap[
Apply@Function[{n1, n2, n3},
Total[
LogGamma[1 + n1 + n2 + n3] - LogGamma[1 + n1] -
LogGamma[1 + n2] - LogGamma[1 + n3] + (n1 + n2 + n3) log10 //
Exp, 2]],
Transpose[
Outer[List, Range[0., lim], Range[0., lim], Range[0., lim]],
{1, 3, 4, 2}]
] // Total // AbsoluteTiming

(*  {1.14043, 1.42857}  *)


Demonstration of Henrik's comment:

lim = 20;
log10 = Log[0.1];
res2 = ParallelMap[
Apply@
Function[{n1, n2, n3},
Total[
LogGamma[1 + n1 + n2 + n3] - LogGamma[1 + n1] -
LogGamma[1 + n2] - LogGamma[1 + n3] + (n1 + n2 + n3) log10 //
Exp, 2]],
Transpose[
Outer[List, Range[0., lim], Range[0., lim], Range[0., lim]],
{1, 3, 4, 2}]
] // Total // AbsoluteTiming
res - res2 // Last

(*
{0.01348, 1.42857}
-6.66134*10^-16
*)

• Thanks a lot! Indeed, Henrik's comment is useful for this toy sum, but Michael's method is useful for a general hypergeometric series. Dec 3, 2020 at 19:12

$$\Gamma (n+1)=\int_0^{\infty } t^n e^{-t} \, dt$$

So, your original expression can be rewritten as (assuming the order of integration and summation can be interchanged):

$$\int_0^{\infty } e^{-t} \left(\sum _{n=0}^L \frac{a^n t^n}{n!}\right){}^3 \, \ dt$$

The sum can be computed:

Sum[(a^n t^n)/n!, {n, 0, L}]


(E^(a t) Gamma[1 + L, a t])/Gamma[1 + L]

So, your original expression is equivalent to the following integral:

$$\int_0^{\infty } e^{-t} \left(\frac{e^{a t} \Gamma (L+1,a t)}{\Gamma (L+1)}\right)^3 \, dt$$

A function that does this computation:

s[a_, L_, opts:OptionsPattern[NIntegrate]] := NIntegrate[
Exp[-t] (Exp[a t] Gamma[L+1, a t]/Gamma[L+1])^3,
{t, 0, Infinity},
opts
]


s[1/10, 150, WorkingPrecision -> 50] //AbsoluteTiming