One can compute the amount of twin primes below a positive integer $n$ by using the Mathematica command (taken from OEIS A001097):

Length[Select[Prime[Range[n]], PrimeQ[# + 2] &]]

The twin prime conjecture states that this value should approach $1.320323632\ldots\times\int_2^n \frac{dt}{\log^2 t}$. So I tried using


However, I got vastly different results than I expected. For instance, for $n=1000$ and $n=10000$ I get, respectively, $45.8\ldots$ and $214.21\ldots$, while the real values are $174$ and $1270$. Obviously, there is something wrong with the Mathematica command above. But what is it?

  • 2
    $\begingroup$ You've posted this in a blog already? mathematics.filegala.com/2018/09/14/… $\endgroup$ – JimB Sep 14 '18 at 13:10
  • 1
    $\begingroup$ It would be helpful if you posted the exact link at OEIS. $\endgroup$ – JimB Sep 14 '18 at 13:12
  • 2
    $\begingroup$ @JimB No, that looks like someone took my post and copy-pasted it somewhere... That is highly alarming. $\endgroup$ – Klangen Sep 14 '18 at 13:21
  • 1
    $\begingroup$ The "real values" are 35 and 205 and @AnxonPues gives you the right answer. The number of "pairs" under n might be better given as (Length[Select[ Prime[Range[n]], (PrimeQ[# - 2] || PrimeQ[# + 2]) && # <= n & ]] + 1)/2. It is necessary to add in the condition # <= n. $\endgroup$ – JimB Sep 14 '18 at 13:30
  • 2
    $\begingroup$ I agree. That is more than highly alarming. $\endgroup$ – JimB Sep 14 '18 at 13:32

This I think fails because Prime[Range[n]] gives the first n primes not the primes below n. For this purpose, there might be another function, read help and you will find.

  • $\begingroup$ The code I write in my question is the one on the OEIS, to calculate the twin primes under a positive integer $n$. I doubt it is wrong. $\endgroup$ – Klangen Sep 14 '18 at 12:58

Anxon Pués is absolutely right. What you compute is $\pi_2(p(n))$, so the integral has also to run from $2$ to $\pi_2(p(n))$:

f[n_] := Rest[Accumulate[Subtract[1, Unitize[Differences[Prime[Range[n]]] - 2]]]];
g[n_] = Integrate[Log[t]^(-2), {t, 2, n}, Assumptions -> n >= 2];
n = 100000
ListLinePlot[1.320323632 g@N[Prime[Range[2, n - 1]]]/f[n], 
 PlotRange -> All,
 AxesLabel -> {"p[n]", "Ratio"}]

enter image description here

  • $\begingroup$ Thank you. Then the code in the OEIS is wrong? $\endgroup$ – Klangen Sep 14 '18 at 13:03
  • $\begingroup$ That depends on what the programs are suppose to do. A there is no such program description on OEIS. I guess they are just meant to generate the first couple of elements of the sequence, not the function $\pi_2(n)$. $\endgroup$ – Henrik Schumacher Sep 14 '18 at 13:09
  • $\begingroup$ Alternatively, could you please help me fix it so that it returns $\pi_2(n)$? $\endgroup$ – Klangen Sep 14 '18 at 13:09
  • $\begingroup$ See function f in my post: f[n] returns the list of $\pi_2(2),\pi_2(3), \dotsc,\pi_2(n)$. That's a bit more efficient than searching for the first few twins over and over again. Moreover, it avoids performing prime tests. $\endgroup$ – Henrik Schumacher Sep 14 '18 at 13:11
  • $\begingroup$ Thank you, I understand now. $\endgroup$ – Klangen Sep 14 '18 at 13:23

Not the answer you're looking for? Browse other questions tagged or ask your own question.