Published August 20, 2020|3 minute read
Inspired by Matt Parker’s recent video, I decided to search for primes with , the first of which is the 46-digit
How do you go about finding more?
Well, we want to be big — very big. From high-school trigonometry, we know that this occurs when is just a tiny bit less than a half-integer multiple of . In other words, we want
which rearranges to
Moreover, for small , we have the excellent approximation . Requiring then imposes the condition , or
Summing up the above, we want rational approximations to that are a bit too small, but exceedingly close.[1]How do we go about looking for such approximations? The continued fraction expansion, of course!
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def continued_fraction_coefficients():
"""Read the continued fraction coefficients of pi from a file."""
# Left as an exercise to the reader
def rational_approximations():
"""Generate the sequence of best rational approximants of pi."""
# We use the isomorphism between SL_2(R) and the fractional linear transformations.
# We start with the function 1/x, and repeatedly pre-compose
# with 1/(a + x), where a is the continued fraction coefficient.
# We return the limit of the current rational function as x → ∞.
n0, n1 = 0, 1
d0, d1 = 1, 0
for a in continued_fraction_coefficients():
(n0, n1) = (n1, a * n1 + n0)
(d0, d1) = (d1, a * d1 + d0)
yield (n1, d1)
For each rational approximation with even and odd, we’d like to check if . But since we’re going to be working with huge numbers, we want to avoid floating point math if possible.
One way to accomplish this is to replace the in this check with an even better rational approximation . For this check to be accurate, we’d like . The condition then becomes
which only involves integer math. Summing up the above:
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
gen1 = rational_approximations()
gen2 = rational_approximations()
for n, d in gen1:
# We want n/d with n even, d odd.
if n % 2 != 0 or d % 2 == 0: continue
p = n // 2
# We'd like d' >> pd; I found that d' > 4pd works fine.
threshold = 4 * p * d
while True:
nn, dd = next(gen2)
if dd > threshold:
break
# Our integer-math-only check for an 'interesting' value of p.
if 0 <= p * (d*nn - 2*p*dd) < 2*dd:
# At this point, we should have tan(p) > p.
# All that remains is a primality check:
if is_probable_prime(p):
print(p)
If you’re following along, you might want to add more logging just to see what the code is doing. All that remains is to define a reasonable probabilistic primality test is_probable_prime
. Personally, I got the best performance by doing trial division on primes up to a million, then calling PFGW via subprocess
. If you want to stay within Python, the gmpy2 library is also a reasonable option:
1
2
3
4
5
6
7
from gmpy2 import is_strong_prp, is_prime
PRIMES = [k for k in range(1, 1000000) if is_prime(k)]
def is_probable_prime(n):
if any(n % p == 0 for p in PRIMES):
return False
return is_strong_prp(n, 3)
The verdict? In less than a second, we reproduce the 46-digit prime mentioned at the start of the blog post. In two or three seconds, we find a 1017-digit prime solution (which Matt mentioned), which begins and ends as
And after a few hours, we find the monstrous 35,085-digit prime
which, at the time of writing, I believe Matt (and, perhaps, everybody else in the world) is unaware of.
I think that's enough for today.