Ramanujan primes/twins

From Rosetta Code
Task
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



C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <vector>

int32_t ramanujan_maximum(const int32_t& number) {
	return ceil(4 * number * log(4 * number));
}

std::vector<int32_t> initialise_prime_pi(const int32_t& limit) {
	std::vector<int32_t> result(limit, 1);
	result[0] = 0;
	result[1] = 0;
	for ( int32_t i = 4; i < limit; i += 2 ) {
		result[i] = 0;
	}
	for ( int32_t p = 3, square = 9; square < limit; p += 2 ) {
		if ( result[p] != 0 ) {
			for ( int32_t q = square; q < limit; q += p << 1 ) {
				result[q] = 0;
			}
		}
		square += ( p + 1 ) << 2;
	}

	for ( uint64_t i = 1; i < result.size(); ++i ) {
		result[i] += result[i - 1];
	}
	return result;
}

int32_t ramanujan_prime(const std::vector<int32_t>& prime_pi, const int32_t& number) {
	int32_t maximum = ramanujan_maximum(number);
	if ( ( maximum & 1) == 1 ) {
		maximum -= 1;
	}

	int32_t index = maximum;
	while ( prime_pi[index] - prime_pi[index / 2] >= number ) {
		index -= 1;
	}
	return index + 1;
}

std::vector<int32_t> list_primes_less_than(const int32_t& limit) {
	std::vector<bool> composite(limit, false);
	int32_t n = 3;
	int32_t nSquared = 9;
	while ( nSquared <= limit ) {
		if ( ! composite[n] ) {
			for ( int32_t k = nSquared; k < limit; k += 2 * n ) {
				composite[k] = true;
			}
		}
		nSquared += ( n + 1 ) << 2;
		n += 2;
	}

	std::vector<int32_t> result;
	result.emplace_back(2);
	for ( int32_t i = 3; i < limit; i += 2 ) {
		if ( ! composite[i] ) {
			result.emplace_back(i);
		}
	}
	return result;
}

int main() {
	const int32_t limit = 1'000'000;
	const std::vector<int32_t> prime_pi = initialise_prime_pi(ramanujan_maximum(limit) + 1);
	const int32_t millionth_ramanujan_prime = ramanujan_prime(prime_pi, limit);
	std::cout << "The 1_000_000th Ramanujan prime is " << millionth_ramanujan_prime << std::endl;

	std::vector<int32_t> primes = list_primes_less_than(millionth_ramanujan_prime);
	std::vector<int32_t> ramanujan_prime_indexes(primes.size());
	for ( uint64_t i = 0; i < ramanujan_prime_indexes.size(); ++i ) {
		ramanujan_prime_indexes[i] = prime_pi[primes[i]] - prime_pi[primes[i] / 2];
	}

	int32_t lowerLimit = ramanujan_prime_indexes[ramanujan_prime_indexes.size() - 1];
	for ( int32_t i = ramanujan_prime_indexes.size() - 2; i >= 0; --i ) {
		if ( ramanujan_prime_indexes[i] < lowerLimit ) {
			lowerLimit = ramanujan_prime_indexes[i];
		} else {
			ramanujan_prime_indexes[i] = 0;
		}
	}

	std::vector<int32_t> ramanujan_primes;
	for ( uint64_t i = 0; i < ramanujan_prime_indexes.size(); ++i ) {
		if ( ramanujan_prime_indexes[i] != 0 ) {
			ramanujan_primes.emplace_back(primes[i]);
		}
	}

	int32_t twins_count = 0;
	for ( uint64_t i = 0; i < ramanujan_primes.size() - 1; ++i ) {
		if ( ramanujan_primes[i] + 2 == ramanujan_primes[i + 1] ) {
			twins_count++;
		}
	}
	std::cout << "There are " << twins_count << " twins in the first " << limit << " Ramanujan primes." << std::endl;
}
Output:
The 1_000_000th Ramanujan prime is 34072993
There are 74973 twins in the first 1000000 Ramanujan primes.

EasyLang

global cnt[] .
proc primcnt limit . .
   cnt[] = [ 0 1 1 ]
   for i = 4 step 2 to limit
      cnt[] &= 0
      cnt[] &= 1
   .
   p = 3
   sq = 9
   while sq <= limit
      if cnt[p] <> 0
         for q = sq step p * 2 to limit
            cnt[q] = 0
         .
      .
      sq += (p + 1) * 4
      p += 2
   .
   for i = 2 to limit
      sum += cnt[i]
      cnt[i] = sum
   .
