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:

enter image description here

What am I doing wrong?

  • 3
    $\begingroup$ 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. $\endgroup$ Oct 21, 2018 at 8:09
  • 1
    $\begingroup$ 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 $\endgroup$
    – t-smart
    Oct 21, 2018 at 8:32
  • 1
    $\begingroup$ This evaluates fine if you simply use a really big number, like 10^6 $\endgroup$
    – Carl Lange
    Oct 21, 2018 at 8:36
  • 4
    $\begingroup$ See oeis.org/A005597 for better methods, how to compute this constant. $\endgroup$ Oct 21, 2018 at 8:56
  • $\begingroup$ code-review should not be applied for "wrong" code. $\endgroup$
    – user202729
    Oct 21, 2018 at 14:39

2 Answers 2


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}]


{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}$

  • $\begingroup$ Thank you. So if I want an error of maximum $10^{-x}$, I only need to use the first $10^{x-1}$ primes? $\endgroup$
    – Klangen
    Oct 21, 2018 at 9:47
  • $\begingroup$ @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. $\endgroup$ Oct 21, 2018 at 9:52
  • $\begingroup$ Thank you for your answer and comments, especially given that it is Sunday! $\endgroup$
    – Klangen
    Oct 21, 2018 at 10:19
  • $\begingroup$ 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? $\endgroup$
    – Klangen
    Oct 21, 2018 at 10:25
  • $\begingroup$ 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. $\endgroup$ 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]]];
ListLinePlot[{N[seq], ConstantArray[N[seq[[-1]]], Length[seq]]},
 PlotRange -> {1/2, 1},
 PlotStyle -> {Thick, Dashed}


enter image description here


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.