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) = φ(x, a) + a - 1, where a = π(√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 Number 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).
nthprime(N) -> array:get(N - 1, get(primes)).
pi(N) ->
Primes = get(primes), Last = array:get(array:size(Primes) - 1, Primes), if N =< Last -> small_pi(N); true -> A = pi(imath:sqrt(N)), phi(N, A) + A - 1 end.
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.
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.
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).
main(_) ->
put(primes, array:from_list(primality:sieve(imath:sqrt(?LIMIT)))), ets:new(memphi, [set, named_table, protected]), output(0, 9).
</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
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