Ramanujan primes/twins

Ramanujan primes/twins
You are encouraged to solve this task according to the task description, using any language you may know.

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

Go

Translation of: Wren

Library: Go-rcu

<lang go>package main

import (

   "fmt"
   "math"
   "rcu"
   "time"

)

var count []int

func primeCounter(limit int) {

   count = make([]int, limit)
   for i := 0; i < limit; i++ {
       count[i] = 1
   }
   if limit > 0 {
       count[0] = 0
   }
   if limit > 1 {
       count[1] = 0
   }
   for i := 4; i < limit; i += 2 {
       count[i] = 0
   }
   for p, sq := 3, 9; sq < limit; p += 2 {
       if count[p] != 0 {
           for q := sq; q < limit; q += p << 1 {
               count[q] = 0
           }
       }
       sq += (p + 1) << 2
   }
   sum := 0
   for i := 0; i < limit; i++ {
       sum += count[i]
       count[i] = sum
   }

}

func primeCount(n int) int {

   if n < 1 {
       return 0
   }
   return count[n]

}

func ramanujanMax(n int) int {

   fn := float64(n)
   return int(math.Ceil(4 * fn * math.Log(4*fn)))

}

func ramanujanPrime(n int) int {

   if n == 1 {
       return 2
   }
   for i := ramanujanMax(n); i >= 2*n; i-- {
       if i%2 == 1 {
           continue
       }
       if primeCount(i)-primeCount(i/2) < n {
           return i + 1
       }
   }
   return 0

}

func rpc(p int) int { return primeCount(p) - primeCount(p/2) }

func main() {

   for _, limit := range []int{1e5, 1e6} {
       start := time.Now()
       primeCounter(1 + ramanujanMax(limit))
       rplim := ramanujanPrime(limit)
       climit := rcu.Commatize(limit)
       fmt.Printf("The %sth Ramanujan prime is %s\n", climit, rcu.Commatize(rplim))
       r := rcu.Primes(rplim)
       c := make([]int, len(r))
       for i := 0; i < len(c); i++ {
           c[i] = rpc(r[i])
       }
       ok := c[len(c)-1]
       for i := len(c) - 2; i >= 0; i-- {
           if c[i] < ok {
               ok = c[i]
           } else {
               c[i] = 0
           }
       }
       var fr []int
       for i, r := range r {
           if c[i] != 0 {
               fr = append(fr, r)
           }
       }
       twins := 0
       for i := 0; i < len(fr)-1; i++ {
           if fr[i]+2 == fr[i+1] {
               twins++
           }
       }
       fmt.Printf("There are %s twins in the first %s Ramanujan primes.\n", rcu.Commatize(twins), climit)
       fmt.Println("Took", time.Since(start))
       fmt.Println()
   }

}</lang>

Output:

The 100,000th Ramanujan prime is 2,916,539
There are 8,732 twins in the first 100,000 Ramanujan primes.
Took 50.745239ms

The 1,000,000th Ramanujan prime is 34,072,993
There are 74,973 twins in the first 1,000,000 Ramanujan primes.
Took 699.967994ms

Julia

<lang julia>using Primes

function rajpairs(N, verbose, interval = 2)

   maxpossramanujan(n) = Int(ceil(4n * log(4n) / log(2)))
   prm = primes(maxpossramanujan(N))
   pivec = accumulate(+, primesmask(maxpossramanujan(N)))
   halfrpc = [pivec[p] - pivec[p ÷ 2] for p in prm]
   lastrpc = last(halfrpc) + 1
   for i in length(halfrpc):-1:1
       if halfrpc[i] < lastrpc
           lastrpc = halfrpc[i]
       else
           halfrpc[i] = 0
       end
   end
   rajvec = [prm[i] for i in eachindex(prm) if halfrpc[i] > 0]
   nrajtwins = count(rajvec[i] + interval == rajvec[i + 1] for i in 1:N-1)
   verbose && println("There are $nrajtwins twins in the first $N Ramanujan primes.")
   return nrajtwins

end

rajpairs(100000, false) @time rajpairs(1000000, true)

</lang>

Output:

There are 74973 twins in the first 1000000 Ramanujan primes.
  0.759511 seconds (50 allocations: 839.383 MiB, 5.64% gc time)

Nim

Translation of: Go

We use Phix algorithm in its Go translation. <lang Nim>import math, sequtils, strutils, sugar, times

let t0 = now()

type PrimeCounter = seq[int32]

proc initPrimeCounter(limit: Positive): PrimeCounter {.noInit.} =

 doAssert limit > 1
 result = repeat(1i32, limit)
 result[0] = 0
 result[1] = 0
 for i in countup(4, limit - 1, 2): result[i] = 0
 var p = 3
 var p2 = 9
 while p2 < limit:
   if result[p] != 0:
     for q in countup(p2, limit - 1, p + p):
       result[q] = 0
   inc p, 2
   p2 = p * p
 # Compute partial sums in place.
 var sum = 0i32
 for item in result.mitems:
   sum += item
   item = sum

func ramanujanMax(n: int): int {.inline.} = int(ceil(4 * n.toFloat * ln(4 * n.toFloat)))

func ramanujanPrime(pi: PrimeCounter; n: int): int =

 if n == 1: return 2
 var max = ramanujanMax(n)
 if (max and 1) == 1: dec max
 for i in countdown(max, 2, 2):
   if pi[i] - pi[i div 2] < n:
     return i + 1

