{ "cells": [ { "cell_type": "markdown", "id": "87fa1305", "metadata": {}, "source": [ "# 100001st Prime\n", "> By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.\n", ">\n", "> What is the 10001st prime number?\n", "\n", "Unlike integer factorization, we *can* efficiently determine whether a number is prime. One of the simplest and fastest methods is the [Miller-Rabin primality test](https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test), which is a probabilistic test.\n", "\n", "Of course, SageMath has tools to let us compute the 10001st prime without implementing a primality test ourselves! " ] }, { "cell_type": "code", "execution_count": 1, "id": "3a27f454", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "104743" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Primes.unrank is 0-indexed, so the 10001st prime is at index 10000\n", "Primes().unrank(10000)" ] }, { "cell_type": "markdown", "id": "6f644292", "metadata": {}, "source": [ "Let's look into how the Miller-Rabin test works anyway, though.\n", "\n", "## Fermat's little theorem\n", "First, a little bit of elementary number theory. [Fermat's little theorem](https://en.wikipedia.org/wiki/Fermat%27s_little_theorem) says that if $p$ is a prime number, then for any integer $a$,\n", "$$a^p \\equiv a \\pmod{p}$$\n", "It follows from the [contrapositive](https://en.wikipedia.org/wiki/Contraposition) that if there exists an integer $a$ such that $a^p \\not\\equiv a \\pmod{p}$, then $p$ is not prime. We then call $a$ a *witness* that proves $p$ is composite.\n", "\n", "Consider the fact that $3^{6} \\equiv 3 \\pmod{6}$, but 6 is composite. In cases like this, we say that $p$ is a *base $a$ Fermat pseudoprime*, e.g. 6 is a base 3 Fermat pseudoprime. 6 is also a base 4 Fermat pseudoprime, but 2 and 5 both serve as witnesses to 6 being composite (we need only consider bases modulo $p$).\n", "\n", "This may lead us to believe that we can just check every base $a$, and if there are no witnesses, then $p$ is prime. However, (along with this method being inefficient for large $p$), there are numbers that are Fermat pseudoprimes for *every* base, called [Carmichael numbers](https://en.wikipedia.org/wiki/Carmichael_number). Therefore, the [converse](https://en.wikipedia.org/wiki/Converse_(logic)) of Fermat's little theorem does not hold.\n", "\n", "Nevertheless, Carmichael numbers are relatively rare, so if you pick a random number to test and a few random bases, and find no witnesses, you can be pretty sure the number is in fact prime.\n", "\n", "## Miller-Rabin test\n", "\n", "The Miller-Rabin test employs another useful fact about prime numbers: if $p$ is prime, then the only values of $x$ such that \n", "$$x^2 \\equiv 1 \\pmod{p}$$\n", "are $x \\equiv \\pm 1 \\pmod{p}$. (*Proof: if $x^2 \\equiv 1 \\pmod{p}$, then $p \\vert (x+1)(x-1)$. Then by [Euclid's lemma](https://en.wikipedia.org/wiki/Euclid%27s_lemma), $p$ must divide one of $x+1$ or $x-1$, implying $x \\equiv \\pm 1 \\pmod{p}$.*)\n", "\n", "Similarly to Fermat's little theorem, the converse of this statement doesn't hold. For example, $x^2 \\equiv 1 \\pmod{9}$ only has $x=1$ and $x=-1$ as solutions, but 9 is composite.\n", "\n", "We now have two questions to lead an investigation into the primality of a number $n$:\n", "* Is there an integer $a$ such that $a^n \\not\\equiv a \\pmod{n}$?\n", "* Does 1 have a non-trivial square root modulo $n$?\n", "\n", "If the answer to either of these is yes, then $n$ is *definitely* composite. If both answers are no, this is strong evidence that $n$ is prime, *but not proof*, since we have already demonstrated that both of these conditions can have false positives. A composite number that is erroneously suggested to be prime by these conditions is called a *base $a$ strong pseudoprime*. Fortunately, no number (including Carmichael numbers) is a base $a$ strong pseudoprime for all bases $a$. In fact, Rabin (one of the namesakes of the algorithm) proved that if $n$ is composite, at least 3/4 of the numbers 1 to $n-1$ are witnesses.\n", "\n", "Let's finally look at the actual test. Suppose we want to test if an odd number $n$ is prime (testing whether an even number is prime is trivial). There must exist $t \\geq 1$ and odd $u$ such that $n - 1 = 2^t u$." ] }, { "cell_type": "code", "execution_count": 2, "id": "720f7e8c", "metadata": {}, "outputs": [], "source": [ "def factor_twos(n):\n", " t = 0\n", " u = n - 1\n", " while u % 2 == 0:\n", " u //= 2\n", " t += 1\n", " \n", " return t, u" ] }, { "cell_type": "markdown", "id": "980bd2cf", "metadata": {}, "source": [ "It follows that\n", "$$a^{n-1} \\equiv a^{2^t u} \\pmod{n}$$\n", "for any base $a$, so we'll randomly select a base $1 \\leq a \\leq n-1$.\n", "\n", "The test proceeds by repeatedly squaring $a^u$ modulo $n$ until we reach $a^{n-1}$. If at any point in this squaring, we find that $a^{2^{i-1} u}$ does not equal 1 or -1, but the following squaring gives $a^{2^i u} = 1$, then we have found a non-trivial square root of 1 modulo $n$. Therefore, $n$ must be composite.\n", "\n", "After $t$ squarings, we will reach $a^{n-1}$. If this equals anything other than 1, then by the contrapositive of Fermat's last theorem, this also tells us that $n$ is composite.\n", "\n", "If, after $t$ squarings, we have not found a non-trivial square root of 1 modulo $n$ and $a^{n-1} \\equiv 1 \\pmod{n}$, then either $n$ is prime, or is a base $a$ strong pseudoprime." ] }, { "cell_type": "code", "execution_count": 3, "id": "123d1068", "metadata": {}, "outputs": [], "source": [ "def witness(a, n):\n", " t, u = factor_twos(n)\n", " x = pow(a, u, n)\n", " for _ in range(0, t):\n", " x, y = pow(x, 2, n), x\n", " if x == 1 and y != 1 and y != n - 1:\n", " return True\n", " \n", " return x != 1" ] }, { "cell_type": "markdown", "id": "f8454603", "metadata": {}, "source": [ "You can reduce the chance of a false report of primality by repeating the procedure with different bases. If we test several bases and do not find any witnesses, it is exceedingly unlikely that $n$ is composite." ] }, { "cell_type": "code", "execution_count": 4, "id": "e71ba6ec", "metadata": {}, "outputs": [], "source": [ "from random import randint\n", "\n", "def is_prime(n, k=10):\n", " if n < 2:\n", " return False\n", " elif n == 2:\n", " return True\n", " elif n % 2 == 0:\n", " return False\n", " \n", " for _ in range(0, k):\n", " a = randint(1, n - 1)\n", " if witness(a, n):\n", " return False\n", " \n", " return True" ] } ], "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 }