javascript libraries: mathjax, prism.js

For Project Euler Problem 70 I needed to compute a lot of Euler Totients quickly. Solutions are supposed to be under a minute and I was nowhere near that on the first attempt.

The Euler Totient function \(\phi(n) \) is common in mathematics, particularly in number theory and group theory. For a given integer \( n \), the totient function computes the number of relatively prime positive integers less than \(n\) with no common factors with \(n\). Given a prime factorization, the totient can be computed with the formula
\[\phi(n) = (p_1 - 1)p_1^{k_1 - 1} \cdots (p_r - 1)p_r^{k_r - 1}\]
where \(n = p_1^{k_1} \cdots p_r^{k_r}\). Computing prime factorizations is *slow*, so we will exploit the following formula, which only requires the set of prime divisors: \[ \phi(n) = n \cdot \prod_{p|n}{\left(1 - \frac{1}{p}\right)},\]
where the sum ranges over the prime divisors of \(n\). Unfortunately would still be fairly slow because computing the set of prime divisors directly also takes a lot of time, but there is a trick to avoid computing the set of prime divisors of each \(n\) individually.

To solve Problem 70, we need to compute the totients of all the numbers between \(1\) and \(10^7\). We could try to exploit the fact that \(\phi(nm) = \phi(n)\phi(m) \) if \(n\) and \(m\) are relatively prime and memoize, but that is still fairly slow. Instead, we will compute all the totients at once! The fact that you only need the set of prime divisors and not the full factorization is what makes this possible and quick. Let's write an algorithm in Python to compute the totients for all the numbers between \(1\) and \(n\).

We need the list of primes up to \(n\), so assume you have a function `primes_up_to(n)`

that computes the list of primes, such as the following:

```
def primes_up_to(n):
# http://code.activestate.com/recipes/576640/
nroot = int(math.sqrt(n))
sieve = [True] * (n + 1)# range(n + 1)
sieve[0] = False
sieve[1] = False
for i in range(2, nroot+1):
if sieve[i]:
m = n / i - i
sieve[i * i: n + 1: i] = [False] * (int(m) + 1)
return [i for i in range(n+1) if sieve[i]]
```

This is *fast* -- it takes less than one second to generate all the primes up to \(10^7\) on my machine. We also define a convenience generator that generates the list of multiples of a given number:

```
def multiple_gen(modulus):
"""Generates the list of multiples of modulus."""
count = 1
while True:
yield modulus * count
count += 1
```

Now here is the trick to compute the totients iterating through the list of primes only once, rather than many times to compute the prime factorization for every value less than \(n\). First rewrite the formula \[ \phi(n) = n \cdot \prod_{p|n}{\left(1 - \frac{1}{p}\right)} = n \cdot \prod_{p|n}{\left(\frac{p-1}{p}\right)}. \] Start with a list of the numbers from \(1\) to \(n\). For each prime \(p\) and every entry with index divisible by \(p\), multiply by \(p-1\) and divide by \(p\). This will be an integer because \(n\) will always be divisible by \(p\) for all prime divisors. By using the multiple generator, we can avoid testing for divisibility as well since we want to do this for every divisor.

The code:

```
def fast_compute_totients(bound):
"""Simultaneously compute totients for all numbers
up to and including n."""
# Populate the initial list with the leading factor of n.
results = list(range(1, bound+1))
# Get the list of primes up to the bound.
primes = primes_up_to(bound)
for p in primes:
for m in multiple_gen(p):
if m > bound:
break
results[m-1] = (results[m-1] / p) * (p - 1)
return results
```

It takes about 7 seconds to generate all the totients on my machine. Note that if we had used a formula depending on the full prime factorization, this would not have been so easy because we would need to multiply by multiples of prime powers, making sure to avoid multiples of lower prime powers! Also notice that in the inner loop, once the greatest prime divisor for a particular number has been passed, the corresponding entry of results holds the totient. (So the algorithm does not *quite* compute them all simultaneously since some are completed before others.)

What is the complexity of this algorithm? By the prime number theorem, the number of primes up to \(n\) is asymptotic to \(\frac{n}{\log{n}}\). The number of iterations of the inner loop is bounded above by \[n \sum_{p < n}{ \frac{1}{p}}\] since there are approximately \(n / p\) multiples of \(p\) less than \(n\), for each \(p\). This sum is asymptotic to \(n * \log(\log(n))\) (see here), so the fast_compute_totients function iterates at most \[ \frac{n}{\log{n}} * n \log{\log{n}} = \frac{n^2 \log{\log{n}}}{\log{n}} < n^2,\] not too shabby! To compute the totients individually would require roughly (this is a very loose approximation) \[ n * \sum_{k=1}^{n}{\frac{\sqrt{k}}{\log{\sqrt{k}}}} > \frac{2 n^{5/2}}{\log{n}}\] iterations, with each iteration having additional contributions from the now necessary divisibility tests. We will have to be satisfied with the iteration count comparison -- to calculate the actual computational complexity requires knowing exactly how python implements multiplication, division, and modular arithmetic.

Here's the full code:

```
import math
def primes_up_to(n):
# http://code.activestate.com/recipes/576640/
nroot = int(math.sqrt(n))
sieve = [True] * (n+1)# range(n+1)
sieve[0] = False
sieve[1] = False
for i in range(2, nroot+1):
if sieve[i]:
m = n/i - i
sieve[i*i: n+1:i] = [False] * (int(m)+1)
return [i for i in range(n+1) if sieve[i]]
def multiple_gen(modulus):
"""Generates the list of multiples of modulus."""
count = 1
while True:
yield modulus * count
count += 1
def fast_compute_totients(bound):
"""Simultaneously compute totients for all numbers
up to and including n."""
# Populate the initial list with the leading factor of n.
results = list(range(1, bound+1))
# Get the list of primes up to the bound.
primes = primes_up_to(bound)
for p in primes:
for m in multiple_gen(p):
if m > bound:
break
results[m-1] = (results[m-1] / p) * (p - 1)
return results
if __name__ == '__main__':
fast_compute_totients(10**7)
```