2012年5月

在了解验证与获得素数前，先来了解一下什么是素数：

质数，又称素数，指在一个大于1的自然数中，除了1和此整数自身外，无法被其他自然数整除的数（也可定义为只有1和本身两个因数的数）。比1大但不是素数的数称为合数。1和0既非素数也非合数。素数在数论中有着很重要的地位。最小的素数是2，也是素数中唯一的偶数（双数）；其他素数都是奇数（单数）。质数有无限多个，所以不存在最大的质数。围绕著素数存在很多问题、猜想和定理。著名的有孪生素数猜想和哥德巴赫猜想。素数序列的开头是这样的：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

根据素数的定义想到的一些解决办法

从上面我们可以知道，任何一个大于 1 且除1和自身之外不能被其它整数整除的整数为素数，那么我们首先可能就会想到一个办法来验证一个整数是不是素数：给定一个整数 n ，让 n 依次除以 1 -> n 之间的所有整数，如果都没有能整除的，那么这个数肯定是素数：

def factor(n):
'''Get number's factors'''
end = n+1
for i in xrange(1,end):
if n % i == 0 : yield i
def force_checker(n):
'''Checker weither a number is a prime number'''
if len([i for i in factor(n)]) == 2: return True
return False

上面的代码我们作了什么？其实很简单，根据前面我们的分析，对于任何一个需要验证的数字，我们首先需要得到它的因子，这个工作是由 factor() 函数来完成的，遍历 1 -> n 之间的所有整数，一个一个测试，如果是这个数的因子，则将其返回（ yield i ）这一句，所以 factor() 是一个“生成器”:/code/python-functions-decorators ，有了这个之后，我再写了一个 force_checker() 函数，来暴力验证，办法还是最原始的，先找到整数 n 的所有因数子，然后检查其个数，如果为两个（也就是 1 与其自身），则这个数是素数。

上面的这两个函数确实是已经可以实现素数的验证了，但是对于一整数，还没有任何关系，但是如果整数过大，那就麻烦了，你的电脑将有很长的时间，CPU 处于 100%的状态，我们可以来作一个测试，来查看一下验证需要多少时间：

>>> from time import time
>>> def test(n):
... print(time())
... print force_checker(n)
... print(time())
...
>>> test(104729)
1317222986.09
True
1317222986.1
>>> test(1047290)
1317222997.53
False
1317222997.66
>>> test(10472900)
1317223000.42
False
1317223001.55

从上面我们可以看到，对于6位的这个数字，整个验证花了不到0.01 秒，我在 104729 后面追加了一个 0 之后，花了 0.13 秒，但是当我再追加一个 0 之后，就花了 1.13 秒了，经过我自己电脑上面的验证，对于大于 10 位的数字，这种办法根本就行不通，所以，还得想其它办法。

我再想到是，先对数字进行过滤，我想到的是下面这些数字其实根本就不需要去测的：

  1. 未位为（0，2，4，5，6，8）的数字都能被 2 整除，所以肯定不是素数
  2. 每一位相加得到的数字如果能被 3 整除，那么这个数也能被 3 整除
  3. 任何一个大于一个整数 ½ 的整数，绝对不是该整数的因子

跟着上面的思路，我想到了先得修改一下下 force_checker() 函数：

def force_checker(n):
'''Checker weither a number is a prime number'''
if n in (2,3,5): return True
if str(n)[-1:] in ('0','2','4','5', '6','8') or eval('+'.join(str(n))) % 3 == 0: return False
if len([i for i in factor(n)]) == 2: return True
return False

来看看现在的效果：

>>> test(10472900)1317223764.19
False
1317223764.19
>>> test(104729001)
1317223822.58
False
1317223822.58

几乎是不花时间就能检测出我所提供的这两个数字是否为素数，但是你可能也看到了，我其实是取了一个特别的数字，一个是 0 结尾的，一个是各位数字之和为 3 的整数倍的，但是我们还是可以通过另一个办法来看看加了上面的一行判断之后的实际效果，我们来生成一个素数列，看看需要花多少时间：

我们先来写一个用来生成素数数列的函数： prime_generator() ，它是一个生成器，同时，为了以后的测试方法，我还写了一个检测运行时间的函数 test() ，现在我的所有函数是下面这样的：

