# Evaluation of a triple sum does not finish in reasonable time

I'm trying to compute the following triple sum, but no result is produced within a reasonable amount of time. What to do?

Sum[1/( i j k (i + j + k + 1)), {i, 1, Infinity}, {j, 1, Infinity}, {k, 1, Infinity}]


I use Mathematica 8.0.

• This works... ParallelSum[Sum[Sum[1/(i j k (i + j + k + 1)), {i, 1, Infinity}] /. j -> j1, {k, 1, Infinity}], {j1, 1, 50}]  but note that it stops at j=50 May 8, 2013 at 9:39
• FYI its value is 5.35132... does not seem to converge quickly though... May 8, 2013 at 10:30
• note that the value I gave you is for jmax=50 not infinity! May 8, 2013 at 10:38
• @chris You can get a closed form for one summation, it might help convergence. May 8, 2013 at 12:39
• @chris and are the two of you siblings ? May 8, 2013 at 12:40

@whuber gave me another idea : first change the summand to an exponential and then do the sums using the symmetry.

Integrate[Exp[-t i j k (i + j + k + 1)], {t, 0, Infinity}, Assumptions -> {i > 0, j > 0, k > 0}]
(* 1/(i j k (1 + i + j + k)) *)


so we need to do the sums of the exponential and then integrate at the end. Now we can make a change of variables t i j k = y which will simplify the exponent and bring a factor 1/(i j k) from the Jacobian. This makes it easier to see that

$$\sum_{i,j,k=1}^{\infty} \frac{1}{i j k} \exp(-y\ (i+j+k+1)) = \exp(-y) \left( \sum_{i=1}^{\infty} \frac{1}{i}\exp(-y\ i)\right)^3$$

Simplify[Sum[Exp[-y i]/i, {i, 1, Infinity}, Assumptions -> {y > 0}],
Assumptions -> {y > 0}]
(* -Log[1 - E^-y] *)


and we are left with the integration :

Integrate[Exp[-y] (-Log[1 - E^-y])^3, {y, 0, Infinity}]
(* 6 *)


It's also straightforward to generalize to the sum over n integers and the result is :

Integrate[Exp[-y] (-Log[1 - E^-y])^n, {y, 0, Infinity},
Assumptions -> {n \[Element] Integers, n >= 1}]
(* Gamma[1 + n] *)

• Aside from justifying switching the order of integration and summation and a typo (missing $i$ in second summation), it looks good. Numerically I got an upper bound of about 6.08 on the sum. (+1) May 8, 2013 at 19:26
• @MichaelE2 Fixed, thanks a lot. May 8, 2013 at 19:28
• Very nice idea! @Michael E2 has clearly identified my mistake, which now squares our two results. (I promised only to find a lower bound for the sum, so a numerical value of the integral around $3.5$ is just fine). I apologize for making you work so hard :-) and will make up for it by transferring all the points from my answer to yours, which deserves it. (There's a day's wait before I can do that, though...) May 8, 2013 at 19:37
• @whuber You had the key insight and deserve most of the credit. May 8, 2013 at 19:44
• btw, congrats on your 10k!
– rm -rf
May 11, 2013 at 22:20

Here is another variant along the same lines as @whuber proposed. Consider function $$f(x,y,z,t)=\sum _{i=1}^{\infty } \sum _{j=1}^{\infty } \sum _{k=1}^{\infty } x^{i-1} y^{j-1} z^{k-1} t^{i+j+k}.$$ Mma gives

f[x_, y_, z_, t_] =
Sum[x^(i - 1) y^(j - 1) z^(k - 1) t^(i + j + k), {i,
1, \[Infinity]}, {j, 1, \[Infinity]}, {k, 1, \[Infinity]}]


$$-\frac{t^3}{(t x-1) (t y-1) (t z-1)}.$$ Then the required sum is equal to $$\int _0^1\int _0^1\int _0^1\int _0^1f(x,y,z,t)\,dzdydxdt,$$

Integrate[f[x, y, z, t], {t, 0, 1}, {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]


which gives 6 as in the answer of @b.gatessucks.

If I take the challenge to be to get Mathematica to do the sum in the form given (as distinct from applying some mathematics insightfully), then here is a way to get the answer. The only trick I use is to break the triple sum into a double and single sum.

s2 = Sum[1/(i j (i + j + k + 1)), {i, 1, Infinity}, {j, 1, Infinity}]; // Timing
(s2 = FullSimplify[s2]) // Timing
Sum[s2/k, {k, 1, Infinity}] // Timing

{93.696694, Null}
{0.234598, (\[Pi]^2 + 6 HarmonicNumber[1 + k]^2 - 6 PolyGamma[1, 2 + k])/(6 + 6 k)}
{12.027743, 6}


The intermediate step of simplifying s2 seems to be necessary.

Interestingly Mathematica caches auxiliary expressions in computing the sum. The second time around is faster. (Almost all the time is spent in FullSimplify)

(s2 = Sum[1/(i j (i + j + k + 1)), {i, 1, Infinity}, {j, 1, Infinity}];
s2 = FullSimplify[s2];
Sum[s2/k, {k, 1, Infinity}]) // Timing

{0.228539, 6}


Not an answer, but it might help. You can get one summation in closed form :

sum1 = Sum[1/(i j k (i + j + k + 1)), {i, 1, Infinity}]
(* (EulerGamma + PolyGamma[0, 2 + j + k])/(j k (1 + j + k)) *)


You can also do the easy bit :

sum2 = Sum[Sum[EulerGamma/(j k (1 + j + k)), {j, 1, Infinity}], {k, 1, Infinity}]
(* 2 EulerGamma *)


so you're left with a double sum

sum = 2 EulerGamma +
Sum[PolyGamma[0, 2 + j + k]/(j k (1 + j + k)), {j, 1, Infinity}, {k, 1, Infinity}]


You can make some further simplification by using the symmetry of the summand or proceed numerically.