func primesLe(limit: Positive): seq[int] =

 var composite = newSeq[bool](limit + 1)
 var n = 3
 var n2 = 9
 while n2 <= limit:
   if not composite[n]:
     for k in countup(n2, limit, 2 * n):
       composite[k] = true
   n2 += (n + 1) shl 2
   n += 2
 result = @[2]
 for n in countup(3, limit, 2):
   if not composite[n]: result.add n

proc main() =

 const Lim = 1_000_000
 let pi = initPrimeCounter(1 + ramanujanMax(Lim))
 let rpLim = ramanujanPrime(pi, Lim)
 echo "The 1_000_000th Ramanujan prime is $#.".format(($rpLim).insertSep())
 let r = primesLe(rpLim)
 var c = r.mapIt(pi[it] - pi[it div 2])
 var ok = c[^1]
 for i in countdown(c.len - 2, 0):
   if c[i] < ok:
     ok = c[i]
   else:
     c[i] = 0
 let fr = collect(newSeq, for i, p in r: (if c[i] != 0: p))
 var twins = 0
 var prev = -1
 for p in fr:
   if p == prev + 2: inc twins
   prev = p
 echo "There are $1 twins in the first $2 Ramanujan primes.".format(($twins).insertSep(), ($Lim).insertSep)

main() echo "\nElapsed time: ", (now() - t0).inMilliseconds, " ms"</lang>

Output:

The 1_000_000th Ramanujan prime is 34_072_993.
There are 74_973 twins in the first 1_000_000 Ramanujan primes.

Elapsed time: 1139 ms

Phix

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

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)

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
"0.6s"

with the higher limit (desktop/Phix only, alas not possible under pwa/p2js):

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

Raku

Timing is informational only. Will vary by system specs and load. <lang perl6>use ntheory:from<Perl5> <ramanujan_primes>; use Lingua::EN::Numbers;

my @rp = ramanujan_primes(100_000_000)[^1_000_000];

for (1e5, 1e6)».Int -> $limit {

   say "\nThe {comma $limit}th Ramanujan prime is { comma @rp[$limit-1]}";
   say "There are { comma +(^($limit-1)).race.grep: { @rp[$_+1] == @rp[$_]+2 } } twins in the first {comma $limit} Ramanujan primes.";

}

say (now - INIT now).fmt('%.3f') ~ ' seconds';</lang>

Output:

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

The 1,000,000th Ramanujan prime is 34,072,993
There are 74,973 twins in the first 1,000,000 Ramanujan primes.
2.863 seconds

Wren

Translation of: Phix

Library: Wren-trait

Library: Wren-math

Library: Wren-fmt

As calculating the first 1 million Ramanujan primes individually is out of the question, we follow the Phix entry's clever strategy here which brings this task comfortably within our ambit. <lang ecmascript>import "/trait" for Stepped, Indexed import "/math" for Int, Math import "/fmt" for Fmt

var count

var primeCounter = Fn.new { |limit|

   count = List.filled(limit, 1)
   if (limit > 0) count[0] = 0
   if (limit > 1) count[1] = 0
   for (i in Stepped.new(4...limit, 2)) count[i] = 0
   var p = 3
   var sq = 9
   while (sq < limit) {
       if (count[p] != 0) {
           var q = sq
           while (q < limit) {
               count[q] = 0
               q = q + p * 2
           }
       }
       sq = sq + (p + 1) * 4
       p = p + 2
   }
   var sum = 0
   for (i in 0...limit) {
       sum = sum + count[i]
       count[i] = sum
   }

}

var primeCount = Fn.new { |n| (n < 1) ? 0 : count[n] }

var ramanujanMax = Fn.new { |n| (4 * n * (4*n).log).ceil }

var ramanujanPrime = Fn.new { |n|

   if (n == 1) return 2
   for (i in ramanujanMax.call(n)..2) {
       if (i % 2 == 1) continue
       if (primeCount.call(i) - primeCount.call((i/2).floor) < n) return i + 1
   }
   return 0

}

var rpc = Fn.new { |p| primeCount.call(p) - primeCount.call((p/2).floor) }

for (limit in [1e5, 1e6]) {

   var start = System.clock
   primeCounter.call(1 + ramanujanMax.call(limit))
   var rplim = ramanujanPrime.call(limit)
   Fmt.print("The $,dth Ramanujan prime is $,d", limit, rplim)
   var r = Int.primeSieve(rplim)
   var c = r.map { |p| rpc.call(p) }.toList
   var ok = c[-1]
   for (i in c.count - 2..0) {
       if (c[i] < ok) {
           ok = c[i]
       } else {
           c[i] = 0
       }
   }
   var ir = Indexed.new(r).where { |se| c[se.index] != 0 }.toList
   var twins = 0
   for (i in 0...ir.count-1) {
       if (ir[i].value + 2 == ir[i+1].value) twins = twins + 1
   }
   Fmt.print("There are $,d twins in the first $,d Ramanujan primes.", twins, limit)
   System.print("Took %(Math.toPlaces(System.clock -start, 2, 0)) seconds.\n")

}</lang>

Output:

The 100,000th Ramanujan prime is 2,916,539
There are 8,732 twins in the first 100,000 Ramanujan primes.
Took 1.33 seconds.

The 1,000,000th Ramanujan prime is 34,072,993
There are 74,973 twins in the first 1,000,000 Ramanujan primes.
Took 15.75 seconds.