import sys
def factor(n):
'''Get number's factors'''
end = n+1
for i in xrange(1,end):
if n % i == 0 : yield i
def force_checker_old(n):
if len([i for i in factor(n)]) == 2: return True
return False
def force_checker(n):
if n in (2,3,5): return True
if str(n)[-1:] in ('0','2','4','5', '6','8') or eval('+'.join(str(n))) % 3 == 0: return False
if len([i for i in factor(n)]) == 2: return True
return False
def prime_generator(limit = None, start = 1, end = None, checker = force_checker):
seek, counter = start, 0
while True:
if limit is not None and counter >= limit: break
if end is not None and seek >= end: break
if checker(seek):
yield seek
seek += 1
counter += 1
else:
seek += 1
def test(limit = None, start = 1, end = None, checker = force_checker):
import sys
from time import time
print('Generate prime numbers with {} checker.'.format(checker.__name__))
time_start = time()
column, counter = 1, 1
for prime in prime_generator(limit,start,end,checker):
sys.stdout.write('{:>5} : {:<10}'.format(counter, prime))
counter += 1
if column % 5 == 0:
sys.stdout.write('n')
column += 1
time_end = time()
print('nnI got {} results in {} seconds(form {} to {}).'.format(
counter-1,time_end - time_start, time_start, time_end))
def main():
test(1000)
test(1000,checker=force_checker_old)
if __name__ == '__main__':
main()

上面代码中 prime_generator() 就是一个素数生成器，我们可以在 for 循环中由 limit 限制的数量个素数，如果 limit 设置为 None，则将不限制素数个数（在 for 循环中，不建议这么干）， start 表示我们从哪一个数开始生成，默认为 1 ，表示需要生成 大于 1 的素数， end 默认为 None ，表示不限制最大素数（不是最大素数个数），由于我们需要测试不通的素数测试方法的性能，所以我这里使用了一个 checker 参数，默认使用 force_checker() ，但是我们可以指定不同的验证函数，只要该函数是同样的使用方法（传入一个数，返回 True 或者 False）。

test() 函数是我写的一个测试函数，我们将由它来调用 prime_generator() 函数，并且记录整个程序的运行时间与素数个数，我还将上面的所有代码都保存在一个名为 prime.py 的文件中，这样我可以直接使用我的 TextMate 编辑器运行它，或者你也可以使用 ：

$python /path/to/prime.py

命令来运行，我们来测试一下吧，在我的机器上面，同时用两个验证器来验证素数，生成 1000 条记录的时间如下：

Generate prime numbers with force_checker checker.
1 : 2 2 : 3 3 : 5 4 : 7 5 : 11
...... 中间部分省略 ......
996 : 7879 997 : 7883 998 : 7901 999 : 7907 1000 : 7919
I got 1000 results in 1.05513501167 seconds(form 1317225770.58 to 1317225771.64).

Generate prime numbers with force_checker_old checker.
1 : 2 2 : 3 3 : 5 4 : 7 5 : 11
...... 中间部分省略 ......
996 : 7879 997 : 7883 998 : 7901 999 : 7907 1000 : 7919
I got 1000 results in 3.41558718681 seconds(form 1317225771.64 to 1317225775.05).

上面的结果一目了然了吧，整整快了2.4秒（机器不同，运行时间与环境不同，可能时间有些许不一样，你的结果也可能和我的不一样，但是绝对是快，而且快很多），在我这里测试获取 2000 个素数的结果是：优化之后的函数使用 4.46425414085 秒，未优化的使用了 15.689232111 秒，获取个数不断增加，它们之间的时间差将成倍的增加。当然，上面优化过后的函数也不会是最好的，为什么我给这个函数取名叫作 force_checker() ？那是因为它是最笨的办法，暴力验证，就是不管你是什么数，我一个一个试，只到试完为止，这是一种没有任何技巧可言的方法。

Miller Rabin 验证法

上面的方法对于小整数而言，速度已经很快了，但是在我们现在的现实应用中，基本上不可能使用上在的这个办法，比如目前素数应用最广的应用领域公共密钥体系中，一般选择的素数都相当的大，至少都在100位以上，而上面的这个办法，对于10位的数字就已经力不从心了，而对于100位的数字，可能需要耗尽你一生的时间才能验证成功，所以，还得选择其它的方法，比如最有名的一个验证素数的办法：Rabin Miller，按我自己的想法，这个属性概率验证法吧，如果证明这个方法的有效性不是我所研究的对象，我也没这个能力，但是至少它是被所有人都认同的快速的验证法。

