Published August 20, 20203 minute read

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

p=1169809367327212570704813632106852886389036911.p = 1169809367327212570704813632106852886389036911.

p=116980936732886389036911.p = 116980936732 \dots 886389036911.

How do you go about finding more?

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

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

which rearranges to

π=2p2k+1+2ϵ2k+1.\pi = \frac{2p}{2k + 1} + \frac{2\epsilon}{2k+1}.

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

p[(2k+1)π2p][0,2).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.[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/dn/d with n=2pn = 2p even and d=2k+1d = 2k+1 odd, we’d like to check if p(dπ2p)[0,2)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/dn' / d'. For this check to be accurate, we’d like dpdd' \gg p \cdot d. The condition then becomes

0p(dn2pd)<2d,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=230835870782610326005069.p = 230835870782\dots 610326005069.

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

p=409461998988829370185991,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. [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).