我正在用Python解决Project Euler问题。我见过一些求解器使用isPrime()函数,这些函数只是测试2 to x ** 0.5中的所有y是否为x % y == 0。这效率不高,我想基于num % 30测试编写一个更好的isPrime()函数。这是我想出来的:
primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
primality = [1, 7, 11, 13, 17, 19, 23, 29]
def isPrime(num):
if not type(num) in (int, long):
raise ValueError
if num in primes:
return True
elif num < 2:
return False
elif num % 30 not in primality:
return False
else:
for prime in primes[3:]:
if num % prime == 0:
return False
seed, sqrt, tryfactor = 1, getIntSquareRoot(num), 1
while tryfactor < sqrt:
for trymod in primality:
tryfactor = seed * 30 + trymod
if num % tryfactor == 0 and not(num == tryfactor):
return False
seed += 1
return True问题7是寻找第1000个素数。因此,我决定让代码将所有这些素数添加到一个列表中,以便后续问题可以参考。我认为给定一个5位数的数字num,num in primes会比重复num % tryfactor的操作快得多。对于其中包含primes[-1] < num < (primes[-1] ** 0.2)的参数,从列表中获取tryfactor值应该比通过tryfactor = seed * 30 + trymod重复生成它们更快。
因此,我想出了以下几点:
def problem7():
count, seed = len(primes), 1
while True:
for modulus in primality:
num = seed * 30 + modulus
if isPrime(num):
count += 1
primes.append(num)
if count > 10000:
return num
seed += 1
def isPrimeB(num):
if not type(num) in (int, long):
raise ValueError
if num in primes:
return True
elif num < 2:
return False
elif num % 30 not in primality:
return False
else:
for prime in primes[3:]:
if num % prime == 0:
return False
seed, sqrt, tryfactor = 1, getIntSquareRoot(num), 1
while tryfactor < sqrt:
for trymod in primality:
tryfactor = seed * 30 + trymod
if num % tryfactor == 0 and not(num == tryfactor):
return False
seed += 1
return True当然,我预计问题7的代码会慢得多,因为生成质数列表需要几秒钟。我还期望针对以后调用isPrime() (例如10、27、35、41和58)的问题的代码会运行得更快。
然而,当问题27、35、41和58的代码变得更慢时,我感到震惊。有人能解释一下为什么在列表中查找值比计算它们慢得多吗?或者是我的代码中有错误?我还能做些什么来提高isPrime()函数的效率呢?
发布于 2014-05-21 20:39:10
它较慢的原因是因为列表查找是O(n)。不使用列表,而是使用集合:
primes = set()
primes.add(num)num in primes校验现在将为O(1)。
同时也忘了这个“优化”:primes[3:]。它实际上减慢了你的代码,因为它重新创建了列表(请注意,如果你切换到集合,它无论如何都不会工作)。
最后,您可以实现the Sieve of Eratosthenes (尽管有更复杂的筛子)来快速生成素数。
发布于 2014-05-21 22:03:07
@Freakish回答了你关于为什么isPrimeB很慢的问题。让我对您所写的内容提出几个替代方案。
下面是我使用2,3,5轮子进行素数检查的版本,它与您的isPrime函数的算法相同,但声明方式有所不同:
def isPrime(n):
d, w, wheel = 2, 0, [1,2,2,4,2,4,2,4,6,2,6]
while d*d <= n:
if n%d == 0: return False
d = d + wheel[w]
w = 3 if w == 10 else w+1
return True有几种方法可以计算第n个素数。一种方法是使用Eratosthenes筛子。数论告诉我们,第n个素数总是小于n(log + log ),以e为底的对数,所以你可以筛选到极限并丢弃多余的素数。下面是筛子:
def primes(n):
m = (n-1)//2; b = [True] * m
i, p, ps = 0, 3, [2]
while i < m:
if b[i]:
ps.append(p)
for j in range(2*i*i+6*i+3, m, p):
b[j] = False
i, p = i+1, p+2
return ps所以,为了得到第1000个素数:
>>> 10000 * ( log(10000) + log(log(10000)) )
114306.67178344031
>>> (primes(114306))[10000]
104743另一种选择是使用2,3,5,7轮生成候选素数,并通过对3个基数2,7,61进行Miller-Rabin检验来确认它们的素性,这对于小于2^32 (实际上,稍大一点)的素数是足够的:
def genPrimes(): # valid to 2^32
def isPrime(n):
def isSpsp(n, a):
d, s = n-1, 0
while d%2 == 0:
d /= 2; s += 1
t = pow(a,d,n)
if t == 1: return True
while s > 0:
if t == n-1: return True
t = (t*t) % n; s -= 1
return False
for p in [2, 7, 61]:
if n % p == 0: return n == p
if not isSpsp(n, p): return False
return True
w, wheel = 0, [1,2,2,4,2,4,2,4,6,2,6,4,2,4,\
6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,\
2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10]
p = 2; yield p
while True:
p = p + wheel[w]
w = 4 if w == 51 else w + 1
if isPrime(p): yield p然后,可以通过表达式next(itertools.islice(genPrimes(), n, n+1))来计算第n个素数
>>> next(itertools.islice(genPrimes(), 10000, 10001))
104743这两种方法都会在您按下enter键后立即返回第1000个素数。
如果你对质数编程感兴趣(或者你只想解决Project Euler中的质数问题),你可能会对我博客上的this essay或these entries感兴趣。
https://stackoverflow.com/questions/23783236
复制相似问题