.
func log n .
   e = 2.7182818284590452354
   return log10 n / log10 e
.
func ramamax n .
   return floor (4 * n * log (4 * n))
.
func ramaprim_twins n .
   i = ramamax n
   i -= i mod 2
   while i >= 2
      if cnt[i] - cnt[i / 2] < n
         p1 = p
         p = i + 1
         if p1 - p = 2
            cnt += 1
         .
         n -= 1
      .
      i -= 2
   .
   return cnt
.
primcnt ramamax 1000000
print ramaprim_twins 1000000
Output:
74973

F#

This task uses Ramanujan primes (F#)

// 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"
Output:
There are 74973 twins in the first million Ramanujan primes

FreeBASIC

#define floor(x) ((x*2.0-0.5) Shr 1)

Dim Shared pi() As Integer

Sub primeCounter(limit As Integer)
    Dim As Integer i, q, p, sq, total
    Redim pi(limit)
    pi(0) = 0
    pi(1) = 0
    For i = 2 To limit
        pi(i) = 1
    Next
    If limit > 2 Then
        For i = 4 To limit Step 2
            pi(i) = 0
        Next i
        p = 3
        sq = 9
        While sq <= limit
            If pi(p) <> 0 Then
                For q = sq To limit Step p*2
                    pi(q) = 0
                Next q
            End If
            sq += (p + 1) * 4
            p += 2
        Wend
        total = 0
        For i = 2 To limit
            total += pi(i)
            pi(i) = total
        Next i
    End If
End Sub

Function ramanujanMax(n As Integer) As Integer
    Return floor(4 * n * Log(4*n))
End Function

Function ramanujanPrimeTwins(n As Integer) As Integer
    Dim As Integer maxposs, p1, p, cnt
    maxposs = ramanujanMax(n)
    maxposs -= maxposs Mod 2
    While maxposs >= 2
        If pi(maxposs) - pi(maxposs \ 2) < n Then
            p1 = p
            p = maxposs + 1            
            If p1 - p = 2 Then cnt += 1
            n -= 1
        End If    
        maxposs -= 2
    Wend
    Return cnt
End Function

Dim As Integer limit = 1e6
Dim As Double t0 = Timer
primeCounter ramanujanMax(limit)
Print Using "There are & twins in the first & Ramanujan primes"; ramanujanPrimeTwins(limit); limit
Print Using "##.##sec."; Timer - t0

Sleep
Output:
There are 74973 twins in the first 1000000 Ramanujan primes
 1.53sec.

Go

Translation of: Wren
Library: Go-rcu
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()
    }
}
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

J

   _2 +/@:= 2 -/\ (((] - _1&(33 b.)@:>:@[ { ]) _1&p:) i. 35e6) i: i. 1e6
74973

Java

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

public final class RamanujanPrimesTwins {

	public static void main(String[] aArgs) {
		final int limit = 1_000_000;
		final int[] primePi = initialisePrimePi(ramanujanMaximum(limit) + 1);
		final int millionthRamanujanPrime = ramanujanPrime(primePi, limit);
		System.out.println("The 1_000_000th Ramanujan prime is " + millionthRamanujanPrime);
		
		List<Integer> primes = listPrimesLessThan(millionthRamanujanPrime);		
		int[] ramanujanPrimeIndexes = new int[primes.size()];
		for ( int i = 0; i < ramanujanPrimeIndexes.length; i++ ) {
			ramanujanPrimeIndexes[i] = primePi[primes.get(i)] - primePi[primes.get(i) / 2];
		}
		int lowerLimit = ramanujanPrimeIndexes[ramanujanPrimeIndexes.length - 1];
		for ( int i = ramanujanPrimeIndexes.length - 2; i >= 0; i-- ) {
			if ( ramanujanPrimeIndexes[i] < lowerLimit ) {
				lowerLimit = ramanujanPrimeIndexes[i];
			} else {
				ramanujanPrimeIndexes[i] = 0;
			}
		}
		
		List<Integer> ramanujanPrimes = new ArrayList<Integer>();
		for ( int i = 0; i < ramanujanPrimeIndexes.length; i++ ) {
			if ( ramanujanPrimeIndexes[i] != 0 ) {
				ramanujanPrimes.add(primes.get(i));
			}
		}
		
		int twinsCount = 0;
        for ( int i = 0; i < ramanujanPrimes.size() - 1; i++ ) {
            if ( ramanujanPrimes.get(i) + 2 == ramanujanPrimes.get(i + 1) ) {
                twinsCount += 1;
            }
        }
		System.out.println("There are " + twinsCount + " twins in the first " + limit + " Ramanujan primes.");
	}
	
