152 lines
5.2 KiB
Plaintext
152 lines
5.2 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "dea82444",
|
|
"metadata": {},
|
|
"source": [
|
|
"# [Amicable Numbers](https://projecteuler.net/problem=21)\n",
|
|
"\n",
|
|
"The sum of the proper divisors of a number is called the [aliquot sum](https://en.wikipedia.org/wiki/Aliquot_sum).\n",
|
|
"\n",
|
|
"SageMath [provides the `sigma` function](https://doc.sagemath.org/html/en/reference/rings_standard/sage/arith/misc.html#sage.arith.misc.Sigma), which can compute the sum of the divisors of $n$. The aliquot sum is then just `sigma(n) - n`.\n",
|
|
"\n",
|
|
"If a number is equal to its own aliquot sum, it's called a [perfect number](https://en.wikipedia.org/wiki/Perfect_number), and we exclude those numbers from our total."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"id": "144e1f54",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"31626\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"aliquot_sum = lambda n: sigma(n) - n\n",
|
|
"\n",
|
|
"amicables = set()\n",
|
|
"for a in range(2, 10000):\n",
|
|
" if a in amicables:\n",
|
|
" continue\n",
|
|
" \n",
|
|
" b = aliquot_sum(a)\n",
|
|
" if a == b:\n",
|
|
" continue\n",
|
|
" \n",
|
|
" if a == aliquot_sum(b):\n",
|
|
" amicables.update({a, b})\n",
|
|
"\n",
|
|
"print(sum(amicables))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "ca7d5f36",
|
|
"metadata": {},
|
|
"source": [
|
|
"Funny enough, there's only five pairs of amicable numbers below 10,000. If you looked up [amicable numbers](https://en.wikipedia.org/wiki/Amicable_numbers), you may have stumbled on all the numbers you need to answer the problem!\n",
|
|
"\n",
|
|
"## Sum of divisors\n",
|
|
"Of course, you could implement your own [divisor sum function](https://en.wikipedia.org/wiki/Divisor_function). In [problem 12](https://projecteuler.net/problem=12), we implemented a divisor *counting* function, which is related.\n",
|
|
"\n",
|
|
"One important property of the divisor sum function $\\sigma(n)$ is that it is [multiplicative](https://en.wikipedia.org/wiki/Multiplicative_function). This means that if $a$ and $b$ are [coprime](https://en.wikipedia.org/wiki/Coprime_integers), $\\sigma(ab) = \\sigma(a)\\sigma(b)$. It follows that if $n = 2^{r_1}3^{r_2}5^{r_3}7^{r_4}\\cdots$, then\n",
|
|
"$$\\sigma(n) = \\sigma(2^{r_1})\\sigma(3^{r_2})\\sigma(5^{r_3})\\sigma(7^{r_4})\\cdots$$\n",
|
|
"\n",
|
|
"Furthermore, for a prime $p$, $\\sigma(p^k) = 1 + p + p^2 + \\cdots + p^k$. This expression is actually a partial sum of a [geometric series](https://en.wikipedia.org/wiki/Geometric_series), which has a closed formula:\n",
|
|
"$$1 + p + p^2 + \\cdots + p^k = \\frac{p^{k+1}-1}{p-1}$$\n",
|
|
"\n",
|
|
"Putting it all together, we can compute $\\sigma(n)$ as\n",
|
|
"$$\\sigma(n) = \\frac{2^{r_1+1}-1}{2-1} \\cdot \\frac{3^{r_2+1}-1}{3-1} \\cdot \\frac{5^{r_3+1}-1}{5-1} \\cdot \\frac{7^{r_4+1}-1}{7-1} \\cdot \\cdots$$\n",
|
|
"\n",
|
|
"Therefore, if you have the number's factorization (see [problem 3](https://projecteuler.net/problem=3)), you can use it to compute the sum of its divisors.\n",
|
|
"\n",
|
|
"## Sieving the sums of divisors\n",
|
|
"Since we need all the sums of divisors up to 10000, instead of factoring each number individually, we could sieve the values of $\\sigma$. This is possible because the sum of divisors function is multiplicative."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"id": "4cbc68db",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def sum_of_divisors_sieve(limit):\n",
|
|
" dsum = [1 for _ in range(0, limit)]\n",
|
|
" \n",
|
|
" for n in range(0, limit):\n",
|
|
" if n == 0 or n == 1 or dsum[n] != 1:\n",
|
|
" yield dsum[n]\n",
|
|
" continue\n",
|
|
" \n",
|
|
" m = 1\n",
|
|
" while m * n < limit:\n",
|
|
" dsum[m*n] = dsum[m] + m*n\n",
|
|
" for k in range(2 * m * n, limit, m * n):\n",
|
|
" if k % (m*n^2) != 0:\n",
|
|
" dsum[k] *= dsum[m*n]\n",
|
|
" \n",
|
|
" m *= n\n",
|
|
" \n",
|
|
" yield dsum[n]"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "fdbabc18",
|
|
"metadata": {},
|
|
"source": [
|
|
"Here's a check to make sure the sieve matches the values from SageMath."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"id": "0539d1fc",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"s = list(sum_of_divisors_sieve(10000))\n",
|
|
"assert all(s[n] == sigma(n) for n in range(1, 10000))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "a847f754",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Relevent sequences\n",
|
|
"* Amicable numbers: [A063990](https://oeis.org/A063990)"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "SageMath 9.5",
|
|
"language": "sage",
|
|
"name": "sagemath"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.11.2"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 5
|
|
}
|