# A prime-generating formula and SymPy

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.

## 2 thoughts on “A prime-generating formula and SymPy”

1. Maybe I’m thinking of another thing, but isn’t this true for all i only if the Riemann Hypothesis is true?

2. I don’t think that Mill’s result depends on RH. However, the value of Mill’s constant, or at least the algorithm used to compute it, does depend on RH.