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
with javascript_semantics
sequence pi = {}
 
procedure primeCounter(integer limit)
    pi = repeat(1,limit)
    if limit > 1 then
        pi[1] = 0
        for i=4 to limit by 2 do
            pi[i] = 0
        end for
        integer p = 3, sq = 9
        while sq<=limit do
            if pi[p]!=0 then
                for q=sq to limit by p*2 do
                    pi[q] = 0
                end for
            end if
            sq += (p + 1)*4
            p += 2
        end while
        atom total = 0
        for i=2 to limit do
            total += pi[i]
            pi[i] = total
        end for
    end if
end procedure
 
function ramanujanMax(integer n)
    return floor(4*n*log(4*n))
end function
 
function ramanujanPrime(integer n)
    if n=1 then return 2 end if
    integer maxposs = ramanujanMax(n)
    for i=maxposs-odd(maxposs) to 1 by -2 do
        if pi[i]-pi[floor(i/2)] < n then
            return i + 1
        end if
    end for
    return 0
end function

constant lim = 1e5      -- 0.6s
--constant lim = 1e6    -- 4.7s -- (not pwa/p2js)
atom t0 = time()
primeCounter(ramanujanMax(lim))
integer rplim = ramanujanPrime(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)
