Primes $p$ with $\tan{p} > p$

Published August 20, 2020｜3 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

which rearranges to

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

Summing up the above, we want rational approximations to $\pi$ 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 $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

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.