{
 "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`."
   ]
  },
  {
   "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": [
    {
     "data": {
      "text/plain": [
       "31626"
      ]
     },
     "execution_count": 3,
     "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 a == aliquot_sum(b):\n",
    "        amicables.update({a, b})\n",
    "        \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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca7d5f36",
   "metadata": {},
   "source": [
    "## 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",
    "## 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": 5,
   "id": "4cbc68db",
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_multiples(dsum, p, limit):\n",
    "    q = p\n",
    "    while True:\n",
    "        # sigma(a*b) = sigma(a) * sigma(b) if gcd(a, b) = 1\n",
    "        for k in range(2 * q, limit, q):\n",
    "            if k % (p*q) != 0:\n",
    "                dsum[k] *= dsum[q]\n",
    "\n",
    "        if p * q >= limit:\n",
    "            break\n",
    "\n",
    "        # sigma(p^k) = p^k + sigma(p^(k-1))\n",
    "        dsum[p*q] = p * q + dsum[q]\n",
    "        q *= p\n",
    "        \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 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",
    "        update_multiples(dsum, n, limit)\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",
   "metadata": {},
   "source": [
    "## Relevant sequences\n",
    "* Sums of divisors: [A000203](https://oeis.org/A000203)\n",
    "* Amicable numbers: [A063990](https://oeis.org/A063990)\n",
    "\n",
    "#### Copyright (C) 2025 filifa\n",
    "\n",
    "This work is licensed under the [Creative Commons Attribution-ShareAlike 4.0 International license](https://creativecommons.org/licenses/by-sa/4.0/) and the [BSD Zero Clause license](https://spdx.org/licenses/0BSD.html)."
   ]
  }
 ],
 "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
}