7

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?

Klangen
  • 1,009
  • 5
  • 14

2 Answers2

10

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

SonerAlbayrak
  • 2,088
  • 10
  • 11
  • 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 09: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. – SonerAlbayrak Oct 21 '18 at 09: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. – SonerAlbayrak Oct 21 '18 at 18:15
7

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

enter image description here

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309