diff --git a/notebooks/problem0021.ipynb b/notebooks/problem0021.ipynb index 8825113..e87c765 100644 --- a/notebooks/problem0021.ipynb +++ b/notebooks/problem0021.ipynb @@ -67,6 +67,61 @@ "\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)" ]