filifa/eulerbooks
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 3 Wiki Activity

simplify dirichlet convolution approach

Browse Source
This commit is contained in:
filifa
2026-05-11 22:31:08 -04:00
parent ef94059ff6
commit 171a0e106c
1 changed files with 24 additions and 110 deletions
Download Patch File Download Diff File Expand all files Collapse all files

134
notebooks/problem0072.ipynb
View File

@@ -130,137 +130,51 @@
"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", "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", "$$(f * g)(n) = \\sum_{xy=n} f(x) g(y)$$\n",
"\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", "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 * 1 = \\mathrm{Id}$, where $1$ is just the constant function $1(n) = 1$ and $\\mathrm{Id}(n) = n$ is the [identity function](https://en.wikipedia.org/wiki/Identity_function) (note that this convolution identity is the same as Gauss's identity presented above).\n",
"\n", "\n",
"## Dirichlet's hyperbola method\n", "## Dirichlet's hyperbola method\n",
"By applying the above Dirichlet convolution to\n", "As mentioned above, it's well-known that\n",
"$$\\Phi(n) = \\sum_{k=1}^n \\phi(n)$$\n", "$$\\sum_{k=1}^n k = \\frac{n(n+1)}{2}$$\n",
"we can derive an equivalent summation:\n", "We can go further by applying the Dirichlet convolution above to obtain\n",
"$$\\Phi(n) = \\sum_{k=1}^n \\sum_{xy=k} x \\mu(y)$$\n", "$$\\frac{n(n+1)}{2} = \\sum_{k=1}^n k = \\sum_{k=1}^n \\sum_{xy=k} \\phi(x)$$\n",
"\n", "\n",
"We're going to transform this summation by applying \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", "[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", "$$\\frac{n(n+1)}{2} = \\sum_{x=1}^a \\sum_{y=1}^{n/x} \\phi(x) + \\sum_{y=1}^b \\sum_{x=1}^{n/y} \\phi(x) - \\sum_{x=1}^a \\sum_{y=1}^b \\phi(x)$$\n",
"\n", "\n",
"Then with some algebra, we can transform this into\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", "$$\\frac{n(n+1)}{2} = \\sum_{x=1}^a \\phi(x)\\left\\lfloor\\frac{n}{x}\\right\\rfloor + \\sum_{y=1}^b \\Phi\\left(\\left\\lfloor\\frac{n}{y}\\right\\rfloor\\right) - \\lfloor b \\rfloor\\Phi(\\lfloor a \\rfloor)$$\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", "\n",
"By using this formula for $\\Phi(n)$, instead of calculating $n$ different totients, our sums only run to $a$ and $b$.\n", "This might seem complicated, but note that in the second summation on the right, when $y=1$, we have $\\Phi(n)$. We can extract that term and rearrange to get a recursive formula for $\\Phi(n)$:\n",
"$$\\Phi(n) = \\frac{n(n+1)}{2} + \\lfloor b \\rfloor\\Phi(\\lfloor a \\rfloor) - \\sum_{x=1}^a \\phi(x)\\left\\lfloor\\frac{n}{x}\\right\\rfloor - \\sum_{y=2}^b \\Phi\\left(\\left\\lfloor\\frac{n}{y}\\right\\rfloor\\right)$$\n",
"\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", "The major advantage of this formula is that we only need the values of the totient function up to $a$, instead of up to $n$. We can set $a = \\sqrt{1000000} = 1000$.\n",
"\n", "\n",
"$$\\sum_{k=1}^n \\epsilon(k) = \\sum_{k=1}^n \\sum_{xy=k} \\mu(x)$$\n", "This all leads to this code, yet another implementation of $\\Phi(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. By taking advantage of $\\mu$ being [multiplicative](https://en.wikipedia.org/wiki/Multiplicative_function), we can compute values using a similar strategy to the sieve of Eratosthenes - see [problem 10](https://projecteuler.net/problem=10)."
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 4, "execution_count": 4,
"id": "48c9275b",
"metadata": {},
"outputs": [],
"source": [
"def mobius_range(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."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c2071ba1",
"metadata": {},
"outputs": [],
"source": [
"from math import isqrt\n",
"mu = list(mobius_range(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", "id": "394d604a",
"metadata": { "metadata": {
"scrolled": true "scrolled": true
}, },
"outputs": [], "outputs": [],
"source": [ "source": [
"@cache\n",
"def totient_sum(n):\n", "def totient_sum(n):\n",
" if n == 1:\n",
" return 1\n",
"\n",
" a = sqrt(n)\n", " a = sqrt(n)\n",
" b = n / a\n", " b = n / a\n",
"\n", "\n",
" total = sum(x * M(floor(n/x)) for x in range(1, floor(a) + 1))\n", " W = (n*(n+1)) // 2\n",
" total += sum(polygonal_number(3, floor(n/y)) * mu[y] for y in range(1, floor(b) + 1))\n", " X = floor(b) * totient_sum(floor(a))\n",
" total -= M(floor(b)) * polygonal_number(3, floor(a))\n", " Y = sum(euler_phi(x) * (n//x) for x in range(1, floor(a)+1))\n",
"\n", " Z = sum(totient_sum(n//y) for y in range(2, floor(b)+1))\n",
" return total" " return W + X - Y - Z"
] ]
}, },
{ {
@@ -273,7 +187,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 8, "execution_count": 5,
"id": "55f0719b", "id": "55f0719b",
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
@@ -283,7 +197,7 @@
"303963552391" "303963552391"
] ]
}, },
"execution_count": 8, "execution_count": 5,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "execute_result"
} }
@@ -309,7 +223,7 @@
], ],
"metadata": { "metadata": {
"kernelspec": { "kernelspec": {
"display_name": "SageMath 9.5", "display_name": "SageMath 10.8",
"language": "sage", "language": "sage",
"name": "sagemath" "name": "sagemath"
}, },
@@ -323,7 +237,7 @@
"name": "python", "name": "python",
"nbconvert_exporter": "python", "nbconvert_exporter": "python",
"pygments_lexer": "ipython3", "pygments_lexer": "ipython3",
"version": "3.11.2" "version": "3.12.5"
} }
}, },
"nbformat": 4, "nbformat": 4,