{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cab3c13d",
   "metadata": {},
   "source": [
    "[Home](Home.ipynb) <br/>\n",
    "[Crypto](Crypto.ipynb)\n",
    "\n",
    "# Fun with Fermat \n",
    "\n",
    "![Pierre de Fermat](https://ingrangel.net/wp-content/uploads/2019/01/image-338.jpg)\n",
    "\n",
    "Fermat is best known for [his famous last theorem](https://www.nsf.gov/discoveries/disc_summ.jsp?cntn_id=100029&org=NSF), which has been regarded as proved only recently, in 1994 by Andrew Wiles and Richard Taylor, meaning Fermat's proof, if he ever had one, is lost to history.  \n",
    "\n",
    "In light of the complexity of the existing proof, most are skeptical that he ever had one. He claimed he did, but it wouldn't fit in the margin. \"Yeah right\" right?\n",
    "\n",
    "### Fermat's Little Theorem\n",
    "\n",
    "However his so-called Little Theorem is likewise easy to state and, unlike the Last Theorem, [relatively easy to prove](YoutubeGallery.ipynb#Number-Theory).\n",
    "\n",
    "Statement:\n",
    "\n",
    "If $p$ is a prime and $a$ is any integer not divisible by p, then $a^{p − 1} − 1$ is divisible by p.\n",
    "\n",
    "Put another way:\n",
    "\n",
    "$a^{p − 1} \\equiv 1\\mod{p}$, where $a \\nmid p$.\n",
    "\n",
    "Yet another way:\n",
    "\n",
    "$(a^{p} - a) \\mid p$\n",
    "\n",
    "As long as $a$ and $p$ are strangers, and $p$ is prime, raising $a$ to the $p - 1$ power begets a number divisible by $p$ but for a remainder of 1.\n",
    "\n",
    "Where do we get all these fun little math symbols?  From [$\\LaTeX$](https://oeis.org/wiki/List_of_LaTeX_mathematical_symbols)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "135be4f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import primes\n",
    "from math import gcd\n",
    "\n",
    "p = primes.primes_gen.PrimeNumbers()\n",
    "some_primes = [next(p) for _ in range(20)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6efddeab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "some_primes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5ad187c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def strangers(a, b):\n",
    "    \"\"\"returns True when a, b have nothing in common\"\"\"\n",
    "    return gcd(a, b) == 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa1dc6c3",
   "metadata": {},
   "source": [
    "$p$ divides $a^{p-1} -1$ with no remainder, hence all these zeroes:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "096fdaef",
   "metadata": {},
   "source": [
    "The code cell below continues with list comprehension syntax, by adding an optional if clause, a filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "def488bd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = 17\n",
    "[(pow(a, p-1) - 1) % p for p in some_primes if strangers(a, p)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "032b73d1",
   "metadata": {},
   "source": [
    "Computing $a^{p-1} \\mod{p}$ directly is faster, using ```pow``` (which takes a 3rd argument, the modulus) and nets us that residue of one, every time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e24597f7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[pow(a, p-1, p) for p in some_primes]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "085219ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[pow(a, p-1, p) for p in some_primes if strangers(a,p)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c1aa8ac",
   "metadata": {},
   "source": [
    "### Fermat's Primality Test \n",
    "\n",
    "Given all of the above, i.e. the complete reliability of this theorem for primes (no proof yet), one might think plugging in arbitrary values for p and keeping only those for which the \"Fermat Test\" holds true, would be a recipie for getting primes.\n",
    "\n",
    "Lets get all the odd numbers between 3 and 101 and see which, according to the Fermat Test, might be prime."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8f5383d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 65, 67, 69, 71, 73, 75, 77, 79, 81, 83, 85, 87, 89, 91, 93, 95, 97, 99, 101]\n"
     ]
    }
   ],
   "source": [
    "candidates = [2] + [c for c in range(3,102,2)]\n",
    "print(candidates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5b39ae70",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101]\n"
     ]
    }
   ],
   "source": [
    "def fermat_test(c):\n",
    "    bases = [a for a in range(2,c) if strangers(a, c)]\n",
    "    for a in bases:\n",
    "        if not pow(a, c-1, c) == 1:\n",
    "            return False\n",
    "    return True\n",
    "\n",
    "passed = list(filter(fermat_test, candidates))\n",
    "print(passed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "57dc0767",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True]\n"
     ]
    }
   ],
   "source": [
    "print([primes.isprime(c) for c in passed])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85dd03f9",
   "metadata": {},
   "source": [
    "OMG, we got every one of these right!  Does that mean the Fermat test is always reliable?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7b67287",
   "metadata": {},
   "source": [
    "### Carmichael Numbers\n",
    "\n",
    "Enter the Carmichael Numbers..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ea433065",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[501, 503, 507, 509, 511, 513, 517, 519, 521, 523, 527, 529, 531, 533, 537, 539, 541, 543, 547, 549, 551, 553, 557, 559, 561, 563, 567, 569, 571, 573, 577, 579, 581, 583, 587, 589, 591, 593, 597, 599, 601]\n"
     ]
    }
   ],
   "source": [
    "candidates = [c for c in range(501,602,2) if strangers(c, 5)]\n",
    "print(candidates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c1c4a114",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[503, 509, 521, 523, 541, 547, 557, 561, 563, 569, 571, 577, 587, 593, 599, 601]\n"
     ]
    }
   ],
   "source": [
    "passed = list(filter(fermat_test, candidates))\n",
    "print(passed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d1b29fcb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[True, True, True, True, True, True, True, False, True, True, True, True, True, True, True, True]\n"
     ]
    }
   ],
   "source": [
    "print([primes.isprime(c) for c in passed])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56501a0d",
   "metadata": {},
   "source": [
    "Oh no, we got a False positive.  The number in that location is not actually prime. What is it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "06e2b107",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "561"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "passed[7]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "68834c9c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 3, 11, 17)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "primes.factors(561)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9ec3262",
   "metadata": {},
   "source": [
    "561 turns out to be the first of a special set of positive integers known as Carmichael Numbers and/or Fermat Pseudo Primes.  No matter the base, these numbers always pass the Fermat Test.\n",
    "\n",
    "Does this mean that some composites might pass the Fermat Test for one base, but not another, always preserving the stipulation that the base and the candidate are strangers?  Investigate.  Searching the internet is fine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "32e6afec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def fermat_test(c):\n",
    "    global liar, witness, gcd_witness  # postmortem \n",
    "    liar, witness, gcd_witness = [], None, None\n",
    "    for a in [a for a in range(2,c)]:\n",
    "        if not strangers(a, c):\n",
    "            gcd_witness = a\n",
    "            return False\n",
    "        if not pow(a, c-1, c) == 1:\n",
    "            witness = a  # Fermat witness\n",
    "            return False\n",
    "        else: \n",
    "            liar.append(a) # not liar if c is prime\n",
    "    return True\n",
    "\n",
    "fermat_test(341)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "506ff6bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "gcd_witness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "4d4cbe5b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "witness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5a1b03f7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "liar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "f8c44751",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fermat_test(541)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "c6c90f3e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "539"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(liar)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e845b772",
   "metadata": {},
   "source": [
    "Questions / Activities:\n",
    "\n",
    "* Do some composites pass the Fermat Test that are not Carmichael Numbers?\n",
    "* What are some of the Carmichael Numbers after 561?  \n",
    "* Explore between 1000 and 2000.\n",
    "* Check OEIS for [a listing](http://oeis.org/A002997).\n",
    "* Reading: [Teenager Solves Stubborn Riddle About Prime Number Look-Alikes](https://www.quantamagazine.org/teenager-solves-stubborn-riddle-about-prime-number-look-alikes-20221013/)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}