Published August 20, 20203 minute read

Inspired by Matt Parker’s recent video, I decided to search for primes $p$ with $\tan{p} > p$, the first of which is the 46-digit

$p = 1169809367327212570704813632106852886389036911.$

$p = 116980936732 \dots 886389036911.$

How do you go about finding more?

Well, we want $\tan{p}$ to be big — very big. From high-school trigonometry, we know that this occurs when $p$ is just a tiny bit less than a half-integer multiple of $\pi$. In other words, we want

$p + \epsilon = \left(k + \frac{1}{2}\right) \pi,$

which rearranges to

$\pi = \frac{2p}{2k + 1} + \frac{2\epsilon}{2k+1}.$

Moreover, for small $\epsilon$, we have the excellent approximation $\tan{p} \sim \epsilon^{-1}$. Requiring $\tan{p} > p$ then imposes the condition $p \epsilon \in [0, 1)$, or

$p \cdot \left[(2k + 1) \, \pi - 2p\right] \in [0, 2).$

Summing up the above, we want rational approximations to $\pi$ that are a bit too small, but exceedingly close.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.
# 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 $n/d$ with $n = 2p$ even and $d = 2k+1$ odd, we’d like to check if $p \cdot (d \pi - 2p) \in [0, 2)$. 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 $\pi$ in this check with an even better rational approximation $n' / d'$. For this check to be accurate, we’d like $d' \gg p \cdot d$. The condition then becomes

$0 \leq p \cdot (dn' - 2pd') < 2d',$

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

$p = 230835870782\dots 610326005069.$

And after a few hours, we find the monstrous 35,085-digit prime

$p = 409461998988\dots 829370185991,$

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.

## Footnotes

1. More precisely, we want an approximation accurate to the order of denominator squared. The constants are a bit too small to apply Dirichlet’s approximation theorem directly, so we aren’t guaranteed infinitely many solutions (although I suspect there are).