From 301628bc448133ef8668c02f1e05d3ca87df0803 Mon Sep 17 00:00:00 2001
From: filifa <filifa.breath652@8alias.com>
Date: Thu, 10 Apr 2025 00:37:48 -0400
Subject: [PATCH] add problem 21

---
 problem0021.ipynb | 96 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 96 insertions(+)
 create mode 100644 problem0021.ipynb

diff --git a/problem0021.ipynb b/problem0021.ipynb
new file mode 100644
index 0000000..8825113
--- /dev/null
+++ b/problem0021.ipynb
@@ -0,0 +1,96 @@
+{
+ "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",
+    "## 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
+}