用 python 计算 Brun 常数的最佳方法是什么（包括素性测试和长期总和）？

Question

布伦常数： https://en.wikipedia.org/wiki/Brun%27s_theorem http://numbers.computation.free.fr/Constants/Primes/twin.html

如何使用 python 计算高达 10^20 的 Brun 常数 知道素性检查的成本很高，并且将结果求和到 10^20 是一项漫长的任务

这是我的 2 美分尝试：

IsPrime：我知道检查数字是否为素数的最快方法 digital_root ：计算数字的数字根

如果有人知道如何改进才能达到 10^20 的计算量，不客气

import numpy as np
import math

#Brun's constant
#B_2(p) = SUM[P:P+2€P] (1/p + 1/(p+2))
#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...

#Twin prime constant
C_2 = np.float64(0.6601618158468695739278121100145557784326233602847334133194484233354056423)

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
    
    #exclude number not ending with 1 3 7 9
    if (n != 2 and n != 7 and n != 5 and str(n)[len(str(n))-1] not in ("1","3","7","9")): 
        return False

    #check if the number n has a factor contained in the 25 first primes
    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

#approx De Brun's
B_2 = np.float64(0)
B_2Aster = np.float64(0)
lastPrime = 2
lastNonPrime = 1
for p in range(3, 100000000000000000000):
    if IsPrime(p):
        if lastPrime == p-2 and lastNonPrime == p-1: #identify Twin primes
            B_2 = B_2 + (np.float64(np.float64(1)/np.float64(p)) + np.float64(np.float64(1)/np.float64(lastPrime)))
            B_2Aster = B_2 + np.float64(np.float64(np.float64(4)*C_2)/np.float64(np.log(p-2)))
        lastPrime = p
    else:
        lastNonPrime = p
    print(f'p:{p} B_2:{B_2Aster:.52f}',end="\r")

Answer 1

现在，如果您还记得数论课程，有一种叫做“埃拉托斯特尼筛法”的东西，它是一种查找任意给定限制内的所有素数的算法：

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