Miller Rabin 验证法的细节是这样的：

  1. 选择一个代测的随机数 p ，计算 b ， b 是 2 整除 p - 1 的次数，然后计算 m ，使用 n = 1 + (2**m) 。
  2. 选择一个小于 p 的随机数 a 。
  3. 设 j = 0 且 z = (a**m) mod p 。 # 如果 z = 1@ 或者 z = p - 1 ，那么 p 通过测试，可能是素数
  4. 如果 j > 0 且 z = 1 ，那么 p 不是素数
  5. 设 j = j + 1 如果 j < b 且 z != p - 1 ， 设 z = (z**m) mod p ，然后咽到上面步骤，如果 z = p - 1 ，那么 p 通过测试，可能为素数
  6. 如果 j = b 且 z != p - 1 ，不是素数

数 a 被当成证据的概率为 75% ，这意味着我们在上面的测试过程中，迭代次数为 t 时，产生一个假素数所花费的时间不会超过 1/4**t ，当我们迭代 5 次时，产生素数的概率已经高达 99.999% 了， 这几乎可以说是一个素数了，而国际上通用的作用是，必须迭代50次，下面我们来在我们的prime.py 中实现它：

def miller_rabin_checker(n):
'''Checker wiether a number is a prime number via rabin miller method'''
if n <= 1: return False
if n in (2,3,5): return True
if str(n)[-1:] in ('0','2','4','5', '6','8') or eval('+'.join(str(n))) % 3 == 0: return False
import sys, random
def to_binary(n):
r = []
while n > 0:
r.append(n % 2)
n = n / 2
return r
def test(a, n):
'''test(a, n) -> bool Test whether n is complex
Returns:
- True : if n is complex
- False: if n is probably prime
'''
b = to_binary(n-1)
d = 1
for i in xrange(len(b) - 1, -1, -1):
x = d
d = (d * d) % n
if d == 1 and x != 1 and x != n - 1:
return True
if b[i] == 1:
d = (d * a) % n
if d != 1:
return True
return False
def checker(n, s = 50):
'''checker(n, s = 50) -> bool checks whether n is prime or not
Returns:
- True : if n is probably prime.
- False: if n is complex.
'''
for j in xrange(1, s + 1):
a = random.randint(1, n - 1)
if test(a,n):
return False
return True
return checker(n)

我们同样可以使用 test() 函数来验证上面的这个函数的效率：

Generate prime numbers with miller_rabin_checker checker.
1 : 2 2 : 3 3 : 5 4 : 7 5 : 11
...... 中间部分省略 ......
96 : 503 97 : 509 98 : 521 99 : 523 100 : 541
I got 100 results in 0.0682229995728 seconds(form 1317228097.28 to 1317228097.35).

Generate prime numbers with force_checker checker.
1 : 2 2 : 3 3 : 5 4 : 7 5 : 11
...... 中间部分省略 ......
96 : 503 97 : 509 98 : 521 99 : 523 100 : 541
I got 100 results in 0.0126042366028 seconds(form 1317228097.35 to 1317228097.38).

Generate prime numbers with force_checker_old checker.
1 : 2 2 : 3 3 : 5 4 : 7 5 : 11
...... 中间部分省略 ......
96 : 503 97 : 509 98 : 521 99 : 523 100 : 541
I got 100 results in 0.0203368663788 seconds(form 1317228097.38 to 1317228097.39).

上面的数据显示运行时间由短到长的验证方法分别为：

force_checker < force_checker_old < miller_rabin_checker

这是怎么回事儿？其实很容易解释，对于小整数，force_checker_old 所作的计算量更少而已，但是如果对于大整数，由于 force_checker_old 或者 force_checker 所用的方法都是试除法，就是所有小于 n 的数字都拿去除一次，所以，待测整数越大，其时间就越多，而且时间是成倍的增加，来看下面的验证（下面的测试我将不再验证 force_checker_old() 函数）：

生成1000个素数

I got 1000 results in 0.895807027817 seconds(form 1317228892.57 to 1317228893.47).

I got 1000 results in 3.35119318962 seconds(form 1317228898.47 to 1317228901.82).

生成10000个素数

Generate prime numbers with miller_rabin_checker checker.
1 : 2 2 : 3 3 : 5 4 : 7 5 : 11
...... 中间部分省略 ......
9996 : 104707 9997 : 104711 9998 : 104717 9999 : 104723 10000 : 104729
I got 10000 results in 11.3178238869 seconds(form 1317228989.72 to 1317229001.03).

