Counting Fractions¶
Like problem 71, we're looking at Farey sequences. This time we're interested in the cardinality of $F_{1000000}$.
To begin, first note that $F_1 = \{0, 1\}$, so $|F_1| = 2$ (this problem isn't counting 0 and 1 in its totals - we'll handle that at the end). Then consider that for any Farey sequence $F_n$, the next sequence $F_{n+1}$ will contain all the terms from $F_n$, along with all irreducible fractions $\frac{k}{n+1}$, (since any reducible fraction would already be in $F_n$).
How many new fractions does this get us? Well, the fraction only reduces if $k$ and $n+1$ have a common factor - in other words, if $k$ and $n+1$ are coprime, the fraction will not reduce. How many number less than $n+1$ are coprime to $n+1$? The totient function will tell us! So the number of irreducible fractions with denominator $n+1$ is simply $\phi(n+1)$. This gives us $$|F_{n+1}| = |F_n| + \phi(n+1)$$ From this, we can derive a non-recursive formula: $$|F_n| = 1 + \sum_{k=1}^n \phi(k)$$
The sum of totients is denoted $\Phi(n)$. As mentioned before, we'll actually subtract two from this total, since the problem isn't counting 0 or 1. So this problem boils down to computing $\Phi(1000000) - 1$.
Using SageMath¶
Since SageMath implements a decently fast totient function, we could opt for the most straightforward approach.
sum(euler_phi(n) for n in range(1, 1000001)) - 1
303963552391
If you try to implement the totient function yourself, you might find it difficult to make this approach fast enough. An alternative is to "sieve" the values of totient - see problem 70. However, there's a few more very interesting methods to compute totient sums.
Recursive totient sum¶
You might think a fast totient function is key to solving this problem quickly, but it turns out there's a crazy recursive formula for $\Phi(n)$ that does not require implementing the totient function at all! Since it's a little hard to believe it works, here's a derivation of it.
We start with the following identity, proven by Gauss: $$\sum_{d|k}\phi(d) = k$$
Then, we can apply the formula for triangular numbers: $$\sum_{k=1}^n \sum_{d|k}\phi(d) = \sum_{k=1}^n k = \frac{n(n+1)}{2}$$
Now we'll rewrite the bounds. $$\sum_{1 \leq k \leq n} \sum_{xy=k}\phi(x) = \frac{n(n+1)}{2}$$
This makes it easier to see that we can condense the two nested summations into one. $$\sum_{1 \leq xy \leq n} \phi(x) = \frac{n(n+1)}{2}$$
We can then split this single summation into a sum-of-summations, letting $y$ vary from 1 to $n$. $$\sum_{1 \leq x \leq n} \phi(x) + \sum_{1 \leq 2x \leq n} \phi(x) + \sum_{1 \leq 3x \leq n} \phi(x) + \cdots + \sum_{1 \leq xn \leq n} \phi(x) = \frac{n(n+1)}{2}$$
If we apply the definition of $\Phi(n)$, we get $$\Phi\left(\left\lfloor\frac{n}{1}\right\rfloor\right) + \Phi\left(\left\lfloor\frac{n}{2}\right\rfloor\right) + \Phi\left(\left\lfloor\frac{n}{3}\right\rfloor\right) + \cdots + \Phi\left(\left\lfloor\frac{n}{n}\right\rfloor\right) = \frac{n(n+1)}{2}$$
Then we can rearrange and simplify to get $$\Phi(n) = \frac{1}{2}n(n+1) - \sum_{d=2}^n \Phi\left(\left\lfloor \frac{n}{d} \right\rfloor\right)$$
Not a super intuitive definition for $\Phi(n)$, but pretty simple to implement, and surprisingly quick!
from functools import cache
@cache
def totient_sum(n):
if n == 1:
return 1
return n*(n+1)/2 - sum(totient_sum(n//d) for d in range(2, n + 1))
totient_sum(1000000) - 1
303963552391
Both of the above strategies are simple and fast enough to get the job done on a modern computer, but it turns out we can go even deeper. But first, we need to become familiar with some important formulas.
Dirichlet convolution¶
The Dirichlet convolution of two arithmetic functions $f$ and $g$ is $$(f * g)(n) = \sum_{xy=n} f(x) g(y)$$
Several important mathematical functions have some representation as a Dirichlet convolution - see the above page for a list. An important one for our purposes is $\phi = \mathrm{Id} * \mu$, where $\mathrm{Id(n) = n}$ is the identity function and $\mu$ is the Möbius function - if you've never heard of the latter, don't worry about it yet.
Dirichlet's hyperbola method¶
By applying the above Dirichlet convolution to $$\Phi(n) = \sum_{k=1}^n \phi(n)$$ we can derive an equivalent summation: $$\Phi(n) = \sum_{k=1}^n \sum_{xy=k} x \mu(y)$$
We're going to transform this summation by applying Dirichlet's hyperbola method. Let $1 < a < n$, and let $b = n / a$. Then $$\Phi(n) = \sum_{x=1}^a \sum_{y=1}^{n/x} x \mu(y) + \sum_{y=1}^b \sum_{x=1}^{n/y} x \mu(y) - \sum_{x=1}^a \sum_{y=1}^b x \mu(y)$$
Then with some algebra, we can transform this into $$\Phi(n) = \sum_{x=1}^a x M\left(\left\lfloor \frac{n}{x}\right\rfloor\right) + \sum_{y=1}^b \mu(y) \frac{\left\lfloor \frac{n}{y} \right\rfloor \left(\left\lfloor \frac{n}{y}\right\rfloor + 1\right)}{2} - M(\lfloor b \rfloor) \frac{\lfloor a \rfloor(\lfloor a\rfloor +1)}{2}$$ where $M(n)$ is the Mertens function, which is just the partial sums of the Möbius function: $$M(n) = \sum_{k=1}^n \mu(k)$$
By using this formula for $\Phi(n)$, instead of calculating $n$ different totients, our sums only run to $a$ and $b$.
However, for this formula to be effective, we also need an efficient way to calculate $M(n)$. Fortunately, we can just apply the hyperbola method again, with the convolution $\epsilon = \mu * 1$, where $\epsilon$ is the unit identity. Once again, we choose $1 < a < n$ and let $b = n/a$.
$$\sum_{k=1}^n \epsilon(k) = \sum_{k=1}^n \sum_{xy=k} \mu(x)$$$$1 = \sum_{x=1}^a \sum_{y=1}^{n/x} \mu(x) + \sum_{y=1}^b \sum_{x=1}^{n/y} \mu(x) - \sum_{x=1}^a \sum_{y=1}^b \mu(x)$$$$1 = \sum_{x=1}^a \mu(x)\left\lfloor \frac{n}{x}\right\rfloor + \sum_{y=1}^b M\left(\left\lfloor \frac{n}{y}\right\rfloor\right) - \lfloor b \rfloor M(\lfloor a\rfloor)$$$$M(n) = 1 + \lfloor b\rfloor M(\lfloor a\rfloor) - \sum_{x=1}^a \mu(x)\left\lfloor \frac{n}{x}\right\rfloor - \sum_{y=2}^b M\left(\left\lfloor \frac{n}{y}\right\rfloor\right)$$Now we have a recursive implementation of $M(n)$. All that's left is to calculate the values of $\mu(n)$ that we need. By taking advantage of $\mu$ being multiplicative, we can compute values using a similar strategy to the sieve of Eratosthenes - see problem 10.
def mobius_range(limit):
ms = [n for n in range(0, limit)]
for n in range(0, limit):
if n == 0 or n == 1 or ms[n] != n:
yield ms[n]
continue
yield -1
for k in range(2 * n, limit, n):
if k % n ^ 2 == 0:
ms[k] = 0
ms[k] //= -n
We'll use $a = \sqrt{1000000} = 1000$ as our upper bound.
from math import isqrt
mu = list(mobius_range(isqrt(1000000) + 1))
Now we can implement $M(n)$ with our recursive formula:
@cache
def M(n):
if n == 0 or n == 1:
return n
a = sqrt(n)
b = n / a
total = 1 + floor(b) * M(floor(a))
total -= sum(mu[x] * floor(n/x) for x in range(1, floor(a) + 1))
total -= sum(M(floor(n/y)) for y in range(2, floor(b) + 1))
return total
And finally, yet another implementation of $\Phi(n)$:
def totient_sum(n):
a = sqrt(n)
b = n / a
total = sum(x * M(floor(n/x)) for x in range(1, floor(a) + 1))
total += sum(polygonal_number(3, floor(n/y)) * mu[y] for y in range(1, floor(b) + 1))
total -= M(floor(b)) * polygonal_number(3, floor(a))
return total
This approach is significantly faster than the others.
totient_sum(1000000) - 1
303963552391
Relevant sequences¶
Copyright (C) 2025 filifa¶
This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International license and the BSD Zero Clause license.