From 973da7e99d8fd41035092148b402727db2679821 Mon Sep 17 00:00:00 2001 From: filifa Date: Wed, 23 Jul 2025 20:48:58 -0400 Subject: [PATCH] reorganize and elaborate --- notebooks/problem0021.ipynb | 181 +++++++++++++++++++++++++++++++----- 1 file changed, 156 insertions(+), 25 deletions(-) diff --git a/notebooks/problem0021.ipynb b/notebooks/problem0021.ipynb index 4abf6ca..64b2ef0 100644 --- a/notebooks/problem0021.ipynb +++ b/notebooks/problem0021.ipynb @@ -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",