Generate prime numbers with force_checker_old checker.
1 : 2 2 : 3 3 : 5 4 : 7 5 : 11
...... 中间部分省略 ......
9996 : 104707 9997 : 104711 9998 : 104717 9999 : 104723 10000 : 104729
I got 10000 results in 571.780455112 seconds(form 1317229006.04 to 1317229577.82).

这个差距就有点儿太大了哈，不过，这就是一个事实。

完整的程序文件： prime.py

#!/usr/bin/env python
# encoding: utf-8
"""
File: prime.py
Created by Pan Tao on 2011-09-28.
Copyright (c) 2011 CosTony.Com. All rights reserved.
"""
import sys
def factor(n):
'''Get number's factors'''
end = n+1
for i in xrange(1,end):
if n % i == 0 : yield i
def force_checker_old(n):
'''Checker weither a number is a prime number'''
if len([i for i in factor(n)]) == 2: return True
return False
def force_checker(n):
'''Checker weither a number is a prime number'''
if n <= 1: return False
if n in (2,3,5): return True
if str(n)[-1:] in ('0','2','4','5', '6','8') or eval('+'.join(str(n))) % 3 == 0: return False
if len([i for i in factor(n)]) == 2: return True
return False
def miller_rabin_checker(n):
'''Checker wiether a number is a prime number via rabin miller method'''
if n <= 1: return False
if n in (2,3,5): return True
if str(n)[-1:] in ('0','2','4','5', '6','8') or eval('+'.join(str(n))) % 3 == 0: return False
import sys, random
def to_binary(n):
r = []
while n > 0:
r.append(n % 2)
n = n / 2
return r
def test(a, n):
'''test(a, n) -> bool Test whether n is complex
Returns:
- True : if n is complex
- False: if n is probably prime
'''
b = to_binary(n-1)
d = 1
for i in xrange(len(b) - 1, -1, -1):
x = d
d = (d * d) % n
if d == 1 and x != 1 and x != n - 1:
return True
if b[i] == 1:
d = (d * a) % n
if d != 1:
return True
return False
def checker(n, s = 50):
'''checker(n, s = 50) -> bool checks whether n is prime or not
Returns:
- True : if n is probably prime.
- False: if n is complex.
'''
for j in xrange(1, s + 1):
a = random.randint(1, n - 1)
if test(a,n):
return False
return True
return checker(n)
def prime_generator(limit = None, start = 2, end = None, checker = force_checker):
'''Generator of prime numbers
varibles:
- limit : how many prime numbers
- start : where to start
- end : when to stop
- checker : the prime number checker, default to force_checker
'''
seek, counter = start, 0
while True:
if limit is not None and counter >= limit: break
if end is not None and seek >= end: break
if checker(seek):
yield seek
seek += 1
counter += 1
else:
seek += 1
def test(limit = None, start = 2, end = None, checker = force_checker):
import sys
from time import time
print('Generate prime numbers with {} checker.n'.format(checker.__name__))
time_start = time()
column, counter = 1, 1
for prime in prime_generator(limit,start,end,checker):
sys.stdout.write('{:>5} : {:<10}'.format(counter, prime))
counter += 1
if column % 5 == 0:
sys.stdout.write('n')
column += 1
time_end = time()
print('nnI got {} results in {} seconds(form {} to {}).'.format(
counter-1,time_end - time_start, time_start, time_end))
def main():
import sys
checkers = {'fo': force_checker_old, 'f': force_checker, 'mr': miller_rabin_checker, 'e': ''}
while True:
checker_help = '''
Which checker you want to use?
type 'fo' to choose force_checker_old checker
type 'f' to choose force_checker checker
type 'mr' to chose miller_rabin_checker
type 'e' to exit
'''
checker = raw_input(checker_help)
if checker not in checkers:
print('You can only type "fo", "f" or "rm" to choose checker, type again :n')
continue
if checker == 'e':
exit()
action_help = '''
What are you want to do?
type 'g' to generate prime numbers
type 'c' to check whether a number is a primen
'''
action = raw_input(action_help)
if action == 'g':
limit = raw_input('How many numbers do you want to generate?n')
start = raw_input('The min-number of the prime list is:n')
end = raw_input('The max number of the prime list is :n')
if int(end) <= int(start):
print('The end must bigger than start')
continue
if int(limit) < 1:
print('limit must bigger or equal to 1')
continue
test(int(limit), int(start), int(end), checkers[checker])
elif action == 'c':
if checkers[checker](int(raw_input('Enter the number you want to check:n'))):
print('This number is prime number')
else:
print('This number is not a prime number')
if __name__ == '__main__':
main()