[arhynpov] Statistical infinitude of twin primes

Assuming that primes are distributed randomly with density 1/log(n), we can estimate the number of primes with Mathematica:

In[41]:= Integrate[1/Log[x],x]
Out[41]= LogIntegral[x]

How many total primes are there, probably?

In[149]:= Limit[Integrate[1/Log[x],{x,2,n}],n -> Infinity]
Out[149]= Infinity

(Is this the worst "proof" of there being an infinite number of primes?  Among other terrible things, we're using the Prime Number Theorem as a lemma.)

Similarly, we can estimate the number of twin primes.  We assume that each number is an independent random sample, so you just need to multiply probabilities to get the probability that two numbers about the same size are both prime.

In[43]:= InputForm[Integrate[1/Log[x]^2,x]]
Out[43]//InputForm= -(x/Log[x]) + LogIntegral[x]

Both terms are O(x/ln x), according to Wikipedia.  I was unable to coax the output into a form that makes it clear the asymptotic behavior as x tends to infinity.  However, Mathematica can evaluate the limit:

In[150]:= Limit[Integrate[1/Log[x]^2,{x,2,n}],n -> Infinity]
Out[150]= Infinity

We can do the same with prime triplets and beyond:

In[106]:= Integrate[1/Log[x]^3,x]//InputForm
Out[106]//InputForm= -x/(2*Log[x]^2) - x/(2*Log[x]) + LogIntegral[x]/2

In[151]:= Limit[Integrate[1/Log[x]^3,{x,2,n}],n -> Infinity]
Out[151]= Infinity

Therefore, if the Twin Prime Conjecture is true, it would not be statistically surprising.  It would be very surprising if it is false.  Similarly -- more astoundingly -- we statistically expect there to exist an infinite number of any admissible prime constellation, even though for some of them -- many of them, not a single example is known.  Extremely dense clusters of primes keep occuring (over and over again) even very far out on the number line where primes are very sparse.

This seems vaguely philosophically similar to the Boltzmann Brain.

* * *

Some Mathematica version 11.3.0 "out-takes", illustrating it's easy to push this software past its limits:

In[97]:= Integrate[1/Log[x]^5,{x,E,2^64}]

Throw::sysexc: Uncaught SystemException returned to top level. Can be caught with Catch[…, _SystemException].

Out[97]= SystemException[MemoryAllocationFailure]

but the following work

In[1]:= Integrate[1/Log[x]^5,{x,E,2^64+1}] // InputForm
Out[1]//InputForm= (10*E - (18446744073709551617*(6 + Log[18446744073709551617]* (2 + Log[18446744073709551617] + Log[18446744073709551617]^2)))/ Log[18446744073709551617]^4 + LogIntegral[18446744073709551617] - LogIntegral[E])/24

In[2]:= Integrate[1/Log[x]^5,{x,E,2^64-1}] // InputForm
Out[2]//InputForm= (10*E - (18446744073709551615*(6 + Log[18446744073709551615]* (2 + Log[18446744073709551615] + Log[18446744073709551615]^2)))/ Log[18446744073709551615]^4 + LogIntegral[18446744073709551615] - LogIntegral[E])/24

Another bug:

In[1]:= N[Integrate[1/Log[x]^13,{x,E,2^126}]] // InputForm
Out[1]//InputForm= -8.318336974075375*^12

The integrand is strictly positive, so there's no way the result should be negative.  Note that the output is Mathematica's somewhat bizarre notation for -8.318336974075375e+12.

In[2]:= $MachinePrecision
Out[2]= 15.9546

In[3]:= N[Integrate[1/Log[x]^13,{x,E,2^126}],MachinePrecision] // InputForm
Out[3]//InputForm= -8.318336974075375*^12

In[4]:= N[Integrate[1/Log[x]^13,{x,E,2^126}],$MachinePrecision] // InputForm
Out[4]//InputForm= 5.8248146183321134960348204749329`15.954589770191003*^12

In[5]:= N[Integrate[1/Log[x]^13,{x,E,2^126}],15] // InputForm
Out[5]//InputForm= 5.8248146183321134960348204749329`15.*^12

also

In[1]:= N[Integrate[1/Log[x]^13,{x,E,2^126}],1] // InputForm

                                       1
  Divide::infy: Infinite expression ------- encountered.
                                         24
                                    0. 10

Out[1]//InputForm= 5.824814618332`1.*^12

The warning goes away if you evaluate the same expression again.  It's also annoying you can't get non-pretty printed error messages.

In[2]:= N[Integrate[1/Log[x]^13,{x,E,2^126}],1] // InputForm
Out[2]//InputForm= 5.824814618332`1.*^12

One can also suppress the warning if you first evaluate it to a higher precision (e.g., 20) then try with "1".