# Computing the twin prime constant with Mathematica

The twin prime constant is defined as

$$C _2 = { \prod \limits_{p \geqslant 3 }} \, \frac{p(p-2)}{(p-1)^2}.$$

Translating this directly in Mathematica code gives:

N[Product[Prime[i] (Prime[i] - 2)/(Prime[i] - 1)^2, {i, 2, Infinity}]]


However, when I execute this code I get the following error messages:

What am I doing wrong?

• It's just that some of the algorithms used internally by NProduct[] sometimes try to evaluate terms for inexact arguments, so a more clever approach is necessary. – J. M. will be back soon Oct 21 '18 at 8:09
• I didn't get this message. However I suggest not using Infinity inside Product, as there will possibly be no results. Try a big number instead, like 100000 – t-smart Oct 21 '18 at 8:32
• This evaluates fine if you simply use a really big number, like 10^6 – Carl Lange Oct 21 '18 at 8:36
• See oeis.org/A005597 for better methods, how to compute this constant. – Vaclav Kotesovec Oct 21 '18 at 8:56
• code-review should not be applied for "wrong" code. – user202729 Oct 21 '18 at 14:39

The problem with your code is that N[Product[...]] calls NProduct[...] which turns integers into reals before Prime can evaluate, causing the error: See Details and Options section of documentation of NProduct.

For your calculation, you actually do not need these. It looks like the result converges around $$0.660162$$:

Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, {i, 2, 10^k}]], {k, 5}]


{0.665138, 0.66033, 0.66017, 0.660162, 0.660162}

Since the product is monotonic and you are using N at the end, it is not necessary to carry out whole multiplication all the way upto infinity: Depending on the interested accuracy, you can go to higher orders as well. For example,

Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, {i, 2, 10^k}], 8], {k, 5}]


yields

{0.66513840, 0.66033029, 0.66017020, 0.66016232, 0.66016185}

which reveals that using only first $$10^5$$ primes, which is all primes upto $$1299709$$, is sufficient to get a result with error of the order $$10^{-6}$$.

For the record: Mathematica fails to do the calculation analytically:

Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, {i, 2, Infinity}]


$$\prod _{i=2}^{\infty } \frac{\left(p_i-2\right) p_i}{\left(p_i-1\right){}^2}$$

• Thank you. So if I want an error of maximum $10^{-x}$, I only need to use the first $10^{x-1}$ primes? – Klangen Oct 21 '18 at 9:47
• @PierreTheFermented No, not necessarily. I was just trying to say that you do not need all primes and depending on the accuracy needed, different number of primes will be sufficient. I am not sure if there is an analytical way to predict that number though. In any case, using more than first 10^6 primes will probably be fairly slow. – Soner Oct 21 '18 at 9:52
• Thank you for your answer and comments, especially given that it is Sunday! – Klangen Oct 21 '18 at 10:19
• As a final remark, if N[Product[...]] calls NProduct[...], can't we use Rationalize[] for the i within the product so that they don't get evaluated as reals? – Klangen Oct 21 '18 at 10:25
• I tried it out of curiosity but it gave the same errors, so apparently it does not help. Honestly, if you try an infinite sum, it makes much more sense to avoid mathematically complicated commands such as Prime or Rationalize to begin with. – Soner Oct 21 '18 at 18:15

As not all primes are known, Mathematica will have a problem to compute the constant exactly. But the sequence appears to converge quickly to some number $$0.660162...$$. (Soner posted that already.)

So let me give a way to compute the sequence of partial products more efficientlt.

n = 10000;
seq = FoldList[#1 #2 (#2 - 2)/(#2 - 1)^2 &, 1, Prime[Range[2, n]]];
N[seq[[-1]]]
ListLinePlot[{N[seq], ConstantArray[N[seq[[-1]]], Length[seq]]},
PlotRange -> {1/2, 1},
PlotStyle -> {Thick, Dashed}
]


0.660162