请注意,您的主要测试
def prime(x):
if (2**(x-1))%x ==1:
return x
是错误的,它声明了 eg341 = 11*31
和2047 = 23*89
prime。
此外,对于较大x
的 ,产生非常大的中间值,更有效的是
pow(2,x-1,x)
这减少了中间值。
更强的素性检查的适度有效实施:
# The primes below 200 used as bases for the strong Fermat test,
# prime bases have more discriminatory power than composite bases,
# therefore we use prime bases first
bases = [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, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199]
# The strong Fermat test for base b, where n is odd and large,
# n - 1 = m * 2^s with odd m
# the strong test checks whether b**m % n == 1 or
# b**(2**j) % n == n-1 for a 0 <= j < s
# If the test returns False, n is definitely composite, otherwise probably prime
def strongFermat(b,n,m,s):
a = pow(b,m,n)
if a == 1:
return True
n1 = n-1
for i in xrange(s):
if a == n1:
return True
a = (a*a) % n
return False
# Multiple strong Fermat tests, by default use 10 bases
# The probability that a composite passes is less than 0.25**iters
def sppTest(n, iters = 10):
# Assumes n > 1 and with no prime divisors < 200
m = n-1
s = 0
while (m & 1) == 0:
m >>= 1
s += 1
pbases = iters if iters < 47 else 46
for i in xrange(pbases):
if not strongFermat(bases[i],n,m,s):
return False
if pbases < iters:
for i in xrange(iters-pbases):
if not strongFermat(211 + 2*i,n,m,s):
return False
return True
# Trial division to weed out most composites fast
def trialDivisionPrime(n):
if n < 2:
return 0 # Not a prime
if n < 4:
return 2 # definitely prime
if n % 2 == 0 or n % 3 == 0:
return 0 # definitely composite
for d in xrange(5, 200, 6):
if d*d > n:
return 2 # definitely prime
if n % d == 0 or n % (d+2) == 0:
return 0 # composite
return 1 # not yet decided
# The prime test, first do trial division by numbers < 200,
# if that doesn't decide the matter, use some strong Fermat tests
# using 20 tests is the paranoid setting for largish numbers,
# for numbers in 64-bit range, ten or fewer suffice
def prime(n):
td = trialDivisionPrime(n)
return td > 1 or (td == 1 and sppTest(n,20))
# just check a couple of larger numbers
for c in xrange(100000):
if prime(c + 10**25):
print c