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.





with javascript_semantics
--
-- While a phix dictionary can handle keys of {x,a}, for this 
-- task performance was dreadful (7,612,479 entries, >3mins),
-- so instead memophix maps known x to an index to memophia 
-- which holds the full [1..a] for each x, dropping to a much
-- more respectable (albeit not super-fast) 14s
--
constant memophix = new_dict()
sequence memophia = {}  -- 1..a (max 3401) for each x

function phi(integer x, a)
    if a=0 then return x end if
    integer adx = getd(x,memophix), res
    if adx=NULL then
        memophia = append(memophia,repeat(-1,a))
        adx = length(memophia)
        setd(x,adx,memophix)
    else
        object ma = memophia[adx]
        integer l = length(ma)
        if a>l then
            memophia[adx] = 0   -- kill refcount
            memophia[adx] = ma & repeat(-1,a-l)
        else
            res = ma[a]
            if res>=0 then return res end if
        end if
        ma = 0                  -- kill refcount
    end if
    res = phi(x, a-1) - phi(floor(x/get_prime(a)), a-1)
    memophia[adx][a] = res
    return res
end function
 
function pi(integer n)
    if n<2 then return 0 end if
    integer a = pi(floor(sqrt(n)))
    return phi(n, a) + a - 1
end function
 
atom t0 = time()
for i=0 to iff(platform()=JS?8:9) do
    printf(1,"10^%d  %d\n",{i,pi(power(10,i))})
--  printf(1,"10^%d  %d\n",{i,length(get_primes_le(power(10,i)))})
end for
?elapsed(time()-t0)
