Mills’ constant is a number θ such that the integer part of θ raised to a power of 3 is always a prime. We’ll see if we can verify this computationally with SymPy.

from sympy import floor, isprime from sympy.mpmath import mp, mpf # set working precision to 200 decimal places mp.dps = 200 # Mills' constant http://oeis.org/A051021 theta = mpf("1.30637788386308069046861449260260571291678458515671364436805375996643405376682659882150140370119739570729") for i in range(1, 7): n = floor(theta**3**i) print i, n, isprime(n)

Note that the line of code defining `theta`

wraps Mill’s constant as a string and passes it as an argument to `mpf`

. If we just wrote `theta = 1.30637788386308…`

then `theta`

would be truncated to an ordinary floating point number and most of the digits would be lost. Not that I did that in the first version of the code.

This code says that the first five outputs are prime but that the sixth is not. That’s because we are not using Mills’ constant to enough precision. The integer part of θ^3^6 is 85 digits long, and we don’t have enough precision to compute the last digit correctly.

If we use the value of Mills’ constant given here to 640 decimal places, and increase `mp.dps`

to 1000, we can correctly compute the sixth term, but not the seventh.

Mill’s formula is just a curiosity with no practical use that I’m aware of, except possibly as an illustration of impractical formulas.