	private static List<Integer> listPrimesLessThan(int aLimit) {
		boolean[] composite = new boolean[aLimit + 1];
		int n = 3;
		int nSquared = 9;
		while ( nSquared <= aLimit ) { 
		    if ( ! composite[n] ) {
		    	for ( int k = nSquared; k < aLimit; k += 2 * n ) {
		    		composite[k] = true;
		    	}
		    }
		    nSquared += ( n + 1 ) << 2;
		    n += 2;
		}
		
		List<Integer> result = new ArrayList<Integer>();
		result.add(2);
		for ( int i = 3; i < aLimit; i += 2 ) {
		    if ( ! composite[i] ) {
		    	result.add(i);
		    }
		}
		return result;
	}

	private static int ramanujanPrime(int[] aPrimePi, int aNumber) {
		int maximum = ramanujanMaximum(aNumber);
		if ( ( maximum & 1) == 1 ) {
			maximum -= 1;
		}
		
		int index = maximum;
		while ( aPrimePi[index] - aPrimePi[index / 2] >= aNumber ) {
			index -= 1;
		}
		return index + 1;
	}
	
	private static int[] initialisePrimePi(int aLimit) {
		int[] result = new int[aLimit];
        Arrays.fill(result, 1);
        result[0] = 0;
        result[1] = 0;
        for ( int i = 4; i < aLimit; i += 2 ) {
            result[i] = 0;
        }
        for ( int p = 3, square = 9; square < aLimit; p += 2 ) {
            if ( result[p] != 0 ) {
                for ( int q = square; q < aLimit; q += p << 1 ) {
                    result[q] = 0;
                }
            }
            square += ( p + 1 ) << 2;
        }
        Arrays.parallelPrefix(result, Integer::sum);
        return result;
    }
	
	private static int ramanujanMaximum(int aNumber) {
		return (int) Math.ceil(4 * aNumber * Math.log(4 * aNumber));
	}

}
Output:
The 1_000_000th Ramanujan prime is 34072993
There are 74973 twins in the first 1000000 Ramanujan primes.

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)
Output:
There are 74973 twins in the first 1000000 Ramanujan primes.
  0.759511 seconds (50 allocations: 839.383 MiB, 5.64% gc time)

Mathematica/Wolfram Language

$HistoryLength = 1;
l = PrimePi[Range[35 10^6]] - PrimePi[Range[35 10^6]/2];
ramanujanprimes = GatherBy[Transpose[{Range[2, Length[l] + 1], l}], Last][[All, -1, 1]];
ramanujanprimes = Take[Sort@ramanujanprimes, 10^6];
Count[Differences[ramanujanprimes], 2]
Output:
74973

Nim

Translation of: Go

We use Phix algorithm in its Go translation.

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

Perl

Can't do better than the ntheory module.

Library: ntheory
use strict;
use warnings;
use ntheory <ramanujan_primes nth_ramanujan_prime>;

