2

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?

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309
Klangen
  • 1,009
  • 5
  • 14
  • 2
    You've posted this in a blog already? http://mathematics.filegala.com/2018/09/14/asymptotic-density-of-twin-primes-gives-wrong-result-in-mathematica/ – JimB Sep 14 '18 at 13:10
  • 1
    It would be helpful if you posted the exact link at OEIS. – JimB Sep 14 '18 at 13:12
  • 2
    @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
  • 1
    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
  • 2
    I agree. That is more than highly alarming. – JimB Sep 14 '18 at 13:32
  • @JimB Posted on Meta: https://math.meta.stackexchange.com/questions/29098/it-looks-like-this-website-is-copying-mathematica-se-questions-automatically – Klangen Sep 14 '18 at 13:33
  • @JimB Why not simply (Length[Select[ Prime[Range[n]], (PrimeQ[# + 2]) && # <= n & ]] + 1)? – Klangen Sep 14 '18 at 13:44
  • Also: https://meta.stackexchange.com/questions/200177/a-site-or-scraper-is-copying-content-from-stack-exchange-what-should-i-do – Klangen Sep 14 '18 at 13:46
  • 1
    I think you want p2[n_] := Length[Select[Prime[Range[PrimePi[n]]], PrimeQ[# + 2] &]] to get the upper bound correct. – Daniel Lichtblau Sep 14 '18 at 14:54
  • @PreservedFruit. Yes, your function better. I just wanted to add in the # < = n part to your original equation as the item that gives the result you desired. – JimB Sep 14 '18 at 15:14
  • 1
    The answers to this question give further insight and faster methods to count twin primes. – KennyColnago Sep 14 '18 at 15:14

2 Answers2

5

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.

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309
Anxon Pués
  • 907
  • 5
  • 13
  • 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
3

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

enter image description here

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