Ramanujan primes/twins

Ramanujan primes/twins is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

In a manner similar to twin primes, twin Ramanujan primes may be explored. The task is to determine how many of the first million Ramanujan primes are twins.

Related Task

Twin primes

F#

This task uses Ramanujan primes (F#) <lang fsharp> // Twin Ramanujan primes. Nigel Galloway: September 9th., 2021 printfn $"There are %d{rP 1000000|>Seq.pairwise|>Seq.filter(fun(n,g)->n=g-2)|>Seq.length} twins in the first million Ramanujan primes" </lang>

Output:

There are 74973 twins in the first million Ramanujan primes

Phix

While finding the 1,000,000th Ramanujan prime is reasonably cheap (~70s), repeating that trick to find all 1,000,000 of them individually is over 2 years worth of CPU time. Calculating pi(p) - pi(floor(pi/2) for all primes below that one millionth, and then filtering (normal primes) based on that list is significantly faster. Finally we can than filter/scan that list for twins.

with javascript_semantics
constant lim = 1e5      -- 5.5s
--constant lim = 1e6    -- 1min 11s
atom t0 = time()
sequence picache = {}
function pi(integer n)
    if n=0 then return 0 end if
    if n>length(picache) then
        sequence primes = get_primes_le(n)
        for i=length(picache)+1 to n do
            integer k = binary_search(i,primes)
            if k<0 then k=-k-1 end if
            picache &= k
        end for
    end if
    return picache[n]
end function
 
function Ramanujan_prime(integer n)
    integer maxposs = floor(ceil(4*n*(log(4*n)/log(2))))
    for i=maxposs to 1 by -1 do
        if pi(i) - pi(floor(i/2)) < n then
            return i + 1
        end if
    end for
    return 0
end function

integer rplim = Ramanujan_prime(lim)
printf(1,"The %,dth Ramanujan prime is %,d\n", {lim,rplim})
function rpc(integer p) return pi(p)-pi(floor(p/2)) end function
sequence r = get_primes_le(rplim),
         c = apply(r,rpc)
integer ok = c[$]
for i=length(c)-1 to 1 by -1 do
    if c[i]<ok then
        ok = c[i]
    else
        c[i] = 0
    end if
end for
function nzc(integer p, idx) return c[idx]!=0 end function
r = filter(r,nzc)
integer twins = 0
for i=1 to length(r)-1 do
    if r[i]+2 = r[i+1] then twins += 1 end if
end for
printf(1,"There are %,d twins in the first %,d Ramanujan primes\n", {twins,length(r)})
?elapsed(time()-t0)

Output:

using a smaller limit:

The 100,000th Ramanujan prime is 2,916,539
There are 8,732 twins in the first 100,000 Ramanujan primes
"5.5s"

with the higher limit:

The 1,000,000th Ramanujan prime is 34,072,993
There are 74,973 twins in the first 1,000,000 Ramanujan primes
"1 minute and 11s"