what is the best approach to compute Brun's constant with python (witch include primality test and long run sum)?

Question

Brun's constant : https://en.wikipedia.org/wiki/Brun%27s_theorem http://numbers.computation.free.fr/Constants/Primes/twin.html

How to compute Brun's constant up to 10^20 with python knowing that primality check has a heavy cost and summing result up to 10^20 is a long task

here is my 2 cents attempt :

IsPrime : fastest way to check if the number is prime I know digit_root : compute the digital root of the number

If someone know what could be improve bellow to reach computation to 10^20, you're welcome

import numpy as np
import math
import time

#Brun's constant
#p      B_2(p) 
#10^2   1.330990365719...
#10^4   1.616893557432...
#10^6   1.710776930804...
#10^8   1.758815621067...
#10^10  1.787478502719...
#10^12  1.806592419175...
#10^14  1.820244968130...
#10^15  1.825706013240...
#10^16  1.830484424658...
#B_2 should reach 1.9 at p ~ 10^530 which is far beyond any computational project

#B_2*(p)=B_2(p)+ 4C_2/log(p)
#p  B2*(p) 
#10^2   1.904399633290...
#10^4   1.903598191217...
#10^6   1.901913353327...
#10^8   1.902167937960...
#10^10  1.902160356233...
#10^12  1.902160630437...
#10^14  1.902160577783...
#10^15  1.902160582249...
#10^16  1.902160583104...

def digit_root(number):
    return (number - 1) % 9 + 1

