reorganize and elaborate
This commit is contained in:
filifa
parent
212eb04094
commit
973da7e99d
|
|
@ -9,30 +9,57 @@
|
|||
"\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."
|
||||
"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`."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "49bbab13",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def aliquot_sum(n): return sigma(n) - n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "c57c9f76",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"limit = 10000"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a2fe4975",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"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": 3,
|
||||
"id": "144e1f54",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"31626\n"
|
||||
]
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"31626"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"aliquot_sum = lambda n: sigma(n) - n\n",
|
||||
"\n",
|
||||
"amicables = set()\n",
|
||||
"for a in range(2, 10000):\n",
|
||||
"for a in range(2, limit):\n",
|
||||
" if a in amicables:\n",
|
||||
" continue\n",
|
||||
" \n",
|
||||
|
|
@ -42,8 +69,37 @@
|
|||
" \n",
|
||||
" if a == aliquot_sum(b):\n",
|
||||
" amicables.update({a, b})\n",
|
||||
"\n",
|
||||
"print(sum(amicables))"
|
||||
" \n",
|
||||
"sum(amicables)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "5b540764",
|
||||
"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!"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "f6014a35",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{220, 284, 1184, 1210, 2620, 2924, 5020, 5564, 6232, 6368}"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"amicables"
|
||||
]
|
||||
},
|
||||
{
|
||||
|
|
@ -51,8 +107,6 @@
|
|||
"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",
|
||||
|
|
@ -67,31 +121,108 @@
|
|||
"\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$."
|
||||
"## Avoiding factoring\n",
|
||||
"Since we need all the sums of divisors up to 10000, instead of factoring each number individually, we can compute each value of $\\sigma$ in a loop, once again taking advantage of $\\sigma$ being multiplicative. (In some ways, this method is similar to the sieve of Eratosthenes - see [problem 10](https://projecteuler.net/problem=10))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"execution_count": 5,
|
||||
"id": "4cbc68db",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def sum_of_divisors_sieve(limit):\n",
|
||||
" dsum = [1 for _ in range(0, limit)]\n",
|
||||
" \n",
|
||||
"def sum_of_divisors_range(limit): \n",
|
||||
" dsum = [1 for n in range(0, limit)]\n",
|
||||
" dsum[0] = 0\n",
|
||||
" \n",
|
||||
" for n in range(0, limit):\n",
|
||||
" if n == 0 or n == 1:\n",
|
||||
" if n == 0 or n == 1 or dsum[n] != 1:\n",
|
||||
" # n is 0, 1, or composite\n",
|
||||
" yield dsum[n]\n",
|
||||
" continue\n",
|
||||
"\n",
|
||||
" # n is prime\n",
|
||||
" dsum[n] = n + 1\n",
|
||||
"\n",
|
||||
" m = n\n",
|
||||
" while True:\n",
|
||||
" # sigma(a*b) = sigma(a) * sigma(b) if gcd(a, b) = 1\n",
|
||||
" for k in range(2 * m, limit, m):\n",
|
||||
" if k % (m*n) != 0:\n",
|
||||
" dsum[k] *= dsum[m]\n",
|
||||
"\n",
|
||||
" if m * n >= limit:\n",
|
||||
" break\n",
|
||||
" \n",
|
||||
" # sigma(p^k) = p^k + sigma(p^(k-1))\n",
|
||||
" dsum[m*n] = m*n + dsum[m]\n",
|
||||
" m *= n\n",
|
||||
" \n",
|
||||
" for k in range(n, limit, n):\n",
|
||||
" dsum[k] += n\n",
|
||||
" \n",
|
||||
" yield dsum[n]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "f4dc20fe",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"With this method, we can redefine our aliquot sum function:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"id": "6b5cb5f8",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"divisor_sums = list(sum_of_divisors_range(limit))\n",
|
||||
"def aliquot_sum(n): return divisor_sums[n] - n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "bf501433",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Then we compute the amicable numbers the same as before."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"id": "2f70d28c",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"31626"
|
||||
]
|
||||
},
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"amicables = set()\n",
|
||||
"for a in range(2, limit):\n",
|
||||
" if a in amicables:\n",
|
||||
" continue\n",
|
||||
" \n",
|
||||
" b = aliquot_sum(a)\n",
|
||||
" if a == b:\n",
|
||||
" continue\n",
|
||||
" \n",
|
||||
" # if b is greater than limit, it will cause an IndexError\n",
|
||||
" if b < limit and a == aliquot_sum(b):\n",
|
||||
" amicables.update({a, b})\n",
|
||||
" \n",
|
||||
"sum(amicables)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a847f754",
|
||||
|
|
|
|||
Loading…
Reference in New Issue