# Asymptotic density of twin primes gives wrong result in Mathematica [closed]

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

N[Integrate[Log[t]^(-2),{t,2,n}]]*1.320323632


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?

• You've posted this in a blog already? mathematics.filegala.com/2018/09/14/… – JimB Sep 14 '18 at 13:10
• It would be helpful if you posted the exact link at OEIS. – JimB Sep 14 '18 at 13:12
• @JimB No, that looks like someone took my post and copy-pasted it somewhere... That is highly alarming. – Klangen Sep 14 '18 at 13:21
• 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. – JimB Sep 14 '18 at 13:30
• I agree. That is more than highly alarming. – 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.

• 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. – 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))$:

ClearAll[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"}] • Thank you. Then the code in the OEIS is wrong? – Klangen Sep 14 '18 at 13:03
• 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)$. – Henrik Schumacher Sep 14 '18 at 13:09
• Alternatively, could you please help me fix it so that it returns $\pi_2(n)$? – Klangen Sep 14 '18 at 13:09
• 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. – Henrik Schumacher Sep 14 '18 at 13:11
• Thank you, I understand now. – Klangen Sep 14 '18 at 13:23