Каков наилучший подход к вычислению константы Бруна с помощью Python (ведьма включает тест на простоту и долгосрочную сумму)?

Теперь, если вы помните свои занятия по теории чисел, есть нечто, называемое решетом Эратосфена, которое представляет собой алгоритм поиска всех простых чисел до любого заданного предела:

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.

Теперь этот алгоритм, переведенный на Python, просто

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)

Итак, ваш полный код будет

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)

Обратите внимание, здесь я выбрал 10^6 в качестве примера:

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

Итак, прошло 0,4 секунды.

За 10^8 я получил

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

Выше я не пробовал, потому что работал в google.colab и превысил разрешенный лимит.

Обновлять

Как хорошо замечено nocomment ниже, процесс можно сделать еще быстрее, исключив внутренний цикл в функции sieve_of_erastothenes, и он будет выглядеть так:

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]

С этим изменением:

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

для 10^8

Обновление 2: оценка времени выполнения

Как я уже упоминал выше, я не работаю над машиной, позволяющей мне тестировать 10^20. Итак, я могу дать вам метод оценки времени, используя метод регрессии. Сначала начнём с вычисления времени выполнения для 6,7,8,9,10. Затем аппроксимируйте кривую полиномом:

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")

Который дает

Estimated execution time for limit 10^20: 147470457222.1028 seconds

или 40963960.339500725 часов, или 1706493 дней

Итак, я бы предложил использовать машину с графическим процессором.