100001st Prime¶
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, which is a probabilistic test.
Of course, SageMath has tools to let us compute the 10001st prime without implementing a primality test ourselves!
nth_prime(10001)
104743
Let's look into how the Miller-Rabin test works anyway, though.
Fermat's little theorem¶
First, a little bit of elementary number theory. Fermat's little theorem says that if $p$ is a prime number, then for any integer $a$, $$a^p \equiv a \pmod{p}$$ It follows from the contrapositive 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.
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$).
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. Therefore, the converse of Fermat's little theorem does not hold.
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.
Miller-Rabin test¶
The Miller-Rabin test employs another useful fact about prime numbers: if $p$ is prime, then the only values of $x$ such that $$x^2 \equiv 1 \pmod{p}$$ 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, $p$ must divide one of $x+1$ or $x-1$, implying $x \equiv \pm 1 \pmod{p}$.)
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.
We now have two questions to lead an investigation into the primality of a number $n$:
- Is there an integer $a$ such that $a^n \not\equiv a \pmod{n}$?
- Does 1 have a non-trivial square root modulo $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 both of 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.
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$.
def factor_twos(n):
t = 0
u = n - 1
while u % 2 == 0:
u //= 2
t += 1
return t, u
It follows that $$a^{n-1} \equiv a^{2^t u} \pmod{n}$$ for any base $a$, so we'll randomly select a base $1 \leq a \leq n-1$.
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.
After $t$ squarings, we will reach $a^{n-1}$. If this equals anything other than 1, then by the contrapositive of Fermat's little theorem, this also tells us that $n$ is composite.
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.
def witness(a, n):
t, u = factor_twos(n)
x = pow(a, u, n)
for _ in range(0, t):
x, y = pow(x, 2, n), x
if x == 1 and y != 1 and y != n - 1:
return True
return x != 1
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.
from random import randint
def is_prime(n, k=10):
if n < 2:
return False
elif n == 2:
return True
elif n % 2 == 0:
return False
for _ in range(0, k):
a = randint(1, n - 1)
if witness(a, n):
return False
return True
With this implementation, we can generate the primes by iterating over the positive integers and testing if each is prime.
from itertools import count
def primes(start=1):
for n in count(start):
if is_prime(n):
yield n
Then to solve the problem, we can just count the primes until we reach 10001.
Rosser's theorem¶
This strategy works... but we can go even deeper.
Running Miller-Rabin on every integer until we reach the 10001st prime isn't super efficient. To avoid this, we can take advantage of Rosser's theorem, which states that $p_n > n \log n$. In other words, we can start our primality testing at $10001 \log 10001$ instead of 1.
lower = ceil(10001 * log(10001))
lower
92114
Prime-counting function¶
There's one issue preventing us from jumping ahead to 92114: we don't know how many primes are less than 92114, so when we find the next prime (spoiler: it's 92119), we won't know what number prime that is (another spoiler: it's the 8897th prime).
To answer this question, we can use the prime-counting function. SageMath has the prime-counting function available directly as prime_pi, but let's implement it ourselves. As a first pass, we can use a sieve. Check out problem 10 for a sieve implementation - to avoid duplicating work, we'll use SageMath's sieve here.
def prime_counting_function(n):
return len(prime_range(n+1))
This will tell us where to start our tally if we start checking for primes at 92114.
tally = prime_counting_function(lower)
tally
8896
This runs considerably quicker than starting from 1.
for (i, p) in enumerate(primes(start=lower), start=tally+1):
if i == 10001:
break
p
104743
Relevant sequences¶
Copyright (C) 2025 filifa¶
This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International license and the BSD Zero Clause license.