From 4fa83b3333f780b379d223d7e228a80f1fbe8856 Mon Sep 17 00:00:00 2001 From: filifa Date: Thu, 17 Jul 2025 01:10:47 -0400 Subject: [PATCH] expand significantly with additional approaches --- notebooks/problem0072.ipynb | 257 +++++++++++++++++++++++++++++++++++- 1 file changed, 255 insertions(+), 2 deletions(-) diff --git a/notebooks/problem0072.ipynb b/notebooks/problem0072.ipynb index 340ef76..a57c95a 100644 --- a/notebooks/problem0072.ipynb +++ b/notebooks/problem0072.ipynb @@ -11,12 +11,17 @@ "\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", "\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", + "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", "$$|F_{n+1}| = |F_n| + \\phi(n+1)$$\n", "From this, we can derive a non-recursive formula:\n", "$$|F_n| = 1 + \\sum_{k=1}^n \\phi(k)$$\n", "\n", - "As mentioned before, we'll actually subtract two from this total, since the problem isn't counting 0 or 1." + "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." ] }, { @@ -40,6 +45,254 @@ "sum(euler_phi(n) for n in range(1, 1000001)) - 1" ] }, + { + "cell_type": "markdown", + "id": "8e268820", + "metadata": {}, + "source": [ + "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", + "$$\\sum_{1 \\leq xy \\leq n} \\phi(x) = \\frac{n(n+1)}{2}$$\n", + "\n", + "We can then split this single summation into a sum-of-summations, letting $y$ vary from 1 to $n$.\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}$$\n", + "\n", + "If we apply the definition of $\\Phi(n)$, we get\n", + "$$\\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}$$\n", + "\n", + "Then we can rearrange and simplify to get\n", + "$$\\Phi(n) = \\frac{1}{2}n(n+1) - \\sum_{d=2}^n \\Phi\\left(\\left\\lfloor \\frac{n}{d} \\right\\rfloor\\right)$$\n", + "\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", + "\n", + "$$\\sum_{k=1}^n \\epsilon(k) = \\sum_{k=1}^n \\sum_{xy=k} \\mu(x)$$\n", + "\n", + "$$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)$$\n", + "\n", + "$$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)$$\n", + "\n", + "$$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)$$\n", + "\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." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "55f0719b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "303963552391" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "totient_sum(1000000) - 1" + ] + }, { "cell_type": "markdown", "id": "1c67d8da",