Ramanujan primes/twins: Difference between revisions
SqrtNegInf (talk | contribs) m (→{{header|Raku}}: note use of 'ntheory' module) |
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
||
Line 4: | Line 4: | ||
=={{header|F_Sharp|F#}}== |
=={{header|F_Sharp|F#}}== |
||
This task uses [http://www.rosettacode.org/wiki/Ramanujan_primes#F.23 Ramanujan primes (F#)] |
This task uses [http://www.rosettacode.org/wiki/Ramanujan_primes#F.23 Ramanujan primes (F#)] |
||
< |
<syntaxhighlight lang="fsharp"> |
||
// Twin Ramanujan primes. Nigel Galloway: September 9th., 2021 |
// 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" |
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" |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 16: | Line 16: | ||
{{trans|Wren}} |
{{trans|Wren}} |
||
{{libheader|Go-rcu}} |
{{libheader|Go-rcu}} |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 121: | Line 121: | ||
fmt.Println() |
fmt.Println() |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 135: | Line 135: | ||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
< |
<syntaxhighlight lang="julia">using Primes |
||
function rajpairs(N, verbose, interval = 2) |
function rajpairs(N, verbose, interval = 2) |
||
Line 158: | Line 158: | ||
rajpairs(100000, false) |
rajpairs(100000, false) |
||
@time rajpairs(1000000, true) |
@time rajpairs(1000000, true) |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
There are 74973 twins in the first 1000000 Ramanujan primes. |
There are 74973 twins in the first 1000000 Ramanujan primes. |
||
Line 165: | Line 165: | ||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
||
< |
<syntaxhighlight lang="mathematica">$HistoryLength = 1; |
||
l = PrimePi[Range[35 10^6]] - PrimePi[Range[35 10^6]/2]; |
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 = GatherBy[Transpose[{Range[2, Length[l] + 1], l}], Last][[All, -1, 1]]; |
||
ramanujanprimes = Take[Sort@ramanujanprimes, 10^6]; |
ramanujanprimes = Take[Sort@ramanujanprimes, 10^6]; |
||
Count[Differences[ramanujanprimes], 2]</ |
Count[Differences[ramanujanprimes], 2]</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>74973</pre> |
<pre>74973</pre> |
||
Line 176: | Line 176: | ||
{{trans|Go}} |
{{trans|Go}} |
||
We use Phix algorithm in its Go translation. |
We use Phix algorithm in its Go translation. |
||
< |
<syntaxhighlight lang="nim">import math, sequtils, strutils, sugar, times |
||
let t0 = now() |
let t0 = now() |
||
Line 248: | Line 248: | ||
main() |
main() |
||
echo "\nElapsed time: ", (now() - t0).inMilliseconds, " ms"</ |
echo "\nElapsed time: ", (now() - t0).inMilliseconds, " ms"</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 259: | Line 259: | ||
Can't do better than the <code>ntheory</code> module. |
Can't do better than the <code>ntheory</code> module. |
||
{{libheader|ntheory}} |
{{libheader|ntheory}} |
||
< |
<syntaxhighlight lang="perl">use strict; |
||
use warnings; |
use warnings; |
||
use ntheory <ramanujan_primes nth_ramanujan_prime>; |
use ntheory <ramanujan_primes nth_ramanujan_prime>; |
||
Line 270: | Line 270: | ||
printf "There are %s twins in the first %s Ramanujan primes.\n\n", |
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; |
comma( scalar grep { $rp->[$_+1] == $rp->[$_]+2 } 0 .. $limit-2 ), comma $limit; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>The 100,000th Ramanujan prime is 2,916,539 |
<pre>The 100,000th Ramanujan prime is 2,916,539 |
||
Line 283: | Line 283: | ||
Calculating pi(p) - pi(floor(pi/2) for all (normal) primes below that one millionth, and then |
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. |
filtering them based on that list is significantly faster, and finally we can scan for twins. |
||
<!--< |
<!--<syntaxhighlight lang="phix">(phixonline)--> |
||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pi</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> |
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pi</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> |
||
Line 352: | Line 352: | ||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"There are %,d twins in the first %,d Ramanujan primes\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">twins</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)})</span> |
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"There are %,d twins in the first %,d Ramanujan primes\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">twins</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)})</span> |
||
<span style="color: #0000FF;">?</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)</span> |
<span style="color: #0000FF;">?</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)</span> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
{{out}} |
{{out}} |
||
using a smaller limit: |
using a smaller limit: |
||
Line 370: | Line 370: | ||
Timing is informational only. Will vary by system specs and load. |
Timing is informational only. Will vary by system specs and load. |
||
{{libheader|ntheory}} |
{{libheader|ntheory}} |
||
<lang |
<syntaxhighlight lang="raku" line>use ntheory:from<Perl5> <ramanujan_primes nth_ramanujan_prime>; |
||
use Lingua::EN::Numbers; |
use Lingua::EN::Numbers; |
||
Line 380: | Line 380: | ||
} |
} |
||
say (now - INIT now).fmt('%.3f') ~ ' seconds';</ |
say (now - INIT now).fmt('%.3f') ~ ' seconds';</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>The 100,000th Ramanujan prime is 2,916,539 |
<pre>The 100,000th Ramanujan prime is 2,916,539 |
||
Line 396: | Line 396: | ||
<br> |
<br> |
||
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. |
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. |
||
< |
<syntaxhighlight lang="ecmascript">import "/trait" for Stepped, Indexed |
||
import "/math" for Int, Math |
import "/math" for Int, Math |
||
import "/fmt" for Fmt |
import "/fmt" for Fmt |
||
Line 464: | Line 464: | ||
Fmt.print("There are $,d twins in the first $,d Ramanujan primes.", twins, limit) |
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") |
System.print("Took %(Math.toPlaces(System.clock -start, 2, 0)) seconds.\n") |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
Revision as of 11:48, 28 August 2022
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#)
// 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
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()
}
}
- 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
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
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.
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.
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
Wren
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 "/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")
}
- 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.