Legendre prime counting function

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/p_a⌋, a−1), where p_a is the a^th 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, ... 10⁹.

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