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