Distribution of prime powers

The prime number theorem says that I(x), the number of primes less than or equal to x, is asymptotically x / log x. So it's easy to estimate the number of primes below some number N. But what if we want to estimate the number of prime powers less than N? This is a question that comes up in finite fields, for example, since there is a finite field with n elements if and only if n is a prime power. It's also important in finite simple groups because these groups are often indexed by prime powers.

Riemann's prime-power counting function I(x) counts the number of prime powers less than or equal to x. Clearly I(x) > I(x) for x a 4 since every prime is a prime power, and 4 is a prime power but not a prime. Is I(x) much bigger than I(x)? What is its limiting distribution, i.e. what is the analog of the prime number theorem for prime powers?

Numerical examples

It turns out I(x) equals I(x) asymptotically. That is, even though I(x) is always bigger than I(x), their ratio converges to 1 as x increases.

Why is that? Let's first look at N = 1,000,000. The number of primes less than one million happens to be 78,498. The number of prime powers less than N is 78,734. So the latter includes only 236 prime powers with exponent greater than 1.

If we increase N to 1,000,000,000, there are 50,847,534 primes less than N and 50,851,223 prime powers, a difference of 3,689. Said another way, 99.99% of the prime powers less than a billion have exponent 1.

Equation for I(x)

The number of prime powers less than N with exponent 2 equals the number of primes less than the square root of N. And the number of prime powers less than N with exponent 3 equals the number of primes less than the cube root of N. The number of prime powers with exponent 4 equals the number of primes less than the fourth root of N. Etc.

Even if N is large, these counts start getting small pretty soon. How soon? We're taking roots of order r until the rth root of N is less than 2, because then there are no more primes less than that root. That means we keep going until r > log₂N. And so we have the following equation:

Mathematica and Python code

I looked in Mathematica and SymPy for a function to compute I(x) and didn't see one. Maybe I missed something. But in either case it's easy to implement our own using the equation above.

In Mathematica:

pp[n_] := Sum[PrimePi[n^(1/r)], {r, 1, Log2[n]}]

In Python:

from sympy import primepifrom math import log2def pp(n): top = int(log2(n)) return sum( [primepi(n**(1/r)) for r in range(1, 1+top)] )