An explanation of a closed formula for the nth prime number and why it's not very interesting

closed_form_for_primes.py

# based on "closed form formula" for the nth prime, defined by
# "On Formulae for the nth Prime Number", C. P. Willans
# (https://www.jstor.org/stable/3611701)
# 
# which I discovered through a stackexchange question at
# https://hsm.stackexchange.com/questions/13353/who-discovered-this-closed-form-formula-for-the-n-th-prime-number
from math import ceil, cos, factorial as fac, floor, pi, sqrt

# The direct definition. This produces incorrect answers because of floating
# point inaccuracy, and raises OverflowError for n>7. Let's fix this.
def nth_prime_incorrect(n):
    return 1 + sum(floor(
        (n /
         sum(cos(cos((fac(i-1)+1) / i * pi))
             for i in range(1,k+1))
         ) ** (1/n)
    ) for k in range(1, 1 + 2 ** n))

# The crux of this formula is Wilson's theorem: ((x-1)!+1)/x is an integer iff x
# is prime or 1. Factorial grows so fast, however, that this division yields
# floating-point overflow when x > 172.
def wilson(x): return (fac(x-1) + 1)/x

# Mathematically, F(x) = ⌊cos²(π wilson(x))⌋ satisfies F(x)=1 if x is prime or
# 1, and F(x)=0 for x composite. Because of floating-point imprecision, however,
# we actually compute F_incorrect(17) = 0 because
#
#   cos(pi * wilson(17)) ** 2 = 0.9999999579555099
#
# So using this to implement the closed form formula produces wrong answers.
def F_incorrect(x): return floor(cos(pi * wilson(x)) ** 2)

# Therefore I use the simpler and more readily computable prime_or_one. It's
# still using factorial, so it's too expensive to compute for large x. A more
# efficient prime test, such as simply testing division by all numbers up to
# sqrt(x), or using Pollard-Rho, would improve this.
def prime_or_one(x): return (fac(x-1) + 1) % x == 0

# Since F(x) is 1 exactly where x is prime or 1, we can define
#
#   π(m) = -1 + ∑{x=1 to m} F(x)
#
# which computes the number of primes ≤m. I don't know why the paper doesn't
# just ditch the -1 and use ∑{x=2 to m} F(x).
def count_primes_upto(m):
    return sum(1 if prime_or_one(x) else 0 for x in range(2, 1+m))

# IMO this is why this "closed form formula" for primes is uninteresting. What
# counts as "closed form" is a fairly arbitrary definition, one not of
# fundamental theoretical interest but rather a historical accident. "Test
# divisibility up to sqrt(x)" is _simpler_ than F(x). But "closed form"
# formulae, by convention, are allowed to use numbers, arithmetic, and summation
# over a range, but not booleans, inequality tests, or conjunction over a range.

# The function A_n(a) is another example: A_n(a) = ⌊n√(n/(1+a))⌋ is more
# directly defined by A_n(a) = 1 if a < n else 0. What a stupid way to express
# such a simple concept!
def A(n,a): return floor((n / (1+a)) ** (1/n))

# Putting everything together, for sufficiently large N, eg N=2^m, we have
#
#   prime(n) = 1 + ∑{m=1 to N} A_n(π(m))
#            = 1 + ∑{m=1 to N} 1 if π(m) < n else 0
#            = |{m ∈ ℕ : π(m) < n}|               (NB. 0 ∈ ℕ, so we drop the 1+)
#            = |{m ∈ ℕ : there are fewer than n primes p≤m}|
#            = the number of naturals before the nth prime
#
# which explains both why this is correct and why it's a bad way to compute the
# nth prime: we test all numbers below 2^n for primality using a method that
# involves computing their factorial. There is so much wasted work here it's
# embarrassing.
def nth_prime(n):
    return 1 + sum(
        1 if count_primes_upto(m) < n else 0 # A_n(π(m))
        for m in range(1, 1 + 2 ** n))

# Here's a stupidly obvious, totally unoptimized, but still much faster way to
# find the first n primes: sweep through the naturals ≥2, collecting primes by
# testing if each number is divisible by any prime found so far.
def first_n_primes(n):
    primes = []
    i = 2
    while len(primes) < n:
        if all(i % p != 0 for p in primes): primes.append(i)
        i += 1
    return primes

# And here's a prime sieve, which is faster for large n, eg n=10_000.
def sieve(n):
    primes = [2,3]
    K = 2
    while True:
        # Check [K,2K) for compositeness. Note that any new primes p≥K cannot
        # have multiples in [K,2K) because 2p ≥ 2K.
        K += K
        composite = [False] * K
        for p in primes:
            # We don't need to check composites below p^2 because they must have
            # a prime factor q < p, which we've already checked.
            if K+K <= p * p: break
            next_multiple = K + (-K % p) # next_multiple % p == 0
            for x in range(next_multiple, K+K, p):
                composite[x-K] = True # x % p == 0
        # Now extract the new primes.
        primes.extend(K+i for i,c in enumerate(composite) if not c)
        if len(primes) >= n: return primes[:n]