Ramanujan primes
As the integers get larger, the spacing between prime numbers slowly lengthens, but the spacing between primes increases at a slower rate than the numbers themselves increase. A consequence of this difference in rates of increase is the existence of special primes, called Ramanujan primes.
The `n`th Ramanujan prime is defined to be the least integer for which there are at least n primes between x and x/2 for all x greater or equal to n.
- Task
- Generate and show the first 100 Ramanujan prime numbers.
- Find and show the 1000th Ramanujan prime number.
- Stretch task
- Find and show the 10,000th Ramanujan prime number.
- See also
- The pi prime function, not to be confused with the transcendental number π
- The OEIS entry: OEIS entry
- The Wikipedia entry: Ramanujan_prime.
C++
<lang cpp>#include <chrono>
- include <cmath>
- include <iomanip>
- include <iostream>
- include <numeric>
- include <vector>
class prime_counter { public:
explicit prime_counter(int limit);
int prime_count(int n) const { return n < 1 ? 0 : count_.at(n); }
private:
std::vector<int> count_;
};
prime_counter::prime_counter(int limit) : count_(limit, 1) {
if (limit > 0)
count_[0] = 0;
if (limit > 1)
count_[1] = 0;
for (int i = 4; i < limit; i += 2)
count_[i] = 0;
for (int p = 3, sq = 9; sq < limit; p += 2) {
if (count_[p]) {
for (int q = sq; q < limit; q += p << 1)
count_[q] = 0;
}
sq += (p + 1) << 2;
}
std::partial_sum(count_.begin(), count_.end(), count_.begin());
}
int ramanujan_max(int n) {
return static_cast<int>(std::ceil(4 * n * std::log(4 * n)));
}
int ramanujan_prime(const prime_counter& pc, int n) {
int max = ramanujan_max(n);
for (int i = max; i >= 0; --i) {
if (pc.prime_count(i) - pc.prime_count(i / 2) < n)
return i + 1;
}
return 0;
}
int main() {
std::cout.imbue(std::locale(""));
auto start = std::chrono::high_resolution_clock::now();
prime_counter pc(1 + ramanujan_max(100000));
for (int i = 1; i <= 100; ++i) {
std::cout << std::setw(5) << ramanujan_prime(pc, i)
<< (i % 10 == 0 ? '\n' : ' ');
}
std::cout << '\n';
for (int n = 1000; n <= 100000; n *= 10) {
std::cout << "The " << n << "th Ramanujan prime is " << ramanujan_prime(pc, n)
<< ".\n";
}
auto end = std::chrono::high_resolution_clock::now();
std::cout << "\nElapsed time: "
<< std::chrono::duration<double>(end - start).count() * 1000
<< " milliseconds\n";
}</lang>
- Output:
2 11 17 29 41 47 59 67 71 97 101 107 127 149 151 167 179 181 227 229 233 239 241 263 269 281 307 311 347 349 367 373 401 409 419 431 433 439 461 487 491 503 569 571 587 593 599 601 607 641 643 647 653 659 677 719 727 739 751 769 809 821 823 827 853 857 881 937 941 947 967 983 1,009 1,019 1,021 1,031 1,049 1,051 1,061 1,063 1,087 1,091 1,097 1,103 1,151 1,163 1,187 1,217 1,229 1,249 1,277 1,289 1,297 1,301 1,367 1,373 1,423 1,427 1,429 1,439 The 1,000th Ramanujan prime is 19,403. The 10,000th Ramanujan prime is 242,057. The 100,000th Ramanujan prime is 2,916,539. Elapsed time: 46.0828 milliseconds
Go
A decent time though not as quick as Phix. <lang go>package main
import (
"fmt" "math" "rcu" "sort" "time"
)
var start = time.Now()
var primes = rcu.Primes(700000) // say
func ramanujan(n int) int {
fn := float64(n)
max := int(4 * fn * math.Log(4*fn) / math.Ln2)
pi := sort.SearchInts(primes[2*n:], max) // binary search from min of (2n)th prime
for {
if pi+1-rcu.PrimeCount(primes[pi]/2) <= n {
return primes[pi]
}
pi--
}
return 0
}
func main() {
fmt.Println("The first 100 Ramanujan primes are:")
rams := make([]int, 100)
for n := 0; n < 100; n++ {
rams[n] = ramanujan(n + 1)
}
for i, r := range rams {
fmt.Printf("%5s ", rcu.Commatize(r))
if (i+1)%10 == 0 {
fmt.Println()
}
}
fmt.Printf("\nThe 1,000th Ramanujan prime is %6s\n", rcu.Commatize(ramanujan(1000)))
fmt.Printf("\nThe 10,000th Ramanujan prime is %7s\n", rcu.Commatize(ramanujan(10000)))
fmt.Println("\nTook", time.Since(start))
}</lang>
- Output:
The first 100 Ramanujan primes are:
2 11 17 29 41 47 59 67 71 97
101 107 127 149 151 167 179 181 227 229
233 239 241 263 269 281 307 311 347 349
367 373 401 409 419 431 433 439 461 487
491 503 569 571 587 593 599 601 607 641
643 647 653 659 677 719 727 739 751 769
809 821 823 827 853 857 881 937 941 947
967 983 1,009 1,019 1,021 1,031 1,049 1,051 1,061 1,063
1,087 1,091 1,097 1,103 1,151 1,163 1,187 1,217 1,229 1,249
1,277 1,289 1,297 1,301 1,367 1,373 1,423 1,427 1,429 1,439
The 1,000th Ramanujan prime is 19,403
The 10,000th Ramanujan prime is 242,057
Took 946.193311ms
Julia
<lang julia>using Primes
@time let
MASK = primesmask(625000)
PIVEC = accumulate(+, MASK)
PI(n) = n < 1 ? 0 : PIVEC[n]
function Ramanujan_prime(n)
maxposs = Int(ceil(4n * (log(4n) / log(2))))
for i in maxposs:-1:1
PI(i) - PI(i ÷ 2) < n && return i + 1
end
return 0
end
for i in 1:100
print(lpad(Ramanujan_prime(i), 5), i % 20 == 0 ? "\n" : "")
end
println("\nThe 1000th Ramanujan prime is ", Ramanujan_prime(1000))
println("\nThe 10,000th Ramanujan prime is ", Ramanujan_prime(10000))
end
</lang>
- Output:
2 11 17 29 41 47 59 67 71 97 101 107 127 149 151 167 179 181 227 229 233 239 241 263 269 281 307 311 347 349 367 373 401 409 419 431 433 439 461 487 491 503 569 571 587 593 599 601 607 641 643 647 653 659 677 719 727 739 751 769 809 821 823 827 853 857 881 937 941 947 967 983 1009 1019 1021 1031 1049 1051 1061 1063 1087 1091 1097 1103 1151 1163 1187 1217 1229 1249 1277 1289 1297 1301 1367 1373 1423 1427 1429 1439 The 1000th Ramanujan prime is 19403 The 10,000th Ramanujan prime is 242057 0.272471 seconds (625.44 k allocations: 38.734 MiB, 33.07% compilation time)
Nim
This is a straight translation of Phix version, but I had to add some code to manage prime numbers. Note also that in Nim sequences starts at index 0, not 1.
I compiled using command nim c -d:release -d:lto ramanujan_primes.nim, i.e. with runtime checks on and link time optimization. The program runs in about 35 ms on my laptop (i5-8250U CPU @ 1.60GHz, 8 GB Ram, Linux Manjaro). Fast, but this is normal for native code.
<lang Nim>import algorithm, math, strutils, times
let t0 = now()
const N = 400_000
var composite: array[2..N, bool] for n in 2..N:
let n2 = n * n
if n2 > N: break
if not composite[n]:
for k in countup(n2, N, n):
composite[k] = true
proc primesLe(n: int): seq[int] =
for i, comp in composite: if i > n: break if not comp: result.add i
var piCache: seq[int]
proc pi(n: int): int =
if n == 0: return 0
if n > piCache.len:
let primes = primesLe(n)
for i in piCache.len+1..n:
let k = primes.upperBound(i)
piCache.add k
result = piCache[n-1]
proc ramanujanPrime(n: int): int =
let maxPoss = int(ceil(4 * n.toFloat * ln(4 * n.toFloat)))
for i in countdown(maxPoss, 1):
if pi(i) - pi(i div 2) < n:
return i + 1
for n in 1..100:
stdout.write ($ramanujanPrime(n)).align(4), if n mod 20 == 0: '\n' else: ' '
echo "\nThe 1000th Ramanujan prime is ", ramanujanPrime(1000) echo "The 10000th Ramanujan prime is ", ramanujanPrime(10000)
echo "\nElapsed time: ", (now() - t0).inMilliseconds, " ms"</lang>
- Output:
2 11 17 29 41 47 59 67 71 97 101 107 127 149 151 167 179 181 227 229 233 239 241 263 269 281 307 311 347 349 367 373 401 409 419 431 433 439 461 487 491 503 569 571 587 593 599 601 607 641 643 647 653 659 677 719 727 739 751 769 809 821 823 827 853 857 881 937 941 947 967 983 1009 1019 1021 1031 1049 1051 1061 1063 1087 1091 1097 1103 1151 1163 1187 1217 1229 1249 1277 1289 1297 1301 1367 1373 1423 1427 1429 1439 The 1000th Ramanujan prime is 19403 The 10000th Ramanujan prime is 242057 Elapsed time: 34 ms
Phix
You can run this online here.
with javascript_semantics 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 sequence r = apply(tagset(100),Ramanujan_prime) printf(1,"%s\n",join_by(apply(true,sprintf,{{"%5d"},r}),1,20,"")) printf(1,"The 1000th Ramanujan prime is %d\n", Ramanujan_prime(1000)) printf(1,"The 10,000th Ramanujan prime is %d\n", Ramanujan_prime(10000)) ?elapsed(time()-t0)
- Output:
2 11 17 29 41 47 59 67 71 97 101 107 127 149 151 167 179 181 227 229 233 239 241 263 269 281 307 311 347 349 367 373 401 409 419 431 433 439 461 487 491 503 569 571 587 593 599 601 607 641 643 647 653 659 677 719 727 739 751 769 809 821 823 827 853 857 881 937 941 947 967 983 1009 1019 1021 1031 1049 1051 1061 1063 1087 1091 1097 1103 1151 1163 1187 1217 1229 1249 1277 1289 1297 1301 1367 1373 1423 1427 1429 1439 The 1000th Ramanujan prime is 19403 The 10,000th Ramanujan prime is 242057 "0.6s"
Raku
Not super fast; around 20 seconds on my system. <lang perl6>use Math::Primesieve; use Lingua::EN::Numbers;
my $primes = Math::Primesieve.new;
my @mem;
sub ramanujan-prime (\n) {
1 + (1..(4×n × (4×n).log / 2.log).floor).first: :end, -> \x {
my \y = x div 2;
((@mem[x] //= $primes.count(x)) - (@mem[y] //= $primes.count(y))) < n
}
}
say 'First 100:'; say (1..100).map( &ramanujan-prime ).batch(10)».&comma».fmt("%6s").join: "\n"; say "\n 1,000th: { comma 1000.&ramanujan-prime }"; say "10,000th: { comma 10000.&ramanujan-prime }";</lang>
- Output:
First 100:
2 11 17 29 41 47 59 67 71 97
101 107 127 149 151 167 179 181 227 229
233 239 241 263 269 281 307 311 347 349
367 373 401 409 419 431 433 439 461 487
491 503 569 571 587 593 599 601 607 641
643 647 653 659 677 719 727 739 751 769
809 821 823 827 853 857 881 937 941 947
967 983 1,009 1,019 1,021 1,031 1,049 1,051 1,061 1,063
1,087 1,091 1,097 1,103 1,151 1,163 1,187 1,217 1,229 1,249
1,277 1,289 1,297 1,301 1,367 1,373 1,423 1,427 1,429 1,439
1,000th: 19,403
10,000th: 242,057
Wren
This takes about 28 seconds on my machine. <lang ecmascript>import "/math" for Int import "/seq" for Lst import "/fmt" for Fmt import "/sort" for Find
var primes = Int.primeSieve(700000) // say
var ramanujan = Fn.new { |n|
var max = (4 * n * (4 * n).log / 2.log).floor
var pi = Find.all(primes[2*n..-1], max)[2].from // binary search from min of (2n)th prime
while (true) {
var delta = pi + 1 - Int.primeCount((primes[pi]/2).floor)
if (delta <= n) return primes[pi]
pi = pi - 1
}
}
System.print("The first 100 Ramanujan primes are:") var rams = (1..100).map { |n| ramanujan.call(n) }.toList for (chunk in Lst.chunks(rams, 10)) Fmt.print("$,5d", chunk)
Fmt.print("\nThe 1,000th Ramanujan prime is $,6d", ramanujan.call(1000))
Fmt.print("\nThe 10,000th Ramanujan prime is $,7d", ramanujan.call(10000))</lang>
- Output:
The first 100 Ramanujan primes are:
2 11 17 29 41 47 59 67 71 97
101 107 127 149 151 167 179 181 227 229
233 239 241 263 269 281 307 311 347 349
367 373 401 409 419 431 433 439 461 487
491 503 569 571 587 593 599 601 607 641
643 647 653 659 677 719 727 739 751 769
809 821 823 827 853 857 881 937 941 947
967 983 1,009 1,019 1,021 1,031 1,049 1,051 1,061 1,063
1,087 1,091 1,097 1,103 1,151 1,163 1,187 1,217 1,229 1,249
1,277 1,289 1,297 1,301 1,367 1,373 1,423 1,427 1,429 1,439
The 1,000th Ramanujan prime is 19,403
The 10,000th Ramanujan prime is 242,057