Computing π(n)\pi(n)

Wee Kean here. Happy New Years Eve! Lately, my only interaction with math has been through doing Project Euler. So... here's some cool things I've learnt! Unfortunately, there is hardly any structure to this blogpost so I apologize if this comes off as verbal diarrhea.

$$\newcommand\floor[1]{\left \lfloor #1 \right \rfloor}$$

Orders of Growth

Throughout this post, I will write things like "this algorithm runs in $O(n)$ time" or "this requires $O(\sqrt{n})$ space". If you're not familiar with Big O Notation, the term $O(g(x))$ represents a class of functions $f$ such that for all sufficiently large values of $x$, $|f(x)| \leq Mg(x)$. More simply, you can think of $O(g(x))$ as a loose upper bound, ignoring constants. So when I say "this algorithm runs in $O(n)$ time", what I mean is that the time taken to run this algorithm is at most linear in $n$ for large enough $n$.

Everything Prime

Let's start with a simple task - generating all primes up to $n$. I'm sure most have heard of the Sieve of Eratosthenes. We start with an array of numbers from 2 to $n$, which denotes the numbers which may be prime. Iterating from 2 to $n$, if the current number is marked as composite, we do nothing else this number is prime and we mark all multiples of it as composite.

Let's estimate the run time of this algorithm. At each step, if the current number is composite, we perform 0 operations. If the number is prime, we mark $\frac{n}{p}$ numbers in our array as composite. The total number of operations is roughly $\sum_{p \leq n} \frac{n}{p} \approx n \log \log n$.

For all purposes, this is good enough since $\log \log n$ is small. But recently, I stumbled across the following post which details a linear sieve! The idea is that certain composite numbers get set as composite multiple times and we can choose a unique representation for each composite instead.

Here are some runtimes of the above two algorithms implemented in pypy running on my computer. (This is hardly scientific, it's just for fun.)

We can actually optimize further, cutting down constant factors by considering only odd numbers but this is pretty boring. Let's look at something abit more challenging. The prime-counting function $\pi (n)$ denotes the number of primes less than or equal to some positive integer $n$. We already have an algorithm to generate all primes less than $n$ so finding $\pi (n)$ would take at most $O(n)$ time. But looking at our runtimes, it will take quite a long while (possibly not in my lifetime) before we can find $\pi (n)$ for larger values of $n = 10^{10}, 10^{11}, 10^{12}$, not to mention the limitations in memory.

Okay, let's look back at our original sieve. If we are only concerned with the number of primes in the range [2, $n$], then really we only care about how many numbers get marked as composite. Also if you've ever done the sieve by hand, you should notice that after sieving out the primes less than $\sqrt{n}$, the remaining numbers are all prime. This is just because every composite number has a prime factor less than $\sqrt{n}$. So to find $\pi (n)$, we really only need to know all primes up to $\sqrt{n}$, then we can use Principle of Inclusion and Exclusion.

Let's look at an example. Say for $n = 10$, the primes up to $\sqrt{n}$ are $p=2,3$. Then, there are $\floor{\frac{10}{2}} + \floor{ \frac{10}{3} } - \floor {\frac{10}{6} } = 7$ numbers marked as composite in the range [2,10]. So there are 2 prime numbers (not including 2, 3) left in the range.

But hold on, doing PIE this way means we have to calculate $2^k$ fractions where $k = \pi (\sqrt{n})$ which is even worse than our linear algorithm before...

$\pi e$ and $\sqrt{\cdot}$ trick

Let's take a slight detour. Let $d(n)$ be the number of divisors of $n$. Can we compute $D(n) = \sum_{x \leq n} d(x)$ efficiently? I think many would have seen the double counting argument that shows $D(n) = \sum \floor { \frac{n}{x} }$. This approach runs in $O(n)$ time. But on closer inspection, there is a large portion of this summation which is awfully inefficient. For example, $\floor {\frac{n}{x} } = 1$ for $x \in (\frac{n}{2}, n]$ where our program will simply be summing 1s over and over. In particular, for all $y < \sqrt{n}$, the number of $x$ with $\floor {\frac{n}{x}} = y$ is just the number of integers in the range $(\frac{n}{y+1}, \frac{n}{y}]$. So, we can split $$D(n) = \sum_{x \leq \sqrt{n}} \floor{\frac{x}{n}} + \sum_{x > \sqrt{n}}  \floor{\frac{x}{n}}$$ which we can calculate both parts in $O(\sqrt{n})$ time.

Right, so where does this leave us? Well, our PIE in the previous section has a similar inefficiency where for large $n$ there are similarly many fractions which evaluate to 0. Our main trouble now is that this form of PIE is very troublesome and it is hard to see how we can optimize the calculation of these "0s".

Instead, let's look at PIE with two sets. Suppose the primes less than $\sqrt{n}$ are $p_1, p_2, ..., p_m$. Let $A$ be the set of numbers divisible by $p_m$, $B$ be the set of numbers divisible by at least one of $p_1, p_2, ..., p_{m-1}$. We want to find $$|A \cup B| = |A| + |B| - |A \cap B|$$ We know $|A| = \floor {\frac{n}{p_m}}$, $|A \cap B|$ is the number of positive integers $x = p_m y \iff y = \frac{x}{p_m} \leq \floor {\frac{n}{p_m}}$. Do you see a pattern? If $S(n, k)$ denotes the number of integers divisble by the first $k$ primes in the range [2, n], then $S(n, k) = \floor {\frac{n}{p_k}} + S(n, k-1) - S\left(\floor {\frac{n}{p_k}}, k-1\right)$. More succinctly, if $\phi (x, a)$ denotes the number of positive integers less than or equal to $x$ which are not divisible by any of the first $a$ primes then $$\phi (x,a) = \phi (x, a-1) - \phi\left(\floor { \frac{x}{p_a} }, a-1\right)$$ which is Legendre's Formula. Let's look at the runtime of this algorithm without any further optimization. We need to calculate a large table of values of the function $\phi (x, a)$ but since $\floor {\frac{\floor {\frac{n}{a}}}{b}} = \floor {\frac{n}{ab}}$, we have that $x$ takes at most $2\sqrt{n}$ values. We need to calculate for all primes less than $\sqrt{n}$ and there are around $\frac{\sqrt{n}}{\log{\sqrt{n}}}$ of those. Hence our time complexity is at most $O(\frac{n}{\log n})$. Hmm, a great reduction from exponential but not exactly a great improvement from our linear sieve.

We can bring down the algorithm to $O(n^{\frac{3}{4}})$ by omitting the calculation of certain values of $\phi(x, a)$. Recall that we only need to sieve up till $\sqrt{x}$. For any particular value of $x$, if $p_a ^ 2 > x$, then $\phi(x, a) = \phi(x, a-1)$. To analyse our time complexity, for each prime $p_a$, we calculate the number of values $x$ for which we need to calculate $\phi(x, a)$. For primes $p_a < n^{1/4}$, as a loose upper bound, there are at most $2\sqrt{x}$ values which need to be updated. So in total, $O(\frac{n^{\frac{1}{4}}}{\log n} \cdot n^{\frac{1}{2}})$ operations. For primes $p > n^{\frac{1}{4}}$, we only calculate for values of $\floor {\frac{n}{a}} = x > p^2 \iff a < \frac{n}{p^2}$. So at most $\sum_{p>n^{\frac{1}{4}}} \frac{n}{p^2} = O(n^{\frac{3}{4}})$ operations.

More runtimes in pypy:

Multiplicative Functions

One last side note. Recall our divisor function - $d(n)$. This is a member of a larger class of functions called multiplicative functions where $f(mn) = f(m)f(n)$ for any $m, n$ such that $\gcd(m, n) = 1$. This class of functions has very nice and interesting properties. In general, we can apply the same sieving techniques shown above to compute values of any arithmetic function $f(x)$ for all $x \leq n$ in $O(n)$ time.

For two multiplicative functions $f, g$, their dirichlet convolution is defined to be $$(f*g)(n) = \sum_{ab=n} f(a)g(b)$$

For example, $d = 1 * 1$. Then, there is a neat little trick to evaluate $\sum_{x \leq n} (f*g)(n)$ if you can evaluate the partial sums of $f, g$ quickly. In particular, suppose $G(n) = \sum_{x \leq n} g(x), F(n) = \sum_{x \leq n} f(x)$.

$$\begin{aligned} \sum_{x \leq n} (f*g)(x) &= \sum_{x \leq n} \sum_{ab = x} f(a)g(b) \\&= \sum_{a \leq \sqrt{n}} \sum_{b \leq \frac{n}{a}} f(a)g(b) + \sum_{b \leq \sqrt{n}} \sum_{a \leq \frac{n}{b}} f(a)g(b) - \sum_{a \leq \sqrt{n}} \sum_{b \leq \sqrt{n}} f(a)g(b) \\&= \sum_{a \leq \sqrt{n}} f(a)G\left(\frac{n}{a}\right) + \sum_{b \leq \sqrt{n}} F(\frac{n}{a}) g(b) - F(\sqrt{n}) G(\sqrt{n}) \end{aligned}$$

Then, if you can evaluate the partial sums $F, G$ in $O(1)$ time, you can evaluate the above in $O(\sqrt{n})$ time. This is called the Dirichlet Hyperbola Method where the name presumably comes from the fact that we're counting over lattice points under the hyberbola $ab=n$.

Unfortunately, I'm not too familiar with the bigger picture of dirichlet convolutions in number theory. That's all I have! Enjoy the new year!