sub comma { reverse ((reverse shift) =~ s/(.{3})/$1,/gr) =~ s/^,//r }

my $rp = ramanujan_primes nth_ramanujan_prime 1_000_000;
for my $limit (1e5, 1e6) {
    printf "The %sth Ramanujan prime is %s\n", comma($limit), comma $rp->[$limit-1];
    printf  "There are %s twins in the first %s Ramanujan primes.\n\n",
        comma( scalar grep { $rp->[$_+1] == $rp->[$_]+2 } 0 .. $limit-2 ), comma $limit;
}
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.

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.

Library: ntheory
use ntheory:from<Perl5> <ramanujan_primes nth_ramanujan_prime>;
use Lingua::EN::Numbers;

my @rp = ramanujan_primes nth_ramanujan_prime 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';
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.529 seconds

Scala

Translation of: Java
import java.util.{ArrayList, Arrays, List}

object RamanujanPrimesTwins {

  def main(args: Array[String]): Unit = {
    val limit = 1_000_000
    val primePi = initialisePrimePi(ramanujanMaximum(limit) + 1)
    val millionthRamanujanPrime = ramanujanPrime(primePi, limit)
    println(s"The 1_000_000th Ramanujan prime is $millionthRamanujanPrime")

    val primes = listPrimesLessThan(millionthRamanujanPrime)
    val ramanujanPrimeIndexes = new Array[Int](primes.size)
    for (i <- 0 until primes.size by 1) {
      ramanujanPrimeIndexes(i) = primePi(primes.get(i)) - primePi(primes.get(i) / 2)
    }

    var lowerLimit = ramanujanPrimeIndexes(ramanujanPrimeIndexes.length - 1)
    for (i <- ramanujanPrimeIndexes.length - 2 to 0 by -1) {
      if (ramanujanPrimeIndexes(i) < lowerLimit) {
        lowerLimit = ramanujanPrimeIndexes(i)
      } else {
        ramanujanPrimeIndexes(i) = 0
      }
    }

    val ramanujanPrimes = new ArrayList[Integer]()
    for (i <- 0 until ramanujanPrimeIndexes.size by 1) {
      if (ramanujanPrimeIndexes(i) != 0) {
        ramanujanPrimes.add(primes.get(i))
      }
    }

    var twinsCount = 0
    for (i <- 0 until ramanujanPrimes.size - 1) {
      if (ramanujanPrimes.get(i) + 2 == ramanujanPrimes.get(i + 1)) {
        twinsCount += 1
      }
    }
    println(s"There are $twinsCount twins in the first $limit Ramanujan primes.")
  }

  private def listPrimesLessThan(limit: Int): List[Integer] = {
    val composite = new Array[Boolean](limit + 1)
    var n = 3
    var nSquared = 9
    while (nSquared <= limit) {
      if (!composite(n)) {
        for (k <- nSquared until limit by (2*n) ) {
          composite(k) = true
        }
      }
      nSquared += (n + 1) << 2
      n += 2
    }

    var result = new ArrayList[Integer]()
    result.add(2)
    for (i <- 3 until limit by 2) {
      if (!composite(i)) {
        result.add(i)
      }
    }
    result
  }

  private def ramanujanPrime(primePi: Array[Int], number: Int): Int = {
    var maximum = ramanujanMaximum(number)
    if ((maximum & 1) == 1) {
      maximum -= 1
    }

    var index = maximum
    while (primePi(index) - primePi(index / 2) >= number) {
      index -= 1
    }
    index + 1
  }

  private def initialisePrimePi(aLimit : Int): Array[Int] = {
    val result = new Array[Int](aLimit )
    Arrays.fill(result, 1)
    result(0) = 0
    result(1) = 0
    for (i <- 4 until aLimit by 2) {
      result(i) = 0
    }

    var p = 3;var square = 9;
    while(square < aLimit) {
            if ( result(p) != 0 ) {
                for ( q <- square until aLimit by (p << 1) ) {
                    result(q) = 0;
                }
            }
            square += ( p + 1 ) << 2;
            p += 2
    }
    for (i <- 1 until result.length) {
      result(i) += result(i - 1)
    }
    result
  }

  private def ramanujanMaximum(number: Int): Int = {
    Math.ceil(4 * number * Math.log(4 * number)).toInt
  }
}
Output:
The 1_000_000th Ramanujan prime is 34072993
There are 74973 twins in the first 1000000 Ramanujan primes.

Sidef

Library: ntheory
require('ntheory')

var rp = %S<ntheory>.ramanujan_primes(%S<ntheory>.nth_ramanujan_prime(1e6))
for limit in ([1e5, 1e6]) {
    printf("The %sth Ramanujan prime is %s\n", limit.commify, rp[limit-1].commify)
    printf("There are %s twins in the first %s Ramanujan primes.\n\n",
        ^(limit-1) -> count {|i| rp[i+1] == rp[i]+2 }.commify, limit.commify)
}
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.

Wren

Translation of: Phix
Library: Wren-iterate
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.

import "./iterate" 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")
}
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.