ipynb/Beal.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "<div style=\"text-align:right\"><b>Peter Norvig</b> 2000; revised 2015&mdash;2018</div>\n",
    "\n",
    "# Beal's Conjecture Revisited\n",
    "\n",
    "In 1637, Pierre de Fermat wrote in the margin of a book that he had a proof of his famous \"[Last Theorem](https://en.wikipedia.org/wiki/Fermat%27s_Last_Theorem)\":\n",
    "\n",
    "> If $A^n + B^n = C^n$,\n",
    "> <br>where $A, B, C, n$ are positive integers\n",
    "> <br>then $n \\le 2$.\n",
    "\n",
    "Centuries passed before [Andrew Beal](https://en.wikipedia.org/wiki/Andrew_Beal), a businessman and amateur mathematician,\n",
    "made his conjecture in 1993:\n",
    "\n",
    "> If $A^x + B^y = C^z$, \n",
    "> <br>where $A, B, C, x, y, z$ are positive integers \n",
    "> <br>and $x, y, z$ are all greater than $2$, \n",
    "> <br>then $A, B$ and $C$ must have a common prime factor.\n",
    "\n",
    "[Andrew Wiles](https://en.wikipedia.org/wiki/Andrew_Wiles) proved Fermat's theorem in 1995, but Beal's conjecture remains unproved, and Beal has offered [one million dollars](http://www.ams.org/profession/prizes-awards/ams-supported/beal-prize) for a proof or disproof. I don't have the mathematical skills of Wiles, so I could never find a proof, but I can write a program to search for  counterexamples. I first wrote [that program in 2000](http://norvig.com/beal2000.html), and [my name got associated](https://www.google.com/webhp?#q=beal conjecture) with Beal's Conjecture, which means I get a lot of emails with purported proofs or counterexamples (many asking how they can collect their prize money). So far, all the emails have been wrong. This notebook catalogs some of the more common errors, updates my 2000 program, and introduces this tool for verifying counterexamples:\n",
    " \n",
    "\n",
    "| [**Online Beal Counterexample Checker**](http://norvig.com/bealcheck.html)  |\n",
    "|:---:|\n",
    "\n",
    "\n",
    "# How to Not Win A Million Dollars\n",
    "\n",
    "\n",
    "* A proof must show that there are **no** examples that satisfy the conditions. A common error is to show how a certain pattern generates an infinite number of $(A, x, B, y, C, z)$ examples, and that the conjecture holds for this entire infinite collection.  But that's not good enough, unless you can also prove that the conjecture holds for every other possible pattern.\n",
    "\n",
    "\n",
    "* It is valid to use proof by contradiction: assume the conjecture is true, and show that that leads to a contradiction. It is not valid to use proof by circular reasoning: assume the conjecture is true, put in some irrelevant steps, and show that it follows that the conjecture is true.\n",
    "\n",
    "\n",
    "* A valid counterexample needs to satisfy all four conditions&mdash;don't leave one out.\n",
    "\n",
    "\n",
    "* One correspondent claimed that $27^4 + 162 ^ 3 = 9 ^ 7$ was a solution, because the first three conditions hold, and the common factor is 9, which isn't a prime.  But of course, if $A, B, C$ have 9 as a common factor, then they also have 3, and 3 is prime. \"No common prime factor\" means the same thing as \"no common factor greater than 1.\"\n",
    "\n",
    "\n",
    "* Another claimed that $2^3+2^3=2^4$ was a counterexample, because all the bases are 2, which is prime, and prime numbers have no prime factors. But that's not true; a prime number has itself as a factor.\n",
    "\n",
    "\n",
    "* A creative person offered $ 1359072^4 - 940896^4 = 137998080^3$, which fails both because $ 3^3 2^5 11^2 $ is a common factor, and because it has a subtraction rather than an addition (although, as Julius Jacobsen pointed out, it could be rewritten as $ 137998080^3 + 940896^4 = 1359072^4 $).\n",
    "\n",
    "\n",
    "* Mustafa Pehlivan came up with an example involving 76-million-digit numbers, which took some work to prove wrong (using modulo arithmetic). \n",
    "\n",
    "\n",
    "* Another Beal fan started by saying \"Let $C = 43$ and $z = 3$.  Since $43 = 21 + 22$, we have $43^3 = (21^3 + 22^3)$.\" But of course $(a + b)^3 \\ne (a^3 + b^3)$. This fallacy is called [the freshman's dream](https://en.wikipedia.org/wiki/Freshman%27s_dream) (although I remember having different dreams as a freshman).\n",
    "\n",
    "\n",
    "* Multiple people proposed counterexamples similar to this one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from math import gcd #### In Python versions < 3.5, use \"from fractions import gcd\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "button": false,
    "collapsed": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A, B, C = 60000000000000000001, 70000000000000000003, 82376613842809255677\n",
    "x = y = z = 3.\n",
    "\n",
    "A ** x + B ** y == C ** z and gcd(A, B) == gcd(B, C) == 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "**WOW! The result is `True`!** The two sides of the equation are equal, and the greatest common divisor is 1. Is this a real counterexample to Beal? And also a disproof of Fermat's Last Theorem?\n",
    "\n",
    "Alas, it is not.  The decimal point in \"`x = y = z = 3`**`.`**\" indicates a floating point number, with inexact, limited precision. Change the inexact \"`3.`\" to an exact \"`3`\"  and the two sides of the equation are no longer equal. Below we see they are the same for the first 19 digits, but differ starting with the 20th: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "button": false,
    "collapsed": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(559000000000000000054900000000000000002070000000000000000028,\n",
       " 559000000000000000063037470301555182935702892172500189973733)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A ** 3 + B ** 3,\n",
    " C ** 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "They say \"close\" only counts in horseshoes and hand grenades, and if you stood in your yard and threw a horseshoe at a stake on **[Kapteyn-b](https://en.wikipedia.org/wiki/Kapteyn_b)** (an exoplanet 12.8 light years from Earth that is deemed habitable and thus possibly horseshoe-playing) and the flight path differed from the perfect path in the 20th digit, then it would end up \n",
    "[about a millimeter](https://www.google.com/search?q=12.8+light+years+%2F+10%5E20+in+inches&oq=12.8+light+years+%2F+10%5E19+in+inches) \n",
    "from the target. That's really, really close, but close doesn't count in number theory.\n",
    "\n",
    "![Kapteyn-b and Homer Simpson](https://norvig.com/homer.png)\n",
    "Left: [Kapteyn-b](https://www.space.com/26115-oldest-habitable-alien-planet-kapteyn-b.html). &nbsp; &nbsp; Right: [Homer Simpson](https://www.youtube.com/watch?time_continue=1&v=ReOQ300AcSU).\n",
    "\n",
    "\n",
    "# *The Simpsons* and Fermat\n",
    "\n",
    "In two different [episodes of *The Simpsons*](http://www.npr.org/sections/krulwich/2014/05/08/310818693/did-homer-simpson-actually-solve-fermat-s-last-theorem-take-a-look), close counterexamples to Fermat's Last Theorem are shown: \n",
    "$3987^{12} + 4365^{12} = 4472^{12}$ and $1782^{12} + 1841^{12} = 1922^{12}$. These were designed by *Simpsons* writer David X. Cohen to be correct up to the precision of a typical handheld calculator; here we see the two sides of the second equation agree on the first ten digits, `6397665634`, and then differ:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(63976656349698612616236230953154487896987106,\n",
       " 63976656348486725806862358322168575784124416)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "3987 ** 12 + 4365 ** 12, 4472 ** 12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cohen must have found the equations with a program something like this (here `bases` is a sequence of integers to consider for the values of `A` and `B`; the variables `An` and `Bn` hold the `A**n` and `B**n` values; `lhs` is their sum (the left-hand-side of the equation); and the function `Cn` computes the `C**n` that is closest to that sum):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from itertools import combinations\n",
    "\n",
    "def simpsons(bases, n):\n",
    "    \"\"\"Print the (A**n + B**n = C**n) equation that minimizes the relative error,\n",
    "    for a given n and A, B values from the sequence of integers `bases`.\"\"\"\n",
    "    def Cn(lhs):  return iroot(sum(lhs), n) ** n\n",
    "    def err(lhs): return abs(sum(lhs) - Cn(lhs)) / sum(lhs)\n",
    "    def show(Xn): return '{} ** {}'.format(iroot(Xn, n), n)\n",
    "    powers = [b ** n for b in bases]\n",
    "    (An, Bn) = lhs = min(combinations(powers, 2), key=err)\n",
    "    print('{} + {} == {} (with error {:.0g})'\n",
    "          .format(show(An), show(Bn), show(Cn(lhs)), err(lhs)))\n",
    "\n",
    "def iroot(x, n): \"integer nth root\"; return int(round(x ** (1 / n)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1782 ** 12 + 1841 ** 12 == 1922 ** 12 (with error 3e-10)\n",
      "3987 ** 12 + 4365 ** 12 == 4472 ** 12 (with error 2e-11)\n"
     ]
    }
   ],
   "source": [
    "simpsons(range(1000, 2000), 12)\n",
    "simpsons(range(2000, 5000), 12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the same two equations that David X. Cohen found. \n",
    "\n",
    "Can we find other near-misses? I'll try each single-digit exponent. I want A, B, C to be 4 digits each, so I'll limit A and B to 9500 (not 9999), to try to keep C from overflowing to 5 digits. (This takes around 10 minutes to run.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5856 ** 3 + 9036 ** 3 == 9791 ** 3 (with error 1e-12)\n",
      "2396 ** 4 + 4551 ** 4 == 4636 ** 4 (with error 4e-11)\n",
      "3993 ** 5 + 7767 ** 5 == 7822 ** 5 (with error 2e-11)\n",
      "6107 ** 6 + 8919 ** 6 == 9066 ** 6 (with error 8e-13)\n",
      "5592 ** 7 + 9079 ** 7 == 9122 ** 7 (with error 2e-11)\n",
      "4749 ** 8 + 8952 ** 8 == 8959 ** 8 (with error 3e-11)\n",
      "5433 ** 9 + 6725 ** 9 == 6828 ** 9 (with error 4e-11)\n"
     ]
    }
   ],
   "source": [
    "for n in range(3, 10):\n",
    "    simpsons(range(1000, 9500), n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "The equation for *n*=6 has the smallest error yet (in the 12th decimal place).\n",
    "\n",
    "# Back to `beal`\n",
    "\n",
    "In October 2015 I looked back at my original [program from 2000](http://norvig.com/beal2000.html).\n",
    "I ported it from Python 1.5 to  3.5 (`print` is now a function, `long` is `int`). It runs 250 times faster today, a tribute to both computer hardware engineers and the developers of the Python interpreter.\n",
    "\n",
    "I found that I had [misstated the problem](https://web.archive.org/web/19991127081319/http://bealconjecture.com/) in 2000.  I thought that, *by definition*, $A$ and $B$ could not have a common factor, but actually, \n",
    "the conjecture only rules out examples where all three of $A, B, C$ share a common factor. But, as [Mark Tiefenbruck](mailto:mark @tiefenbruck.org) (as well as Edward P. Berlin and Shen Lixing) pointed out, my statement is correct, not by definition, but *by derivation:* if $A$ and $B$ have a common prime factor $p$, then the sum of $A^x + B^y$ must also have that factor $p$, and hence $C^z$, and $C$, must have the factor $p$.\n",
    "\n",
    "Mark Tiefenbruck also suggested another optimization: only consider exponents that are odd primes, or 4.  The idea is that a number like 512 can be expressed as either $2^9$ or $8^3$, and my program doesn't need to consider both. In general, any time we have a composite exponent, such as $b^{qp}$, where $p$ is prime, we should ignore $A=b, x=qp$, and instead consider only $A=b^q, x=p$. There's one complication to this scheme: 2 is a prime, but 2 is not a valid exponent for a Beal counterexample.  So we will allow 4 as an exponent, as well as all odd primes up to `max_x`.\n",
    "\n",
    "Here is the complete, updated program:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "from math      import gcd, log\n",
    "from itertools import combinations, product\n",
    "\n",
    "def beal(max_A, max_x):\n",
    "    \"\"\"See if any A ** x + B ** y equals some C ** z, with gcd(A, B) == 1.\n",
    "    Consider any 1 <= A,B <= max_A and x,y <= max_x, with x,y prime or 4.\"\"\"\n",
    "    Apowers = make_Apowers(max_A, max_x)\n",
    "    Czroots = make_Czroots(Apowers)\n",
    "    for (A, B) in combinations(Apowers, 2):\n",
    "        if gcd(A, B) == 1:\n",
    "            for (Ax, By) in product(Apowers[A], Apowers[B]):       \n",
    "                Cz = Ax + By\n",
    "                if Cz in Czroots:\n",
    "                    C = Czroots[Cz]\n",
    "                    x, y, z = exponent(Ax, A), exponent(By, B), exponent(Cz, C)\n",
    "                    print('{} ** {} + {} ** {} == {} ** {} == {}'\n",
    "                          .format(A, x, B, y, C, z, C ** z))\n",
    "\n",
    "def make_Apowers(max_A, max_x): \n",
    "    \"A dict of {A: [A**3, A**4, ...], ...}.\"\n",
    "    exponents = exponents_upto(max_x)\n",
    "    return {A: [A ** x for x in (exponents if (A != 1) else [3])]\n",
    "            for A in range(1, max_A+1)}\n",
    "\n",
    "def make_Czroots(Apowers): return {Cz: C for C in Apowers for Cz in Apowers[C]}            \n",
    "    \n",
    "def exponents_upto(max_x):\n",
    "    \"Return all odd primes up to max_x, as well as 4.\"\n",
    "    exponents = [3, 4] if max_x >= 4 else [3] if max_x == 3 else []\n",
    "    for x in range(5, max_x, 2):\n",
    "        if not any(x % p == 0 for p in exponents):\n",
    "            exponents.append(x)\n",
    "    return exponents\n",
    "\n",
    "def exponent(Cz, C): \n",
    "    \"\"\"Recover z such that C ** z == Cz (or equivalently z = log Cz base C).\n",
    "    For exponent(1, 1), arbitrarily choose to return 3.\"\"\"\n",
    "    return 3 if (Cz == C == 1) else int(round(log(Cz, C)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "It takes less than a second to verify that there are no counterexamples for combinations up to $100^{100}$, a computation that took Andrew Beal thousands of hours on his 1990s-era computers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "button": false,
    "collapsed": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 353 ms, sys: 4.84 ms, total: 358 ms\n",
      "Wall time: 376 ms\n"
     ]
    }
   ],
   "source": [
    "%time beal(100, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "The execution time goes up roughly with the square of `max_A`, so the following should take about 25 times longer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "button": false,
    "collapsed": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 8.97 s, sys: 56.6 ms, total: 9.03 s\n",
      "Wall time: 9.12 s\n"
     ]
    }
   ],
   "source": [
    "%time beal(500, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# How `beal`  Works\n",
    "\n",
    "The function `beal` first does some precomputation, creating two data structures:\n",
    "* `Apowers`: a dict of the form `{A: [A**3, A**4, ...]}` giving the\n",
    "nonredundant powers (prime and 4th powers) of each base, `A`, from 3 to `max_x`.\n",
    "* `Czroots`: a dict of `{C**z : C}` pairs, giving the zth root of each power in `Apowers`.\n",
    "\n",
    "Here is a very small example Apowers table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "button": false,
    "collapsed": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{1: [1],\n",
       " 2: [8, 16, 32, 128],\n",
       " 3: [27, 81, 243, 2187],\n",
       " 4: [64, 256, 1024, 16384],\n",
       " 5: [125, 625, 3125, 78125],\n",
       " 6: [216, 1296, 7776, 279936]}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Apowers = make_Apowers(6, 10)\n",
    "Apowers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Then we enumerate all combinations of two bases, `A` and `B`, from `Apowers`. Consider the combination where `A` is `3` and `B` is `6`. Of course `gcd(3, 6) == 3`, so the program would not consider them further, but imagine if they did not share a common factor. Then we would look at all possible `Ax + By` sums, for `Ax` in `[27, 81, 243, 2187]` and `By` in `[216, 1296, 7776, 279936].` One of these would be `27 + 216`, which sums to  `243`. We look up `243` in `Czroots`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "button": false,
    "collapsed": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{1: 1,\n",
       " 8: 2,\n",
       " 16: 2,\n",
       " 27: 3,\n",
       " 32: 2,\n",
       " 64: 4,\n",
       " 81: 3,\n",
       " 125: 5,\n",
       " 128: 2,\n",
       " 216: 6,\n",
       " 243: 3,\n",
       " 256: 4,\n",
       " 625: 5,\n",
       " 1024: 4,\n",
       " 1296: 6,\n",
       " 2187: 3,\n",
       " 3125: 5,\n",
       " 7776: 6,\n",
       " 16384: 4,\n",
       " 78125: 5,\n",
       " 279936: 6}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Czroots = make_Czroots(Apowers)\n",
    "Czroots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Czroots[243]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We see that `243` is in `Czroots`, with value `3`, so this would be a counterexample (except for the common factor).  The program uses the `exponent` function to recover the values of  `x, y, z`, and prints the results.\n",
    "\n",
    "# Is the Program Correct?\n",
    "\n",
    "Can we gain confidence in the program?  It is difficult to test `beal`, because the expected output is nothing, for all known inputs.\n",
    "One thing we can do is verify that `beal` finds cases like `3 ** 3 + 6 ** 3 == 3 ** 5 == 243` that would be a counterexample except for the common factor `3`. We can test this by temporarily replacing the `gcd` function with a mock function that always reports no common factors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "button": false,
    "collapsed": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3 ** 3 + 6 ** 3 == 3 ** 5 == 243\n",
      "7 ** 7 + 49 ** 3 == 98 ** 3 == 941192\n",
      "8 ** 4 + 16 ** 3 == 2 ** 13 == 8192\n",
      "8 ** 5 + 32 ** 3 == 16 ** 4 == 65536\n",
      "9 ** 3 + 18 ** 3 == 9 ** 4 == 6561\n",
      "16 ** 5 + 32 ** 4 == 8 ** 7 == 2097152\n",
      "17 ** 4 + 34 ** 4 == 17 ** 5 == 1419857\n",
      "19 ** 4 + 38 ** 3 == 57 ** 3 == 185193\n",
      "27 ** 3 + 54 ** 3 == 3 ** 11 == 177147\n",
      "28 ** 3 + 84 ** 3 == 28 ** 4 == 614656\n",
      "34 ** 5 + 51 ** 4 == 85 ** 4 == 52200625\n"
     ]
    }
   ],
   "source": [
    "def gcd(a, b): return 1\n",
    "\n",
    "beal(100, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let's make sure all those expressions are true:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "button": false,
    "collapsed": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{True}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "{3 ** 3 + 6 ** 3 == 3 ** 5 == 243,\n",
    " 7 ** 7 + 49 ** 3 == 98 ** 3 == 941192,\n",
    " 8 ** 4 + 16 ** 3 == 2 ** 13 == 8192,\n",
    " 8 ** 5 + 32 ** 3 == 16 ** 4 == 65536,\n",
    " 9 ** 3 + 18 ** 3 == 9 ** 4 == 6561,\n",
    " 16 ** 5 + 32 ** 4 == 8 ** 7 == 2097152,\n",
    " 17 ** 4 + 34 ** 4 == 17 ** 5 == 1419857,\n",
    " 19 ** 4 + 38 ** 3 == 57 ** 3 == 185193,\n",
    " 27 ** 3 + 54 ** 3 == 3 ** 11 == 177147,\n",
    " 28 ** 3 + 84 ** 3 == 28 ** 4 == 614656,\n",
    " 34 ** 5 + 51 ** 4 == 85 ** 4 == 52200625}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "I get nervous having an incorrect version of `gcd` around: change it back, quick!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "from math import gcd\n",
    "\n",
    "beal(100, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We can also provide some test cases for the subfunctions of `beal`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "button": false,
    "collapsed": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'tests pass'"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def tests():\n",
    "    assert make_Apowers(6, 10) == {\n",
    "         1: [1],\n",
    "         2: [8, 16, 32, 128],\n",
    "         3: [27, 81, 243, 2187],\n",
    "         4: [64, 256, 1024, 16384],\n",
    "         5: [125, 625, 3125, 78125],\n",
    "         6: [216, 1296, 7776, 279936]}\n",
    "    \n",
    "    assert make_Czroots(make_Apowers(5, 8)) == {\n",
    "        1: 1, 8: 2, 16: 2, 27: 3, 32: 2, 64: 4, 81: 3,\n",
    "        125: 5, 128: 2, 243: 3, 256: 4, 625: 5, 1024: 4,\n",
    "        2187: 3, 3125: 5, 16384: 4, 78125: 5}\n",
    "    Czroots = make_Czroots(make_Apowers(100, 100))\n",
    "    assert 3 ** 3 + 6 ** 3 in Czroots\n",
    "    assert 99 ** 97 in Czroots\n",
    "    assert 101 ** 100 not in Czroots\n",
    "    assert Czroots[99 ** 97] == 99\n",
    "    \n",
    "    assert exponent(10 ** 5, 10) == 5\n",
    "    assert exponent(7 ** 3, 7) == 3\n",
    "    assert exponent(1234 ** 999, 1234) == 999\n",
    "    assert exponent(12345 ** 6789, 12345) == 6789\n",
    "    assert exponent(3 ** 10000, 3) == 10000\n",
    "    assert exponent(1, 1) == 3\n",
    "    \n",
    "    assert exponents_upto(2) == []\n",
    "    assert exponents_upto(3) == [3]\n",
    "    assert exponents_upto(4) == [3, 4]\n",
    "    assert exponents_upto(40) == [3, 4, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]\n",
    "    assert exponents_upto(100) == [\n",
    "        3, 4, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, \n",
    "        67, 71, 73, 79, 83, 89, 97]\n",
    "    \n",
    "    assert gcd(3, 6) == 3\n",
    "    assert gcd(3, 7) == 1\n",
    "    assert gcd(861591083269373931, 94815872265407) == 97\n",
    "    assert gcd(2*3*5*(7**10)*(11**12), 3*(7**5)*(11**13)*17) == 3*(7**5)*(11**12)\n",
    "    \n",
    "    return 'tests pass'\n",
    "    \n",
    "tests()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "The program is mostly straightforward, but relies on the correctness of these arguments: \n",
    "\n",
    "* Are we justified in taking `combinations` without replacements from the `Apowers` table? In other words, are we sure there are no solutions of the form $A^x + A^x = C^z$? Yes, we can be sure, because then  $2\\;A^x = C^z$, and all the factors of $A$ would also be factors of $C$.\n",
    "\n",
    "<p>\n",
    "* Are we justified in having a single value for each key in the `Czroots` table? Consider that $81 = 3^4 = 9^2$. We put `{81: 3}` in the table and discard `{81: 9}`, because any number that has 9 as a factor will always have 3 as a factor as well, so 3 is all we need to know. But what if a number could be formed with two bases where neither was a multiple of the other? For example, what if $2^7 = 5^3 = s$; then wouldn't we have to have both 2 and 5 as values for $s$ in the table? Fortunately, that can never happen, because of the [fundamental theorem of arithmetic](https://en.wikipedia.org/wiki/Fundamental_theorem_of_arithmetic).\n",
    "\n",
    "<p>\n",
    "* Could there be a rounding error involving the `exponent` function that was not caught by the tests? Possibly; but `exponent` is not used to find counterexamples, only to print them, so any such error wouldn't cause us to miss a counterexample.\n",
    "\n",
    "<p>\n",
    "* Are we justified in only considering exponents that are odd primes, or the number 4? In one sense, yes, because when we consider the two terms $A^{(qp)}$ and $(A^q)^p$, we find they are always equal, and always have the same prime factors (the factors of $A$), so for the purposes of the Beal problem, they are equivalent, and we only need consider one of them. In another sense, there is a difference. With this optimization, when we run `beal(6, 10)`, we are no longer testing $512$ as a value of $A$ or $B$, even though $512 = 2^9$ and both $2$ and $9$ are within range, because the program chooses to express $512$ as $8^3$, and $8$ is not in the specified range. So the program is still correctly searching for counterexamples, but the space that it searches for given `max_A` and `max_x` is different with this optimization.\n",
    "\n",
    "<p>\n",
    "* Are we really sure that when $A$ and $B$ have a common factor greater than 1, then\n",
    "$C$ also shares that common factor?  Yes, because if $p$ is a factor of both $A$ and $B$, then it is a factor of $A^x + B^y$, and since we know this is equal to $C^z$, then $p$ must also be a factor of $C^z$, and thus a factor of $C$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "\n",
    "# Faster Arithmetic (mod *p*)\n",
    "\n",
    "Arithmetic is slow with integers that have thousands of digits. If we want to explore much further, we'll have to make the program more efficient. An obvious improvement would be to do all the arithmetic modulo some number $m$. Then we know:\n",
    "\n",
    "$$\\mbox{if} ~~ \n",
    "A^x + B^y = C^z\n",
    "~~ \\mbox{then} ~~\n",
    "(A^x (\\mbox{mod} ~ m) + B^y (\\mbox{mod} ~ m)) (\\mbox{mod} ~ m) = C^z \\;(\\mbox{mod} ~ m)$$\n",
    "\n",
    "\n",
    "So we can do efficient tests modulo $m$, and then do the full arithmetic only for combinations that work modulo $m$. Unfortunately there will be collisions (two numbers that are distinct, but are equal mod $m$), so the tables will have to have lists of values. Here is a simple, unoptimized implementation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "from math        import gcd\n",
    "from itertools   import combinations, product\n",
    "from collections import defaultdict\n",
    "                            \n",
    "def beal_modm(max_A, max_x, m=2**31-1):\n",
    "    \"\"\"See if any A ** x + B ** y equals some C ** z (mod p), with gcd(A, B) == 1.\n",
    "    If so, verify that the equation works without the (mod m).\n",
    "    Consider any 1 <= A,B <= max_A and x,y <= max_x, with x,y prime or 4.\"\"\"\n",
    "    assert m >= max_A\n",
    "    Apowers = make_Apowers_modm(max_A, max_x, m)\n",
    "    Czroots = make_Czroots_modm(Apowers)\n",
    "    for (A, B) in combinations(Apowers, 2):\n",
    "        if gcd(A, B) == 1:\n",
    "            for (Axm, x), (Bym, y) in product(Apowers[A], Apowers[B]):  \n",
    "                Czm = (Axm + Bym) % m\n",
    "                if Czm in Czroots:\n",
    "                    lhs = A ** x + B ** y\n",
    "                    for (C, z) in Czroots[Czm]:\n",
    "                        if lhs == C ** z:\n",
    "                            print('{} ** {} + {} ** {} == {} ** {} == {}'\n",
    "                                  .format(A, x, B, y, C, z, C ** z))                        \n",
    "                    \n",
    "\n",
    "def make_Apowers_modm(max_A, max_x, m): \n",
    "    \"A dict of {A: [(A**3 (mod m), 3), (A**4 (mod m), 4), ...]}.\"\n",
    "    exponents = exponents_upto(max_x)\n",
    "    return {A: [(pow(A, x, m), x) for x in (exponents if (A != 1) else [3])]\n",
    "            for A in range(1, max_A+1)}\n",
    "\n",
    "def make_Czroots_modm(Apowers): \n",
    "    \"A dict of {C**z (mod m): [(C, z),...]}\"\n",
    "    Czroots = defaultdict(list)\n",
    "    for A in Apowers:\n",
    "        for (Axm, x) in Apowers[A]:\n",
    "            Czroots[Axm].append((A, x))\n",
    "    return Czroots "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Here we see that each entry in the `Apowers` table is a list of `(A**x (mod p), x)` pairs.\n",
    "For example, $6^7 = 279,936$, so in our (mod 1000) table we have the pair `(936, 7)` under `6`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{1: [(1, 3)],\n",
       " 2: [(8, 3), (16, 4), (32, 5), (128, 7)],\n",
       " 3: [(27, 3), (81, 4), (243, 5), (187, 7)],\n",
       " 4: [(64, 3), (256, 4), (24, 5), (384, 7)],\n",
       " 5: [(125, 3), (625, 4), (125, 5), (125, 7)],\n",
       " 6: [(216, 3), (296, 4), (776, 5), (936, 7)]}"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Apowers = make_Apowers_modm(6, 10, 1000)\n",
    "Apowers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "And each item in the `Czroots` table is of the form `{C**z (mod m): [(C, z), ...]}`.\n",
    "For example, `936: [(6, 7)]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "button": false,
    "collapsed": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "defaultdict(list,\n",
       "            {1: [(1, 3)],\n",
       "             8: [(2, 3)],\n",
       "             16: [(2, 4)],\n",
       "             24: [(4, 5)],\n",
       "             27: [(3, 3)],\n",
       "             32: [(2, 5)],\n",
       "             64: [(4, 3)],\n",
       "             81: [(3, 4)],\n",
       "             125: [(5, 3), (5, 5), (5, 7)],\n",
       "             128: [(2, 7)],\n",
       "             187: [(3, 7)],\n",
       "             216: [(6, 3)],\n",
       "             243: [(3, 5)],\n",
       "             256: [(4, 4)],\n",
       "             296: [(6, 4)],\n",
       "             384: [(4, 7)],\n",
       "             625: [(5, 4)],\n",
       "             776: [(6, 5)],\n",
       "             936: [(6, 7)]})"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "make_Czroots_modm(Apowers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Let's run the program:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "button": false,
    "collapsed": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 56 s, sys: 436 ms, total: 56.4 s\n",
      "Wall time: 59.2 s\n"
     ]
    }
   ],
   "source": [
    "%time beal_modm(1000, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We don't see a speedup here, but the idea is that as we start dealing with much larger integers, this version should be faster. I could improve this version by caching certain computations, managing the memory layout better, moving some computations out of loops, considering using multiple different numbers as the modulus (as in a Bloom filter), finding a way to parallelize the program, and re-coding in a faster compiled language (such as C++ or Go or Julia). Then I could invest thousands (or millions) of CPU hours searching for counterexamples. \n",
    "\n",
    "But [Witold Jarnicki](https://plus.sandbox.google.com/+WitoldJarnicki/posts) and [David Konerding](http://www.konerding.com/~dek/) already did that: they wrote a C++ program that, in parallel across thousands of machines, searched for $A, B$ up to 200,000 and $x, y$ up to 5,000, but found no counterexamples.  So I don't think it is worthwhile to continue on that path.\n",
    "\n",
    "# Conclusion\n",
    "\n",
    "This was fun, but I can't recommend anyone spend a serious amount of computer time looking for counterexamples to the Beal Conjecture&mdash;the money you would have to spend in computer time would be more than the expected value of your prize winnings.  I suggest you work on a proof rather than a counterexample, or work on some other interesting problem instead!"
   ]
  }
 ],
 "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": 1
}