# 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. Oct 21, 2018 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 Oct 21, 2018 at 8:32
• This evaluates fine if you simply use a really big number, like 10^6 Oct 21, 2018 at 8:36
• See oeis.org/A005597 for better methods, how to compute this constant. Oct 21, 2018 at 8:56
• code-review should not be applied for "wrong" code. Oct 21, 2018 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? Oct 21, 2018 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. Oct 21, 2018 at 9:52
• Thank you for your answer and comments, especially given that it is Sunday! Oct 21, 2018 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? Oct 21, 2018 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. Oct 21, 2018 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

• Good answer, but the graph is not very informative (or, at least, not nearly as informative as it could be). Consider, for example, plotting the partial products minus the known limit value, and in log scale. You can also use the (deprecated?) SequenceLimit, or Chris Chiasson's implementation if necessary. Oct 21, 2018 at 14:13
• @Accidental, not so much "deprecated" as "buried" in a deeper context: NumericalMathNSequenceLimit[]` Oct 21, 2018 at 14:37
• mathworld.wolfram.com/TwinPrimesConstant.html Oct 25, 2018 at 14:28