# Double series over primes

I'm very curious if the following double series over primes has a closed form: $$\sum_{k \in \mathcal{P}}\sum_{n \in \mathcal{P}}\frac{1}{k\;n(k+n)^2}$$

where $\mathcal{P}$ denotes the set of all prime numbers.
Could you help here with Mathematica coding?

• Why did you vote for close? I'm curious about that. – user 1357113 Sep 17 '13 at 18:20
• Chris' sis, whoever it was would like you to provide more details and ideally some Mathematica code. – Verbeia Sep 17 '13 at 20:16
• I voted to close because this looks like a math problem, not a Mathematica one. You're summing over primes, which makes this a special problem which CAS can't be expected to solve symbolically. The command to request the solution is clear (you are familiar with Sum), but you haven't indicated what you are expecting beyond this. (Asking for numerical approximations would be reasonable.) – Szabolcs Sep 17 '13 at 20:40
• Yes, I should have clearly left a comment. To sum up the reasons: 1. Asking to do it with Mathematica doesn't automatically make it a Mathematica problem. The bulk of the work is math here. 2. No effort shown to solve the problem. – Szabolcs Sep 17 '13 at 20:48
• All, given the divided views on whether the question should be closed, I propose to give it the benefit of the doubt and leave it open. If a question has this many upvoted answers, it deserves to be left open. – Verbeia Sep 17 '13 at 23:31

In general Mathematica cannot compute symbolically infinite sums over primes because of the lack of appropriate mathematical tools. However there are infinite products over primes which are basically well understood on the mathematical level. One famous example is the Euler formula for the Riemann zeta function, one of the most beautiful (and mysterious even though its proof can be easily understood) mathematical treasures: for Re[z] > 1

Sum[ 1/n^z, {n, Infinity}] == Defer[ Sum[ 1/n^z, {n, Infinity}]] ==
Defer[Product[ 1/(1 - Prime[i]^-z), {i, Infinity}]] // TraditionalForm


Taking this formula into account we can easily write related symbolic sum over primes e.g.

Sum[ -Log[1 - Prime[i]^-z], {i, Infinity}]

Log[Zeta[z]]


This works also for numerical values, e.g.

Defer[ Sum[ -Log[ 1 - Prime[i]^-2], {i, Infinity}]] ==
Sum[ -Log[ 1 - Prime[i]^-2], {i, Infinity}] // TraditionalForm


I strongly recommend reading Primes - a free book by Barry Mazur and William Stein, discussing relations between primes and the Riemann zeta function. The most important mathematical problem (according to Hilbert the most important at all) is the Riemann Hypothesis, still open, even though generally believed to be true.

By no means Mathematica would have computed infinite sums if there had been no appropriate mechanism recognizing well known patterns leading to more general formulae for symbolic sums as the above definition of the Riemann zeta function.
I doubt Mathematica knows another infinite symbolic sums over primes.

Prime[n] can be computed up to only certain big number, there is an arbitrary cut-off, namely OmegaPrime = 7783516045221;. For details on the issue see What is so special about Prime?.

On the other hand it appears to know asymptotic density of primes when testing certain sum convergence:

SumConvergence[ 1/Prime[n], n]

False


There are related mathematical functions closely related to the issue of distribution of primes: RiemannR, PrimePi, LogIntegral etc.

We can proceed finding only numerical approximations of sums over primes.

One way would be

I

defining:

g[n_Integer, k_Integer] := Boole[PrimeQ[n] && PrimeQ[k]]


We can compute numerical approximation with NSum, e.g. NSum[ g[n, k]/(k n (k + n)^2), {k, ∞}, {n, ∞}] // Quiet, we can get the result more accurate by increasing NSumTerms (for a related discussion see Precision differences):

NSum[ g[n, k]/(k n (k + n)^2), {k, ∞}, {n, ∞}, NSumTerms -> 200,
AccuracyGoal -> 10, PrecisionGoal -> 10] // Quiet

0.0448588


This took roughly 10 minutes on my computer.

There is however another, slightly faster approach:

II

Since the series is monotonical and bounded we could compute exactly finite number of terms (just like NSum does behind the scenes) having quite a good idea of the error. This method took about 3 minutes:

Total[
Array[ 1/(Prime[#1] Prime[#2] ( Prime[#1] + Prime[#2])^2)&, {3000, 3000}], 2]//
N[ #, 10]&

0.04486521704

• Perhaps the prime number theorem could be used to improve the numerical approximation. I'm not sure if this could be feasible or not. If we have more digits, there's a better chance to make use of inverters. – Szabolcs Sep 18 '13 at 1:20
• @Szabolcs It seems that Mathematica sometimes is (partially) aware of the prime number theorem, e.g. evaluation of Sum[1/Prime[n], {n, Infinity}] says that Sum::div: Sum does not converge., on the other hand there may appear another issues where it is not aware even more obvious facts. – Artes Sep 18 '13 at 10:40

## !! WARNING !!

Given the above circumstances the answer below is completely wrong and should not be replicated rather than used for guideline. The problem lies in the fact that Assumptions DON'T affect summation indices hence they are disregarded in the code below.

You can evaluate the series using the following snippet:

Sum[1/(k n (k + n)^2), {k, 1, ∞}, {n, 1, ∞},
Assumptions -> {k ∈ Primes, n ∈ Primes}]


Which evaluates to:

π^4/180

• I upvoted this answer a few hours ago without thinking much about it. Now, it is clear that this is utterly wrong. HOWEVER it is "food for thought", so please instead of deleting it, edit it writing a few lines. Start with a warning for future users and then explain why it is wrong. You can keep my upvote. I learnt something :) – Dr. belisarius Sep 18 '13 at 2:55

For numerical approximation try this

primes = N@Array[Prime, 1000];
res = Sum[1/(k n (k + n)^2), {k, primes}, {n, primes}]


0.0448652

The most straightforward way of evaluating your sum in Mathematica is as follows:

Sum[1/(Prime[k] Prime[n] (Prime[k] + Prime[n])^2), {k, 1, Infinity}, {n, 1, Infinity}]


See for instance the last example on the how-to EvaluateInfiniteSumsAndProducts. This returns unevaluated however, so it seems that if there is a closed form for this series then Mathematica doesn't know it.

Here is a way to get close.

ParallelTable[ Sum[1/(k n (k + n)^2), {k, N@Prime[Range[i]]},
{n, N@Prime[Range[i]]}], {i, 10, 2000, 10}]


This gives the following output:

{0.0445365, 0.0448078, 0.0448455, 0.0448556, 0.0448597, 0.0448617, <<188>>, 0.0448652, 0.0448652, 0.0448652, 0.0448652, 0.0448652, 0.0448652}

Which can be plotted to get a sense of the trend within:

One can then use Rationalize

Rationalize[0.0448652, 0.00000001]


Which gives you:

218/4859

• Heh, I'm not quite certain it makes sense to rationalize here, since the asymptotic value is unknown... the numerator and denominator could very well be anything :) For instance, it seems pretty close to $$\frac{5e}{\pi^5}\left(1+10^{-2}+10^{-4}\right)$$ :D – rm -rf Sep 17 '13 at 23:52
• @rm-rf. You're right. But I only claimed what I showed here is "a way to get close" – RunnyKine Sep 18 '13 at 0:10
• Of course, your answer is helpful. Don't mind me here... I just needed a place to plug in my stupid "approximation" :P – rm -rf Sep 18 '13 at 0:12
• @rm-rf. Actually, it's a good thing you did. That way some random reader don't take my post as "valid". Thanks for sharing the "approximation" – RunnyKine Sep 18 '13 at 0:15