first25Primes=[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
def IsPrime(n) : 
    # Corner cases 
    if (n <= 1) : 
        return False
    if (n <= 3) : 
        return True
  
    # This is checked so that we can skip  
    # middle five numbers in below loop 
    if (n % 2 == 0 or n % 3 == 0) : 
        return False
    #exclude digital root 3 or 6 or 9
    if digit_root(n) in (3,6,9):
        return False
    
    if (n != 2 and n != 7 and n != 5 and str(n)[len(str(n))-1] not in ("1","3","7","9")): #si le nombre ne se termine pas par 1 3 7 9
        return False

    for i in first25Primes:
        if n%i == 0 and i < n:
            return False
            
    if (n>2):
        if (not(((n-1) / 4) % 1 == 0 or ((n+1) / 4) % 1 == 0)):
            return False
    if (n>3):
        if (not(((n-1) / 6) % 1 == 0 or ((n+1) / 6) % 1 == 0)):
            return False
    
    i = 5
    while(i * i <= n) :
        if (n % i == 0 or n % (i + 2) == 0) : 
            return False
        i = i + 6      
  
    return True

def ComputeB_2Aster(B_2val,p):
    return B_2 + (C_2mult4/np.log(p))


start = time.time()
#approx De Brun's
B_2 = np.float64(0)
B_2Aster = np.float64(0)
one = np.float64(1)
#Twin prime constant
C_2 = np.float64(0.6601618158468695739278121100145557784326233602847334133194484233354056423)
C_2mult4 = C_2 * np.float64(4)
lastPrime = 2
lastNonPrime = 1
for p in range(3, 100000000000000000000,2):    
    if IsPrime(p):
        lastNonPrime = p-1
        if lastPrime == p-2 and lastNonPrime == p-1:
            B_2 = B_2 + (one/np.float64(p)) + (one/np.float64(lastPrime))            
        lastPrime = p
    else:
        lastNonPrime = p
    
    if p<10000000000:
        if p%1000001==0:            
            print(f'p:{p} \t\t[elapsed:{time.time()-start}]\nB_2:{B_2:.52f} B_2Aster:{ComputeB_2Aster(B_2,p-2):.52f}\n',end="")
    else:        
        print(f'p:{p} \t\t[elapsed:{time.time()-start}]\nB_2:{B_2:.52f} B_2Aster:{ComputeB_2Aster(B_2,p-2):.52f}\n',end="")
print(f'p:{p} \t\t[elapsed:{time.time()-start}]\nB_2:{B_2:.52f} B_2Aster:{ComputeB_2Aster(B_2,p-2):.52f}\n',end="")

Answer 1

Now, if you remember your classes in number theory there is something called sieve of Eratosthenes which is an algorithm to find all prime numbers up to any given limit:

algorithm Sieve of Eratosthenes is
    input: an integer n > 1.
    output: all prime numbers from 2 through n.

    let A be an array of Boolean values, indexed by integers 2 to n,
    initially all set to true.
    
    for i = 2, 3, 4, ..., not exceeding √n do
        if A[i] is true
            for j = i2, i2+i, i2+2i, i2+3i, ..., not exceeding n do
                set A[j] := false

    return all i such that A[i] is true.

Now, this algorithm, translated into python is simply

def sieve_of_eratosthenes(max_num):
    is_prime = np.full(max_num + 1, True, dtype=bool)
    is_prime[0:2] = False
    p = 2
    while (p * p <= max_num):
        if (is_prime[p] == True):
            for i in range(p * p, max_num + 1, p):
                is_prime[i] = False
        p += 1
    return np.nonzero(is_prime)[0]

Now, to compute Brun's constant, you will need that list of primes and the sum. So,

def brun_constant(primes):
    sum_primes = mpfr('0')
    last_p = mpz('2')
    count = 0
    for p in primes:
        p_gmp = mpz(p)  
        if count > 0:
            if is_prime(p_gmp + 2):
                sum_primes += (1/mpfr(p_gmp) + 1/mpfr(p_gmp + 2))
        last_p = p_gmp
        count += 1
    return float(sum_primes)

So, your full code will be

import numpy as np
import gmpy2
from gmpy2 import mpz, mpfr, is_prime
import time

def sieve_of_eratosthenes(max_num):
    is_prime = np.full(max_num + 1, True, dtype=bool)
    is_prime[0:2] = False
    p = 2
    while (p * p <= max_num):
        if (is_prime[p] == True):
            for i in range(p * p, max_num + 1, p):
                is_prime[i] = False
        p += 1
    return np.nonzero(is_prime)[0]

def brun_constant(primes):
    sum_primes = mpfr('0')
    last_p = mpz('2')
    count = 0
    for p in primes:
        p_gmp = mpz(p)  
        if count > 0:
            if is_prime(p_gmp + 2):
                sum_primes += (1/mpfr(p_gmp) + 1/mpfr(p_gmp + 2))
        last_p = p_gmp
        count += 1
    return float(sum_primes)

start_time = time.time()

limit = 10**6
primes = sieve_of_eratosthenes(limit)
print(primes)
brun_const = brun_constant(primes)

end_time = time.time()
execution_time = end_time - start_time

print(f"Brun's constant up to {limit}: {brun_const}")
print(execution_time)

Note here I chose 10^6 as an example:

[     2      3      5 ... 999961 999979 999983]
Brun's constant up to 1000000: 1.710776930804221
0.4096040725708008

So, it took 0.4 seconds.

For 10^8 I got

[       2        3        5 ... 99999959 99999971 99999989]
Brun's constant up to 100000000: 1.758815621067975
47.44146728515625 seconds

I did not try higher because I was working in google.colab and I exceeded the allowed limit.

Update

As nicely observed by nocomment below, the process can be made even faster by eliminating the inner loop in the function sieve_of_erastothenes so it becomes:

def sieve_of_eratosthenes(max_num):
    is_prime = np.full(max_num + 1, True, dtype=bool)
    is_prime[0:2] = False
    p = 2
    while (p * p <= max_num):
        if (is_prime[p] == True):
            is_prime[p*p :: p] = False
        p += 1
    return np.nonzero(is_prime)[0]

With this change:

[       2        3        5 ... 99999959 99999971 99999989]
Brun's constant up to 100000000: 1.7588156210679355
12.536579847335815 seconds

for 10^8

Update 2: Execution time estimation

As I mentioned above, I am not working on a machine allowing me to test for 10^20. So, I can give you a method to estimate the time, using a regression method. First start of by computing the execution time for a 6,7,8,9,10. Then fit the curve with a polynomial:

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
import matplotlib.pyplot as plt
import time
import numpy as np

limits = [10**exp for exp in range(2, 10)]  
times = []

for limit in limits:
    start_time = time.time()
    primes = sieve_of_eratosthenes(limit)
    brun_const = brun_constant(primes)
    end_time = time.time()
    execution_time = end_time - start_time
    times.append(execution_time)
    print(f"Limit: {limit}, Time: {execution_time}")

X = np.log(limits).reshape(-1, 1)  
y = np.log(times).reshape(-1, 1)   

poly = PolynomialFeatures(degree=2) 
X_poly = poly.fit_transform(X)

model = LinearRegression()
model.fit(X_poly, y)

X_fit = np.linspace(min(X), max(X), 400).reshape(-1, 1)
X_fit_poly = poly.transform(X_fit)
y_fit = model.predict(X_fit_poly)

X_fit = np.exp(X_fit)
y_fit = np.exp(y_fit)

plt.figure(figsize=(10, 6))
plt.scatter(limits, times, color='blue', label='Actual Times')
plt.plot(X_fit, y_fit, color='red', label='Polynomial Model Fit')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Limits')
plt.ylabel('Execution Time (seconds)')
plt.title('Polynomial Regression Fit')
plt.legend()
plt.show()

log_limit_20 = np.log(np.array([10**20], dtype=np.float64)).reshape(-1, 1)
log_limit_20_poly = poly.transform(log_limit_20)
estimated_log_time_20 = model.predict(log_limit_20_poly)
estimated_time_20 = np.exp(estimated_log_time_20)

print(f"Estimated execution time for limit 10^20: {estimated_time_20[0][0]} seconds")

Which gives

Estimated execution time for limit 10^20: 147470457222.1028 seconds

or 40963960.339500725 hours, or 1706493 days

So, I would suggest using a machine with GPU.