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?

  • $\begingroup$ Why did you vote for close? I'm curious about that. $\endgroup$ Sep 17, 2013 at 18:20
  • 1
    $\begingroup$ Chris' sis, whoever it was would like you to provide more details and ideally some Mathematica code. $\endgroup$
    – Verbeia
    Sep 17, 2013 at 20:16
  • 2
    $\begingroup$ 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.) $\endgroup$
    – Szabolcs
    Sep 17, 2013 at 20:40
  • 1
    $\begingroup$ 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. $\endgroup$
    – Szabolcs
    Sep 17, 2013 at 20:48
  • 2
    $\begingroup$ 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. $\endgroup$
    – Verbeia
    Sep 17, 2013 at 23:31

5 Answers 5


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

enter image description here

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

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

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

enter image description here

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]

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



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

This took roughly 10 minutes on my computer.

There is however another, slightly faster approach:


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:

  Array[ 1/(Prime[#1] Prime[#2] ( Prime[#1] + Prime[#2])^2)&, {3000, 3000}], 2]//
N[ #, 10]&
  • 1
    $\begingroup$ 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. $\endgroup$
    – Szabolcs
    Sep 18, 2013 at 1:20
  • $\begingroup$ @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. $\endgroup$
    – Artes
    Sep 18, 2013 at 10:40


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:

  • 8
    $\begingroup$ 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 :) $\endgroup$ Sep 18, 2013 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}]



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:

enter image description here

One can then use Rationalize

Rationalize[0.0448652, 0.00000001]

Which gives you:


  • 4
    $\begingroup$ 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 $\endgroup$
    – rm -rf
    Sep 17, 2013 at 23:52
  • $\begingroup$ @rm-rf. You're right. But I only claimed what I showed here is "a way to get close" $\endgroup$
    – RunnyKine
    Sep 18, 2013 at 0:10
  • 1
    $\begingroup$ Of course, your answer is helpful. Don't mind me here... I just needed a place to plug in my stupid "approximation" :P $\endgroup$
    – rm -rf
    Sep 18, 2013 at 0:12
  • $\begingroup$ @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" $\endgroup$
    – RunnyKine
    Sep 18, 2013 at 0:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.