Ramanujan primes/twins
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 (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 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"
Wren
As calculating the first 1 million Ramanujan primes individually is out of the question, we follow the Phix entry's clever strategy here of which this is a partial translation.
The only reason this is faster than Phix (normally we'd be much slower) is because the basic Ramanujan code is a translation of the very fast C++ code for that task. <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.