# Find the next twin Mersenne prime pair in a Wilson’s primeth recurrence sequence

Define $$p_{x}$$ to be the $$x$$-th prime number (for example, $$p_{15}=47$$). Then define the recurrence $$a_{0}=1$$ and $$a_{n}=p_{{a_{{n-1}}}} \;\;\text{for}\;\; n>0$$. This is Wilson’s primeth recurrence, which results in the sequence $$W_n =\{1, 2, 3, 5, 11, 31, 127, 709, 5381, 52711, 648391, 9737333, 174440041, \dots\}$$ (A007097 in Sloane’s OEIS). Given the prime counting function $$\pi(x)$$, the recurrence should check out thus: $$\pi(a_{n})=a_{{n-1}}$$.

It suffices to mention Euclid’s proof that there are infinitely many primes to show that this recurrence is also infinite. However, the terms of this recurrence quickly become large enough to show the limitations of today’s computational devices.

Robert G. Wilson provided Sloane with just $$15$$ terms. The last of those was shown to be erroneous by Paul Zimmerman, who was able to extend the known sequence by just two more terms. In 2007, David Baugh discovered two more terms.

Now, let's seek Mersenne prime numbers in the sequence $$W_n$$. Curiously, we find out there are a pair of Mersenne prime numbers that seem to be twins. They are $$31$$ and $$127$$. It is easy to see that $$p_{31} = 127$$, and that both are Mersenne numbers, because $$31 = 2^5 -1$$, where exponent $$5$$ is a prime number, and $$127 = 2^7 -1$$, where exponent $$7$$ is a prime number too. Then, a question arises: what is the next twin Mersenne prime pair in Wilson’s primeth recurrence sequence $$W_n$$, if any?

We are seeking pairs of prime numbers $(p,q)$, with $p<q$, such that:

$$2^p-1 = \pi(2^q -1)$$

and where $M_p = 2^p-1$, and $M_q = 2^q-1$, are Mersenne primes. So, we are seeking prime $p$ as a function of prime $q$;

$$p = \frac{\log(\pi(2^q -1)+1)}{\log 2}$$

or we are seeking $q$ as a function of $p$

$$q = \frac{\log(P(2^p -1)+1)}{\log 2}$$

where $P(n$) is inverse function of the prime counting function $\pi(n)$. If you are trying to find pairs $(p,q)$ of that kind, probably you will reach soon to computational limitations. For example, if you use the famous Mathematica software, maybe you are implementing a routine like this one:

Do[m = 2^Prime[n] - 1; p = Rationalize[N[Log[PrimePi[m] + 1]/Log[2], 1000]]; If[IntegerQ[p], Print[{p, Prime[n], 2^p - 1 == PrimePi[m]}]], {n, 1, 20}]

where Prime[n] is n-th prime number, and where PrimePi[m] is $\pi(m)$. The first limitation, beside long time consuming, is found at the computation of PrimePi[2^Prime[16]] and beyond. In order to solve those severe computational problems, we could use logarithmic integral, instead of $\pi(n)$: $$\rm {li}(x)=\int _{0}^{x}{\frac {dt}{\ln t}} \approx \pi(x)$$

but, we will lose accuracy, and probably our wanted pairs could not be chased. I propose the following function (discovered by me), that approaches $\pi(x)$ better and faster than $\rm {li}(x)$

$$\rm {CL}(x)= \rm {li}(x)-\sqrt[3]{x}+1$$

Now, your computational capacity has been improved, and you can compute beyond the standard limitations, if you use this routine:

CL[n_] := IntegerPart[LogIntegral[n] - Power[n, (3)^-1] + 1];

Do[m = 2^Prime[n] - 1; p = Rationalize[N[Log[CL[m] + 1]/Log[2], 1000]]; If[IntegerQ[p], Print[{p, Prime[n], 2^p - 1 == PrimePi[m]}]], {n, 1, 1000}]

The first three elements found in the output sequence are:

{1,2,False}

{2,3,False}

{5,7,True}

...

this third element, {5,7,True} found, means that $2^7 -1 = 127 = P(2^5 -1)$ is true, so it's one of our wanted pairs.