ipynb/Euler's Conjecture.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Euler's Sum of Powers Conjecture\n",
    "\n",
    "In 1769, Leonhard Euler [conjectured that](https://en.wikipedia.org/wiki/Euler%27s_sum_of_powers_conjecture) for all integers *n* and *k* greater than 1, if the sum of *n* *k*th powers of positive integers is itself a *k*th power, then *n* is greater than or equal to *k*. For example, this would mean that no sum of a pair of cubes (a<sup>3</sup> + b<sup>3</sup>) can be equal to another cube (c<sup>3</sup>), but a sum of three cubes can, as in 3<sup>3</sup> + 4<sup>3</sup> + 5<sup>3</sup> = 6<sup>3</sup>. \n",
    "\n",
    "It took 200 years to disprove the conjecture: in 1966 L. J. Lander and T. R. Parkin published a refreshingly short [article](https://projecteuclid.org/download/pdf_1/euclid.bams/1183528522) giving a counterexample of four fifth-powers that summed to another fifth power. They found it via a program that did an exhaustive search.  Can we duplicate their work and find integers greater than 1 such that \n",
    "*a*<sup>5</sup> + *b*<sup>5</sup> + *c*<sup>5</sup> + *d*<sup>5</sup> = *e*<sup>5</sup> ?\n",
    "\n",
    "## Algorithm\n",
    "\n",
    "An exhaustive *O*(*m*<sup>4</sup>) algorithm woud be to look at all values of *a, b, c, d* < *m* and check if *a*<sup>5</sup> + *b*<sup>5</sup> + *c*<sup>5</sup> + *d*<sup>5</sup>  is a fifth power. But we can do better: a sum of four numbers is a sum of two pairs of numbers, so we\n",
    "are looking for\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*pair*<sub>1</sub> + *pair*<sub>2</sub> = *e*<sup>5</sup> &nbsp;&nbsp; **where** &nbsp; *pair*<sub>1</sub> =  *a*<sup>5</sup> + *b*<sup>5</sup> **and**  *pair*<sub>2</sub> = *c*<sup>5</sup> + *d*<sup>5</sup>.\n",
    "\n",
    "We will define *pairs* be a dict of `{`*a*<sup>5</sup> + *b*<sup>5</sup>`: (`*a*<sup>5</sup>`, ` *b*<sup>5</sup>`)}` entries for all *a*  ≤ *b* < *m*; for example, for *a*=2 and *b*=10, the entry is  `{100032: (32, 100000)}`.\n",
    "Then we can ask for each *pair*<sub>1</sub>, and for each *e*, whether there is a *pair*<sub>2</sub> in the `dict` that makes the equation work. There are *O*(*m*<sup>2</sup>) pairs and *O*(*m*) values of *e*, and `dict` lookup is *O*(1), so the whole algorithm is *O*(*m*<sup>3</sup>):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "\n",
    "def euler(m):\n",
    "    \"\"\"Yield tuples (a, b, c, d, e) such that a^5 + b^5 + c^5 + d^5 = e^5,\n",
    "    where all are integers, and 1 < a ≤ b ≤ c ≤ d < e < m.\"\"\"\n",
    "    powers = [e**5 for e in range(2, m)] \n",
    "    pairs  = {sum(pair): pair \n",
    "              for pair in itertools.combinations_with_replacement(powers, 2)}\n",
    "    for pair1 in pairs:\n",
    "        for e5 in powers:\n",
    "            pair2 = e5 - pair1\n",
    "            if pair2 in pairs:\n",
    "                yield fifthroots(pairs[pair1] + pairs[pair2] + (e5,))\n",
    "    \n",
    "def fifthroots(nums): \n",
    "    \"Sorted integer fifth roots of a collection of numbers.\" \n",
    "    return tuple(sorted(int(round(x ** (1/5))) for x in nums))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look for a solution (arbitrarily choosing *m*=500):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.07 s, sys: 21.4 ms, total: 1.09 s\n",
      "Wall time: 1.11 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(27, 84, 110, 133, 144)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time next(euler(500))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That was easy, and it turns out this is the same answer that Lander and Parkin got: 27<sup>5</sup> + 84<sup>5</sup> +  110<sup>5</sup> +  133<sup>5</sup> = 144<sup>5</sup>.\n",
    "\n",
    "We can keep going, collecting all the solutions up to `*m*=1000`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 53s, sys: 706 ms, total: 1min 54s\n",
      "Wall time: 1min 57s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{(27, 84, 110, 133, 144),\n",
       " (54, 168, 220, 266, 288),\n",
       " (81, 252, 330, 399, 432),\n",
       " (108, 336, 440, 532, 576),\n",
       " (135, 420, 550, 665, 720),\n",
       " (162, 504, 660, 798, 864)}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time set(euler(1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the answers are multiples of the first one (this is easiest to see in the middle column: 110, 220, 330, ...).\n",
    "Since 1966 other mathematicians have found [other solutions](https://en.wikipedia.org/wiki/Euler%27s_sum_of_powers_conjecture), but all we need is one to disprove Euler's conjecture."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}