"Like [problem 71](https://projecteuler.net/problem=71), we're looking at [Farey sequences](https://en.wikipedia.org/wiki/Farey_sequence). This time we're interested in the cardinality of $F_{1000000}$.\n",
"\n",
"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$).\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](https://en.wikipedia.org/wiki/Euler%27s_totient_function) will tell us! So the number of irreducible fractions with denominator $n+1$ is simply $\\phi(n+1)$. This gives us\n",
"The [sum of totients](https://en.wikipedia.org/wiki/Totient_summatory_function) is denoted $\\Phi(n)$.\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$.\n",
"\n",
"## Using SageMath\n",
"\n",
"Since SageMath implements a decently fast totient function, we could opt for the most straightforward approach."
"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](https://projecteuler.net/problem=70). However, there's a few more *very* interesting methods to compute totient sums.\n",
"\n",
"## Recursive totient sum\n",
"\n",
"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.\n",
"\n",
"We start with the following identity, proven by Gauss:\n",
"$$\\sum_{d|k}\\phi(d) = k$$\n",
"\n",
"Then, we can apply the formula for [triangular numbers](https://en.wikipedia.org/wiki/Triangular_number):\n",
"$$\\sum_{k=1}^n \\sum_{d|k}\\phi(d) = \\sum_{k=1}^n k = \\frac{n(n+1)}{2}$$\n",
"\n",
"Now we'll rewrite the bounds.\n",
"$$\\sum_{1 \\leq k \\leq n} \\sum_{xy=k}\\phi(x) = \\frac{n(n+1)}{2}$$\n",
"\n",
"This makes it easier to see that we can condense the two nested summations into one.\n",
"Not a super intuitive definition for $\\Phi(n)$, but pretty simple to implement, and surprisingly quick!"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "858a64d2",
"metadata": {},
"outputs": [],
"source": [
"from functools import cache\n",
"\n",
"@cache\n",
"def totient_sum(n):\n",
" if n == 1:\n",
" return 1\n",
" \n",
" return n*(n+1)/2 - sum(totient_sum(n//d) for d in range(2, n + 1))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "3f21311a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"303963552391"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"totient_sum(1000000) - 1"
]
},
{
"cell_type": "markdown",
"id": "6f1085a3",
"metadata": {},
"source": [
"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.\n",
"\n",
"## Dirichlet convolution\n",
"\n",
"The [Dirichlet convolution](https://en.wikipedia.org/wiki/Dirichlet_convolution) of two [arithmetic functions](https://en.wikipedia.org/wiki/Arithmetic_function) $f$ and $g$ is\n",
"$$(f * g)(n) = \\sum_{xy=n} f(x) g(y)$$\n",
"\n",
"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](https://en.wikipedia.org/wiki/Identity_function) and $\\mu$ is the [Möbius function](https://en.wikipedia.org/wiki/M%C3%B6bius_function) - if you've never heard of the latter, don't worry about it yet.\n",
"\n",
"## Dirichlet's hyperbola method\n",
"By applying the above Dirichlet convolution to\n",
"$$\\Phi(n) = \\sum_{k=1}^n \\phi(n)$$\n",
"we can derive an equivalent summation:\n",
"$$\\Phi(n) = \\sum_{k=1}^n \\sum_{xy=k} x \\mu(y)$$\n",
"\n",
"We're going to transform this summation by applying \n",
"[Dirichlet's hyperbola method](https://en.wikipedia.org/wiki/Dirichlet_hyperbola_method). Let $1 < a < n$, and let $b = n / a$. Then\n",
"$$\\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)$$\n",
"\n",
"Then with some algebra, we can transform this into\n",
"$$\\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}$$\n",
"where $M(n)$ is the [Mertens function](https://en.wikipedia.org/wiki/Mertens_function), which is just the partial sums of the Möbius function:\n",
"$$M(n) = \\sum_{k=1}^n \\mu(k)$$\n",
"\n",
"By using this formula for $\\Phi(n)$, instead of calculating $n$ different totients, our sums only run to $a$ and $b$.\n",
"\n",
"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](https://en.wikipedia.org/wiki/Unit_function). Once again, we choose $1 < a < n$ and let $b = n/a$.\n",
"Now we have a recursive implementation of $M(n)$. All that's left is to calculate the values of $\\mu(n)$ that we need. We can do this with a sieve, since $\\mu$ is a [multiplicative function](https://en.wikipedia.org/wiki/Multiplicative_function)."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "48c9275b",
"metadata": {},
"outputs": [],
"source": [
"from math import isqrt\n",
"\n",
"def mobius_sieve(limit):\n",
" ms = [n for n in range(0, limit)]\n",
"\n",
" for n in range(0, limit):\n",
" if n == 0 or n == 1 or ms[n] != n:\n",
" yield ms[n]\n",
" continue\n",
" \n",
" yield -1\n",
" \n",
" for k in range(2 * n, limit, n):\n",
" if k % n ^ 2 == 0:\n",
" ms[k] = 0\n",
"\n",
" ms[k] //= -n"
]
},
{
"cell_type": "markdown",
"id": "f12a8b30",
"metadata": {},
"source": [
"We'll use $a = \\sqrt{1000000} = 1000$ as our upper bound on the sieve."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c2071ba1",
"metadata": {},
"outputs": [],
"source": [
"mu = list(mobius_sieve(isqrt(1000000) + 1))"
]
},
{
"cell_type": "markdown",
"id": "a8db5c70",
"metadata": {},
"source": [
"Now we can implement $M(n)$ with our recursive formula:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "f04ed69a",
"metadata": {},
"outputs": [],
"source": [
"@cache\n",
"def M(n):\n",
" if n == 0 or n == 1:\n",
" return n\n",
" \n",
" a = sqrt(n)\n",
" b = n / a\n",
" \n",
" total = 1 + floor(b) * M(floor(a))\n",
" total -= sum(mu[x] * floor(n/x) for x in range(1, floor(a) + 1))\n",
" total -= sum(M(floor(n/y)) for y in range(2, floor(b) + 1))\n",
" \n",
" return total"
]
},
{
"cell_type": "markdown",
"id": "132dad8d",
"metadata": {},
"source": [
"And finally, yet another implementation of $\\Phi(n)$:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "394d604a",
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"def totient_sum(n):\n",
" a = sqrt(n)\n",
" b = n / a\n",
"\n",
" total = sum(x * M(floor(n/x)) for x in range(1, floor(a) + 1))\n",
" total += sum(polygonal_number(3, floor(n/y)) * mu[y] for y in range(1, floor(b) + 1))\n",
" total -= M(floor(b)) * polygonal_number(3, floor(a))\n",
"\n",
" return total"
]
},
{
"cell_type": "markdown",
"id": "e59a4350",
"metadata": {},
"source": [
"This approach is significantly faster than the others."
