[Python] Optimizing The Exact Prime Number Theorem

For example, given this sequence of the first 499 primes, can you predict the next prime?

2,3,5,7,...,3541,3547,3557,3559


The 500th prime is 3571.

Prime Number Theorem


The Prime Number Theorem (PNT) provides an approximation for the n-th prime:



Computing p_500 ≈ 3107 takes microseconds!

Exact Prime Number Theorem


My experimental Exact Prime Number Theorem (EPNT) computes the exact n-th prime:



Computing p_500 = 3571 takes 25 minutes!

Question


So far, the EPNT correctly predicts the first 500 primes.

Unfortunately, numerically verifying the formula for higher primes is extremely slow!

Are there any optimization tips to improve the EPNT computational speed? Perhaps

• Do not use Python
• Add multiple threads
• Implement a faster math precision library
• Modify the decimal precision mp.dps at runtime
• Use a math computing engine like WolframAlpha

Here's the current Python code:

import time
from mpmath import ceil, ln, mp, mpf, exp, fsum, power, zeta
from sympy import symbols, Eq, pprint, prime

N=500 # <--- Compute the N-th prime.

mp.dps = 20000
primes = []

def vengy_prime(k):
# Compute the k-th prime deterministically
s = ceil(k * ln(k * ln(k)))

# Determine the dynamic Rosser (1941) upper bound
N = int(ceil(k * (ln(k) + ln(ln(k)))))

# Compute finite summation to N
print(f"Computing {N} zeta terms ...")
start_time = time.time()
sum_N = fsum([1 / power(mpf(n), s) for n in range(1, N)])
end_time = time.time()
print(f"Time taken: {end_time - start_time:.6f} seconds")

# Compute the product term involving the first k-1 primes
print(f"Computing product of {k-1} previous primes ...")
start_time = time.time()
prod = exp(fsum([ln(1 - power(p, -s)) for p in primes[:k-1]]))
end_time = time.time()
print(f"Time taken: {end_time - start_time:.6f} seconds")

# Compute next prime p_k
p_k=ceil((1 - 1 / (sum_N * prod)) ** (-1 / s))
return p_k


# Generate the previous known k-1 primes
print("\nListing", N-1, "known primes:")
for k in range(1, N):
p = prime(k)
primes.append(p)
print(primes)

primes.append(vengy_prime(N))
pprint(Eq(symbols(f'p_{N}'), int(primes[-1])))

Update


Wow! Running Jérôme Richard's new optimized code only took 10 seconds!

Computing 4021 zeta terms ...
Time taken: 7.968423 seconds
Computing product of 499 previous primes ...
Time taken: 1.960771 seconds
p₅₀₀ = 3571


The old code timings were 1486 seconds:

Computing 4021 zeta terms ...
Time taken: 1173.899538 seconds
Computing product of 499 previous primes ...
Time taken: 313.833039 seconds
p₅₀₀ = 3571


The optimized code computed the 4000th prime in 45 minutes: N = 4000, precision = 700000



Continue reading...