Legendre prime counting function
You are encouraged to solve this task according to the task description, using any language you may know.
The prime-counting function π(n) computes the number of primes not greater than n. Legendre was the first mathematician to create a formula to compute π(n) based on the inclusion/exclusion principle.
To calculate:
Define
φ(x, 0) = x φ(x, a) = φ(x, a−1) − φ(⌊x/pa⌋, a−1), where pa is the ath prime number.
then
π(n) = 0 when n < 2 π(n) = φ(n, a) + a - 1, where a = π(√n), n ≥ 2
The Legendre formula still requires the use of a sieve to enumerate primes; however it's only required to sieve up to the √n, and for counting primes, the Legendre method is generally much faster than sieving up to n.
- Task
Calculate π(n) for values up to 1 billion. Show π(n) for n = 1, 10, 100, ... 109.
For this task, you may refer to a prime number sieve (such as the Sieve of Eratosthenes or the extensible sieve) in an external library to enumerate the primes required by the formula. Also note that it will be necessary to memoize the results of φ(x, a) in order to have reasonable performance, since the recurrence relation would otherwise take exponential time.
CoffeeScript
<lang CoffeeScript> sorenson = require('sieve').primes # Sorenson's extensible sieve from task: Extensible Prime Generator
- put in outer scope to avoid recomputing the cache
memoPhi = {} primes = []
isqrt = (x) -> Math.floor Math.sqrt x
pi = (n) ->
phi = (x, a) ->
y = memoPhix,a
return y unless y is undefined
memoPhix,a = if a is 0 then x else p = primes[a - 1] throw "You need to generate at least #{a} primes." if p is undefined phi(x, a - 1) - phi(x // p, a - 1)
if n < 2
0
else
a = pi isqrt n
phi(n, a) + a - 1
maxPi = 1e9 gen = sorenson() primes = while (p = gen.next().value) < isqrt maxPi then p
n = 1 for i in [0..9]
console.log "10^#{i}\t#{pi(n)}"
n *= 10
</lang>
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
Erlang
<lang Erlang> -mode(native).
-define(LIMIT, 1000000000).
main(_) ->
put(primes, array:from_list(primality:sieve(floor(math:sqrt(?LIMIT))))), ets:new(memphi, [set, named_table, protected]), output(0, 9).
nthprime(N) -> array:get(N - 1, get(primes)).
output(A, B) -> output(A, B, 1). output(A, B, _) when A > B -> ok; output(A, B, N) ->
io:format("10^~b ~b~n", [A, pi(N)]),
output(A + 1, B, N * 10).
pi(N) ->
Primes = get(primes),
Last = array:get(array:size(Primes) - 1, Primes),
if
N =< Last -> small_pi(N);
true ->
A = pi(floor(math:sqrt(N))),
phi(N, A) + A - 1
end.
phi(X, 0) -> X; phi(X, A) ->
case ets:lookup(memphi, {X, A}) of
[] ->
Phi = phi(X, A-1) - phi(X div nthprime(A), A-1),
ets:insert(memphi, {{X, A}, Phi}),
Phi;
[{{X, A}, Phi}] -> Phi
end.
% Use binary search to count primes that we have already listed. small_pi(N) ->
Primes = get(primes), small_pi(N, Primes, 0, array:size(Primes)).
small_pi(_, _, L, H) when L >= (H - 1) -> L + 1; small_pi(N, Primes, L, H) ->
M = (L + H) div 2,
P = array:get(M, Primes),
if
N > P -> small_pi(N, Primes, M, H);
N < P -> small_pi(N, Primes, 0, M);
true -> M + 1
end.
</lang>
- Output:
10^0 1 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
Go
This runs in about 2.2 seconds which is surprisingly slow for Go. However, Go has the same problem as Wren in not being able to use lists for map keys. I tried using a 2-element array for the key but this was a second slower than the Cantor function.
The alternative of sieving a billion numbers and then counting them takes about 4.8 seconds using the present library function though this could be cut to 1.2 seconds using a third party library such as primegen. <lang go>package main
import (
"fmt" "log" "math" "rcu"
)
var limit = int(math.Sqrt(1e9)) var primes = rcu.Primes(limit) var memoPhi = make(map[int]int)
func cantorPair(x, y int) int {
if x < 0 || y < 0 {
log.Fatal("Arguments must be non-negative integers.")
}
return (x*x + 3*x + 2*x*y + y + y*y) / 2
}
func phi(x, a int) int {
if a == 0 {
return x
}
key := cantorPair(x, a)
if v, ok := memoPhi[key]; ok {
return v
}
pa := primes[a-1]
memoPhi[key] = phi(x, a-1) - phi(x/pa, a-1)
return memoPhi[key]
}
func pi(n int) int {
if n < 2 {
return 0
}
a := pi(int(math.Sqrt(float64(n))))
return phi(n, a) + a - 1
}
func main() {
for i, n := 0, 1; i <= 9; i, n = i+1, n*10 {
fmt.Printf("10^%d %d\n", i, pi(n))
}
}</lang>
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
Julia
Using the Julia Primes library
This would be faster with a custom sieve that only counted instead of making a bitset, instead of the more versatile library function. <lang julia>using Primes
function primepi(N)
delta = round(Int, N^0.8) return sum(i -> count(primesmask(i, min(i + delta - 1, N))), 1:delta:N)
end
@time for power in 0:9
println("10^", rpad(power, 5), primepi(10^power))
end
</lang>
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534 1.935556 seconds (1.64 k allocations: 415.526 MiB, 2.15% gc time)
Using the method given in the task
Run twice to show the benefits of precomputing the library prime functions. <lang julia>using Primes
const maxpi = 1_000_000_000 const memophi = Dict{Vector{Int}, Int}() const _primes = primes(isqrt(maxpi))
function primepi(n)
function phi(x, a)
if !haskey(memophi, [x, a])
memophix, a = a == 0 ? x : phi(x, a - 1) - phi(x ÷ _primes[a], a - 1)
end
return memophix, a
end
n < 2 && return 0 a = primepi(isqrt(n)) return phi(n, a) + a - 1
end
@time for i in 0:9
println("10^", rpad(i, 5), primepi(10^i))
end
@time for i in 0:9
println("10^", rpad(i, 5), primepi(10^i))
end
</lang>
- Output:
10^0 1 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534 12.939211 seconds (69.44 M allocations: 3.860 GiB, 27.98% gc time, 1.58% compilation time) 10^0 1 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534 0.003234 seconds (421 allocations: 14.547 KiB)
Mathematica/Wolfram Language
<lang Mathematica>ClearAll[Phi,pi] $RecursionLimit = 10^6; Phi[x_, 0] := x Phi[x_, a_] := Phi[x, a] = Phi[x, a - 1] - Phi[Floor[x/Prime[a]], a - 1] pi[n_] := Module[{a}, If[n < 2, 0, a = pi[Floor[Sqrt[n]]]; Phi[n, a] + a - 1]] Scan[Print[pi[10^#]] &, Range[0,9]]</lang>
- Output:
0 4 25 168 1229 9592 78498 664579 5761455 50847534
Nim
The program runs in 2.2 seconds on our laptop. Using int32 instead of naturals (64 bits on our 64 bits computer) saves 0.3 second. Using an integer rather than a tuple as key (for instance x shl 32 or a) saves 0.4 second.
<lang Nim>import math, strutils, sugar, tables
const
N = 1_000_000_000 S = sqrt(N.toFloat).int
var composite: array[3..S, bool] for n in countup(3, S, 2):
if n * n > S: break
if not composite[n]:
for k in countup(n * n, S, 2 * n):
composite[k] = true
- Prime list. Add a dummy zero to start at index 1.
let primes = @[0, 2] & collect(newSeq, for n in countup(3, S, 2): (if not composite[n]: n))
var cache: Table[(Natural, Natural), Natural]
proc phi(x, a: Natural): Natural =
if a == 0: return x let pair = (x, a) if pair in cache: return cache[pair] result = phi(x, a - 1) - phi(x div primes[a], a - 1) cache[pair] = result
proc π(n: Natural): Natural =
if n <= 2: return 0 let a = π(sqrt(n.toFloat).Natural) result = phi(n, a) + a - 1
var n = 1 for i in 0..9:
echo "π(10^$1) = $2".format(i, π(n)) n *= 10</lang>
- Output:
π(10^0) = 0 π(10^1) = 4 π(10^2) = 25 π(10^3) = 168 π(10^4) = 1229 π(10^5) = 9592 π(10^6) = 78498 π(10^7) = 664579 π(10^8) = 5761455 π(10^9) = 50847534
Picat
Performance is surprisingly good, considering that the Picat's library SoE is very basic and that I recompute the base primes each time pi() is called. This version takes 23 seconds on my laptop, outperforming the CoffeeScript version which takes 30 seconds. <lang Picat> table(+, +, nt) phi(X, 0, _) = X. phi(X, A, Primes) = phi(X, A - 1, Primes) - phi(X // Primes[A], A - 1, Primes).
pi(N) = Count, N < 2 => Count = 0. pi(N) = Count =>
M = floor(sqrt(N)), A = pi(M), Count = phi(N, A, math.primes(M)) + A - 1.
main =>
N = 1,
foreach (K in 0..9)
writef("10^%w %w%n", K, pi(N)),
N := N * 10.
</lang>
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
Raku
Seems like an awful lot of pointless work. Using prime sieve anyway, why not just use it? <lang perl6>use Math::Primesieve;
my $sieve = Math::Primesieve.new;
say "10^$_\t" ~ $sieve.count: exp($_,10) for ^10;
say (now - INIT now) ~ ' elapsed seconds';</lang>
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534 0.071464489 elapsed seconds
Wren
This runs in about 5.7 seconds which is not too bad for the Wren interpreter. As map keys cannot be lists, the Cantor pairing function has been used to represent [x, a] which is considerably faster than using a string based approach for memoization.
To sieve a billion numbers and then count the primes up to 10^k would take the best part of 3 minutes in Wren so, as expected, the Legendre method represents a considerable speed up. <lang ecmascript>import "/math" for Int
var limit = 1e9.sqrt.floor var primes = Int.primeSieve(limit) var memoPhi = {}
var phi // recursive function phi = Fn.new { |x, a|
if (a == 0) return x var key = Int.cantorPair(x, a) if (memoPhi.containsKey(key)) return memoPhi[key] var pa = primes[a-1] return memoPhi[key] = phi.call(x, a-1) - phi.call((x/pa).floor, a-1)
}
var pi // recursive function pi = Fn.new { |n|
if (n < 2) return 0 var a = pi.call(n.sqrt.floor) return phi.call(n, a) + a - 1
}
var n = 1 for (i in 0..9) {
System.print("10^%(i) %(pi.call(n))")
n = n * 10
}</lang>
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534