在 python 中执行质数对
Performing prime value pairs in python
我目前卡在欧拉计划问题 60 上。问题是这样的：
素数3、7、109、673，很了不起。通过取任意两个质数并以任意顺序连接它们，结果将始终是质数。例如，取 7 和 109，7109 和 1097 都是质数。这四个素数的和，792，代表了一组四个素数的最小和，属性.
找出一组五个素数的最小总和，其中任意两个素数连接产生另一个素数。
在 Python 中，我使用并创建了以下函数。
def isprime(n):
"""Primality test using 6k+-1 optimization."""
if n <= 3:
return n > 1
elif n % 2 == 0 or n % 3 == 0:
return False
i = 5
while i * i <= n:
if n % i == 0 or n % (i + 2) == 0:
return False
i += 6
return True
def prime_pair_sets(n):
heap = [[3,7]]
x = 11
#Phase 1: Create a list of length n.
while getMaxLength(heap)<n:
if isprime(x):
m = [x]
for lst in heap:
tmplst = [x]
for p in lst:
if primePair(x,p):
tmplst.append(p)
if len(tmplst)>len(m):
m=tmplst
heap.append(m)
x+=2
s = sum(maxList(heap))
#Phase 2: Find the lowest sum.
for li in heap:
y=x
while s>(sum(li)+y) and len(li)<n:
b = True
for k in li:
if not primePair(k,y):
b = False
if b == True:
li.append(y)
y+=2
if len(li)>=n:
s = sum(li)
return s
def getMaxLength(h):
m = 0
for s in h:
if len(s) > m:
m = len(s)
return m
def primePair(x, y):
return isprime(int(str(x)+str(y))) and isprime(int(str(y)+str(x)))
def maxList(h):
result = []
for s in h:
if len(s)> len(result):
result = s
return result
我只使用第一阶段执行了 primePairSets 函数。经过大约一个小时的等待，我得到了 74,617 (33647 + 23003 + 16451 + 1493 + 23) 的总和。事实证明，这不是我们要找的总和。我使用这两个阶段再次尝试了这个问题。大约两个小时后，我得到了同样的错误结果。显然，答案小于 74,617。有人可以帮我找到更有效的方法来解决这个问题。请告诉我如何解决。
已找到解决方案
- 素数：13、5197、5701、6733、8389
- 总数：26,033
- 运行 时间减少到约 10 秒，而 OP 解决方案报告的时间为 1-2 小时。
接近
- 基于包含素数列表的扩展路径的回溯算法
- 如果一个素数的成对素数与所有素数都已经在路径中，则可以将其添加到路径中
- 使用埃拉托色尼筛法预先计算我们预期需要的最大素数
- 对于每个素数，预先计算它与哪些其他素数成对素数
代码
import functools
import time
# Decorator for function timing
def timer(func):
"""Print the runtime of the decorated function"""
@functools.wraps(func)
def wrapper_timer(*args, **kwargs):
start_time = time.perf_counter() # 1
value = func(*args, **kwargs)
end_time = time.perf_counter() # 2
run_time = end_time - start_time # 3
print(f"Finished {func.__name__!r} in {run_time:.4f} secs")
return value
return wrapper_timer
#*************************************************************
# Prime number algorithms
#-------------------------------------------------------------
def _try_composite(a, d, n, s):
if pow(a, d, n) == 1:
return False
for i in range(s):
if pow(a, 2**i * d, n) == n-1:
return False
return True # n is definitely composite
def is_prime(n, _precision_for_huge_n=16, _known_primes = set([2, 3, 5, 7])):
if n in _known_primes:
return True
if n > 6 and not (n % 6) in (1, 5):
# Check of form 6*q +/- 1
return False
if any((n % p) == 0 for p in _known_primes) or n in (0, 1):
return False
d, s = n - 1, 0
while not d % 2:
d, s = d >> 1, s + 1
# Returns exact according to http://primes.utm.edu/prove/prove2_3.html
if n < 1373653:
return not any(_try_composite(a, d, n, s) for a in (2, 3))
if n < 25326001:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5))
if n < 118670087467:
if n == 3215031751:
return False
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7))
if n < 2152302898747:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11))
if n < 3474749660383:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13))
if n < 341550071728321:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17))
# otherwise
return not any(_try_composite(a, d, n, s)
for a in _known_primes[:_precision_for_huge_n])
def primes235(limit):
' Prime generator using Sieve of Eratosthenes with factorization wheel of 2, 3, 5 '
# Source: https://rosettacode.org/wiki/Sieve_of_Eratosthenes#Python
yield 2; yield 3; yield 5
if limit < 7: return
modPrms = [7,11,13,17,19,23,29,31]
gaps = [4,2,4,2,4,6,2,6,4,2,4,2,4,6,2,6] # 2 loops for overflow
ndxs = [0,0,0,0,1,1,2,2,2,2,3,3,4,4,4,4,5,5,5,5,5,5,6,6,7,7,7,7,7,7]
lmtbf = (limit + 23) // 30 * 8 - 1 # integral number of wheels rounded up
lmtsqrt = (int(limit ** 0.5) - 7)
lmtsqrt = lmtsqrt // 30 * 8 + ndxs[lmtsqrt % 30] # round down on the wheel
buf = [True] * (lmtbf + 1)
for i in range(lmtsqrt + 1):
if buf[i]:
ci = i & 7; p = 30 * (i >> 3) + modPrms[ci]
s = p * p - 7; p8 = p << 3
for j in range(8):
c = s // 30 * 8 + ndxs[s % 30]
buf[c::p8] = [False] * ((lmtbf - c) // p8 + 1)
s += p * gaps[ci]; ci += 1
for i in range(lmtbf - 6 + (ndxs[(limit - 7) % 30])): # adjust for extras
if buf[i]: yield (30 * (i >> 3) + modPrms[i & 7])
def prime_pair(x, y):
' Checks if two primes are prime pair (i.e. concatenation of two in either order is also a prime)'
return is_prime(int(str(x)+str(y))) and is_prime(int(str(y)+str(x)))
def find_pairs(primes):
' Creates dictionary of what primes can go with others as a pair'
prime_pairs = {}
for i, p in enumerate(primes):
pairs = set()
for j in range(i+1, len(primes)):
if prime_pair(p, primes[j]):
pairs.add(primes[j])
prime_pairs[p] = pairs
return prime_pairs
#*************************************************************
# Main functionality
#-------------------------------------------------------------
@timer
def find_groups(max_prime = 9000, n = 5):
'''
Find group smallest sum of primes that are pairwise prime
max_prime - max prime to consider
n - the size of the group
'''
def fully_connected(p, path):
'''
checks if p is connected to all elements of the path
(i.e. group of primes)
'''
return all(p in prime_pairs.get(path_item, set()) for path_item in path)
def backtracking(prime_pairs, n, path = None):
if path is None:
path = []
if len(path) == n:
yield path[:]
else:
if not path:
for p, v in prime_pairs.items():
if v:
yield from backtracking(prime_pairs, n, path + [p])
else:
p = path[-1]
for t in sorted(prime_pairs[p]):
if t > p and fully_connected(t, path):
yield from backtracking(prime_pairs, n, path + [t])
primes = list(primes235(max_prime)) # Sieve for list of primes up to max_pair
set_primes = set(primes) # set of primes (for easy test if number is prime)
prime_pairs = find_pairs(primes) # Table of primes and set of other primes they can pair with
return next(backtracking(prime_pairs, n), None)
测试运行s
for n, max_prime in [(2, 1000), (3, 1000), (4, 1000), (5, 9000)]:
print(find_groups(max_prime = max_prime, n = n))
运行 次
n max_prime Primes Found Run Time (secs)
2 1000 (3, 7) 0.1280
3 1000 (3, 37, 67) 0.1240
4 1000 (3, 7, 109, 673) 0.1233
5 40,000 (13, 5197, 5701, 6733, 8389) 7.1142
Note: Timings above performed on an older Windows desktop computer,
namely:
~7-year-old HP-Pavilion Desktop i7 CPU 920 @ 2.67 GHz
我根据埃拉托色尼筛法编写了这段代码，它只存储数字 1 mod 6 和 -1 mod 6。速度非常快。
def find_lowest_sum (pmax,num):
n=int(str(pmax)+str(pmax))
sieve5m6 = [True] * (n//6+1)
sieve1m6 = [True] * (n//6+1)
for i in range(1,int((n**0.5+1)/6)+1):
if sieve5m6[i]:
sieve5m6[6*i*i::6*i-1]=[False]*(((n//6+1)-6*i*i-1)//(6*i-1)+1)
sieve1m6[6*i*i-2*i::6*i-1]=[False]*(((n//6+1)-6*i*i+2*i-1)//(6*i-1)+1)
if sieve1m6[i]:
sieve5m6[6*i*i::6*i+1]=[False]*(((n//6+1)-6*i*i-1)//(6*i+1)+1)
sieve1m6[6*i*i+2*i::6*i+1]=[False]*(((n//6+1)-6*i*i-2*i-1)//(6*i+1)+1)
def test_concatenate (p1,p2):
ck=0
p3=int(str(p1)+str(p2))
if sieve1m6[(p3-1)//6] and p3%6==1:
ck=1
elif sieve5m6[(p3+1)//6]and p3%6==5:
ck=1
if ck==1:
p3=int(str(p2)+str(p1))
if sieve1m6[(p3-1)//6] and p3%6==1:
ck=2
elif sieve5m6[(p3+1)//6]and p3%6==5:
ck=2
if ck==2:
return True
else:
return False
kmax=(pmax+1)//6
s1=5*pmax
P1=[]
P2=[]
p=3
kmin=0
while p<s1//5:
if p==3 or(p%6==1 and sieve1m6[(p-1)//6])or(p%6==5 and sieve5m6[(p+1)//6]):
P=[]
if p%6==1:
kmin=p//6
elif p%6==5:
kmin=(p+1)//6
if sieve1m6[kmin]:
if test_concatenate(p,6*kmin+1):
P.append(6*kmin+1)
for k in range(kmin+1,kmax):
if sieve5m6[k]:
if test_concatenate(p,6*k-1):
P.append(6*k-1)
if sieve1m6[k]:
if test_concatenate(p,6*k+1):
P.append(6*k+1)
i1=0
while i1<len(P) and P[i1]<s1:
P1=[p]
s=p
P1.append(P[i1])
s+=P[i1]
for i2 in range(i1+1,len(P)):
for i3 in range(1,len(P1)):
ck=test_concatenate(P[i2],P1[i3])
if ck==False:
break
if len(P1)-1==i3 and ck==True:
P1.append(P[i2])
s+=P[i2]
if len(P1)==num:
if s<s1:
P2=P1
s1=s
break
i1+=1
p+=2
return P2
P=find_lowest_sum(9001,5)
print(P)
s=0
for i in range(0,len(P)):
s+=P[i]
print(s)
我目前卡在欧拉计划问题 60 上。问题是这样的：
素数3、7、109、673，很了不起。通过取任意两个质数并以任意顺序连接它们，结果将始终是质数。例如，取 7 和 109，7109 和 1097 都是质数。这四个素数的和，792，代表了一组四个素数的最小和，属性.
找出一组五个素数的最小总和，其中任意两个素数连接产生另一个素数。
在 Python 中，我使用并创建了以下函数。
def isprime(n):
"""Primality test using 6k+-1 optimization."""
if n <= 3:
return n > 1
elif n % 2 == 0 or n % 3 == 0:
return False
i = 5
while i * i <= n:
if n % i == 0 or n % (i + 2) == 0:
return False
i += 6
return True
def prime_pair_sets(n):
heap = [[3,7]]
x = 11
#Phase 1: Create a list of length n.
while getMaxLength(heap)<n:
if isprime(x):
m = [x]
for lst in heap:
tmplst = [x]
for p in lst:
if primePair(x,p):
tmplst.append(p)
if len(tmplst)>len(m):
m=tmplst
heap.append(m)
x+=2
s = sum(maxList(heap))
#Phase 2: Find the lowest sum.
for li in heap:
y=x
while s>(sum(li)+y) and len(li)<n:
b = True
for k in li:
if not primePair(k,y):
b = False
if b == True:
li.append(y)
y+=2
if len(li)>=n:
s = sum(li)
return s
def getMaxLength(h):
m = 0
for s in h:
if len(s) > m:
m = len(s)
return m
def primePair(x, y):
return isprime(int(str(x)+str(y))) and isprime(int(str(y)+str(x)))
def maxList(h):
result = []
for s in h:
if len(s)> len(result):
result = s
return result
我只使用第一阶段执行了 primePairSets 函数。经过大约一个小时的等待，我得到了 74,617 (33647 + 23003 + 16451 + 1493 + 23) 的总和。事实证明，这不是我们要找的总和。我使用这两个阶段再次尝试了这个问题。大约两个小时后，我得到了同样的错误结果。显然，答案小于 74,617。有人可以帮我找到更有效的方法来解决这个问题。请告诉我如何解决。
已找到解决方案
- 素数：13、5197、5701、6733、8389
- 总数：26,033
- 运行 时间减少到约 10 秒，而 OP 解决方案报告的时间为 1-2 小时。
接近
- 基于包含素数列表的扩展路径的回溯算法
- 如果一个素数的成对素数与所有素数都已经在路径中，则可以将其添加到路径中
- 使用埃拉托色尼筛法预先计算我们预期需要的最大素数
- 对于每个素数，预先计算它与哪些其他素数成对素数
代码
import functools
import time
# Decorator for function timing
def timer(func):
"""Print the runtime of the decorated function"""
@functools.wraps(func)
def wrapper_timer(*args, **kwargs):
start_time = time.perf_counter() # 1
value = func(*args, **kwargs)
end_time = time.perf_counter() # 2
run_time = end_time - start_time # 3
print(f"Finished {func.__name__!r} in {run_time:.4f} secs")
return value
return wrapper_timer
#*************************************************************
# Prime number algorithms
#-------------------------------------------------------------
def _try_composite(a, d, n, s):
if pow(a, d, n) == 1:
return False
for i in range(s):
if pow(a, 2**i * d, n) == n-1:
return False
return True # n is definitely composite
def is_prime(n, _precision_for_huge_n=16, _known_primes = set([2, 3, 5, 7])):
if n in _known_primes:
return True
if n > 6 and not (n % 6) in (1, 5):
# Check of form 6*q +/- 1
return False
if any((n % p) == 0 for p in _known_primes) or n in (0, 1):
return False
d, s = n - 1, 0
while not d % 2:
d, s = d >> 1, s + 1
# Returns exact according to http://primes.utm.edu/prove/prove2_3.html
if n < 1373653:
return not any(_try_composite(a, d, n, s) for a in (2, 3))
if n < 25326001:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5))
if n < 118670087467:
if n == 3215031751:
return False
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7))
if n < 2152302898747:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11))
if n < 3474749660383:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13))
if n < 341550071728321:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17))
# otherwise
return not any(_try_composite(a, d, n, s)
for a in _known_primes[:_precision_for_huge_n])
def primes235(limit):
' Prime generator using Sieve of Eratosthenes with factorization wheel of 2, 3, 5 '
# Source: https://rosettacode.org/wiki/Sieve_of_Eratosthenes#Python
yield 2; yield 3; yield 5
if limit < 7: return
modPrms = [7,11,13,17,19,23,29,31]
gaps = [4,2,4,2,4,6,2,6,4,2,4,2,4,6,2,6] # 2 loops for overflow
ndxs = [0,0,0,0,1,1,2,2,2,2,3,3,4,4,4,4,5,5,5,5,5,5,6,6,7,7,7,7,7,7]
lmtbf = (limit + 23) // 30 * 8 - 1 # integral number of wheels rounded up
lmtsqrt = (int(limit ** 0.5) - 7)
lmtsqrt = lmtsqrt // 30 * 8 + ndxs[lmtsqrt % 30] # round down on the wheel
buf = [True] * (lmtbf + 1)
for i in range(lmtsqrt + 1):
if buf[i]:
ci = i & 7; p = 30 * (i >> 3) + modPrms[ci]
s = p * p - 7; p8 = p << 3
for j in range(8):
c = s // 30 * 8 + ndxs[s % 30]
buf[c::p8] = [False] * ((lmtbf - c) // p8 + 1)
s += p * gaps[ci]; ci += 1
for i in range(lmtbf - 6 + (ndxs[(limit - 7) % 30])): # adjust for extras
if buf[i]: yield (30 * (i >> 3) + modPrms[i & 7])
def prime_pair(x, y):
' Checks if two primes are prime pair (i.e. concatenation of two in either order is also a prime)'
return is_prime(int(str(x)+str(y))) and is_prime(int(str(y)+str(x)))
def find_pairs(primes):
' Creates dictionary of what primes can go with others as a pair'
prime_pairs = {}
for i, p in enumerate(primes):
pairs = set()
for j in range(i+1, len(primes)):
if prime_pair(p, primes[j]):
pairs.add(primes[j])
prime_pairs[p] = pairs
return prime_pairs
#*************************************************************
# Main functionality
#-------------------------------------------------------------
@timer
def find_groups(max_prime = 9000, n = 5):
'''
Find group smallest sum of primes that are pairwise prime
max_prime - max prime to consider
n - the size of the group
'''
def fully_connected(p, path):
'''
checks if p is connected to all elements of the path
(i.e. group of primes)
'''
return all(p in prime_pairs.get(path_item, set()) for path_item in path)
def backtracking(prime_pairs, n, path = None):
if path is None:
path = []
if len(path) == n:
yield path[:]
else:
if not path:
for p, v in prime_pairs.items():
if v:
yield from backtracking(prime_pairs, n, path + [p])
else:
p = path[-1]
for t in sorted(prime_pairs[p]):
if t > p and fully_connected(t, path):
yield from backtracking(prime_pairs, n, path + [t])
primes = list(primes235(max_prime)) # Sieve for list of primes up to max_pair
set_primes = set(primes) # set of primes (for easy test if number is prime)
prime_pairs = find_pairs(primes) # Table of primes and set of other primes they can pair with
return next(backtracking(prime_pairs, n), None)
测试运行s
for n, max_prime in [(2, 1000), (3, 1000), (4, 1000), (5, 9000)]:
print(find_groups(max_prime = max_prime, n = n))
运行 次
n max_prime Primes Found Run Time (secs)
2 1000 (3, 7) 0.1280
3 1000 (3, 37, 67) 0.1240
4 1000 (3, 7, 109, 673) 0.1233
5 40,000 (13, 5197, 5701, 6733, 8389) 7.1142
Note: Timings above performed on an older Windows desktop computer, namely: ~7-year-old HP-Pavilion Desktop i7 CPU 920 @ 2.67 GHz
我根据埃拉托色尼筛法编写了这段代码，它只存储数字 1 mod 6 和 -1 mod 6。速度非常快。
def find_lowest_sum (pmax,num):
n=int(str(pmax)+str(pmax))
sieve5m6 = [True] * (n//6+1)
sieve1m6 = [True] * (n//6+1)
for i in range(1,int((n**0.5+1)/6)+1):
if sieve5m6[i]:
sieve5m6[6*i*i::6*i-1]=[False]*(((n//6+1)-6*i*i-1)//(6*i-1)+1)
sieve1m6[6*i*i-2*i::6*i-1]=[False]*(((n//6+1)-6*i*i+2*i-1)//(6*i-1)+1)
if sieve1m6[i]:
sieve5m6[6*i*i::6*i+1]=[False]*(((n//6+1)-6*i*i-1)//(6*i+1)+1)
sieve1m6[6*i*i+2*i::6*i+1]=[False]*(((n//6+1)-6*i*i-2*i-1)//(6*i+1)+1)
def test_concatenate (p1,p2):
ck=0
p3=int(str(p1)+str(p2))
if sieve1m6[(p3-1)//6] and p3%6==1:
ck=1
elif sieve5m6[(p3+1)//6]and p3%6==5:
ck=1
if ck==1:
p3=int(str(p2)+str(p1))
if sieve1m6[(p3-1)//6] and p3%6==1:
ck=2
elif sieve5m6[(p3+1)//6]and p3%6==5:
ck=2
if ck==2:
return True
else:
return False
kmax=(pmax+1)//6
s1=5*pmax
P1=[]
P2=[]
p=3
kmin=0
while p<s1//5:
if p==3 or(p%6==1 and sieve1m6[(p-1)//6])or(p%6==5 and sieve5m6[(p+1)//6]):
P=[]
if p%6==1:
kmin=p//6
elif p%6==5:
kmin=(p+1)//6
if sieve1m6[kmin]:
if test_concatenate(p,6*kmin+1):
P.append(6*kmin+1)
for k in range(kmin+1,kmax):
if sieve5m6[k]:
if test_concatenate(p,6*k-1):
P.append(6*k-1)
if sieve1m6[k]:
if test_concatenate(p,6*k+1):
P.append(6*k+1)
i1=0
while i1<len(P) and P[i1]<s1:
P1=[p]
s=p
P1.append(P[i1])
s+=P[i1]
for i2 in range(i1+1,len(P)):
for i3 in range(1,len(P1)):
ck=test_concatenate(P[i2],P1[i3])
if ck==False:
break
if len(P1)-1==i3 and ck==True:
P1.append(P[i2])
s+=P[i2]
if len(P1)==num:
if s<s1:
P2=P1
s1=s
break
i1+=1
p+=2
return P2
P=find_lowest_sum(9001,5)
print(P)
s=0
for i in range(0,len(P)):
s+=P[i]
print(s)