大整数质合判断

如何快速判断1~2^64-1之间某个整数是质数还是合数？在《信息安全数学基础》课上，见到了一些同学写的解决这一问题的代码，也听了老师的点评、讲解，觉得很有收获，记录于此。

我是这样思考这个问题的：首先有这样一个判断一个整数是否是质数的定理：设n是一个正整数，若对所有的质数p<根号n，都有p不整除n，则n一定是质数。
只要遍历所有小于根号n的质数，看其中是否存在一个能整除n，就可以判断n是否为质数。由于遍历质数本身是难以完成的，而所有偶数都一定不是质数，所以只用遍历奇数即可。这样便减少了大概一半的计算量。

看到一个同学对此进行了深化：
所有整数都可以表示成6k+i，i=0,1,2,3,4,5。为描述方便，称6为类型基数。而其中:

• 6k，2,3的倍数
• 6k+1，可能为质数
• 6k+2，2的倍数
• 6k+3，3的倍数
• 6k+4，2的倍数
• 6k+5，可能为质数

故对于比6大的整数，只用判断6k+1和6k+5型的整数即可，这样就减少了三分之二的计算量！

沿着这个思路深继续向前走，考察以24为类型基数划分整数，有：

可能为质数的有8项，8/24=1/3，依旧是排除了三分之二。那么是否最多只能排除排除三分之二呢？显然不是。考察以30为类型基数划分整数，则容易知道，可能是质数的有：

• 30k+1
• 30k+7
• 30k+11
• 30k+13
• 30k+17
• 30k+19
• 30k+23
• 30k+29

共8个。8/30 < 1/3。可以看到+后面的数字分别是是1，再加上30以内的除了30的因子2,3,5外的所有质数。
这一规律对6k+i和24k+i同样成立。以此类推，再结合质数分布会越来越稀，所以取更多的质数相乘，所得的数为类型基数，对整数进行划分，可以排除掉更多的数。
当然这只是纯数学的讨论。

老师说到关于此问题编程的三个层次：

1. 简单地实现算法
2. 多线程，充分利用CPU，内存等硬件资源
3. 多机协作

而我们只停留在第一层。老师又讲了自己的计算方法：

1. 先用普通的算法计算2^16内的质数；(计算全部只需一瞬间)
2. 利用上一步的结果计算2^32内的质数；(计算全部大约10分钟)
3. 利用上一步的结果计算2^64内的质数。

这样经过多次计算就可以较为快速的判断大整数是否是质数了。

下面是我按照这个思路做的过程和结果：

先用普通的算法计算2^16内的质数：

#coding:utf-8
import math
def SimpleIsPrime(n):
    '''不借助外力的、简单的判断n是否是质数，
    是质数返回True，否则返回False'''
    #判断n的类型，如果不是整数，则一定不是素数，返回False
    if not isinstance(n, int):
        return False        
    #若n为-1,0,1则不是整数，返回False
    if n==1 or n==-1 or n==0 :
        return False
    #若n为负数，则将其取反
    if n<0 :
        n = -n
    #若n为2，则是素数，返回True
    if n==2 :
        return True
    #若n为偶数，则返回False
    if n%2==0 :
        return False
    #遍历判断小于根号n的奇数，判断是否整除n
    i = 3
    sqn = math.sqrt(n)
    while i<sqn :
        if n%i==0 :
           return False
        i += 2
    #若无一整除n，则n必为素数，返回True
    return True

f = open("2^16prime.txt", "w")
n = 2
maxn = 2**16
while n<=maxn :
    if SimpleIsPrime(n):
        f.write(str(n)+"\n")
    n += 1
f.close()

为何用Python呢？顺手就用它写了，习惯于用Python快速验证、实现自己的想法。
2^16等于65535是个挺小的数，一瞬间就算完了。输出的2^16prime.txt有30多k大，
用“ wc -l 2^16prime.txt ”命令查看，总共有6595行。

利用上一步的结果计算2^32内的质数：

#coding:utf-8
import math
def IsPrimeByList(n, filename):
    """通过filename给定的质数表判断整数n是否是质数"""
    #判断n的类型，如果不是整数，则一定不是素数，返回False
    if not isinstance(n, int):
        return False        
    #若n为-1,0,1则不是整数，返回False
    if n==1 or n==-1 or n==0 :
        return False
    #若n为负数，则将其取反
    if n<0 :
        n = -n
    #若n为2，则是素数，返回True
    if n==2 :
        return True
    #若n为偶数，则返回False
    if n%2==0 :
        return False
    #遍历判断小于根号n的质数（质数表中给出），判断是否整除n
    sqn = math.sqrt(n)
    flist = open(filename, "r")
    for i in flist:
        i = int(i)
        if n%i==0 :
            flist.close()
            return False
        if i > sqn :
            flist.close()
            return True
    return True

f = open("2^32prime.txt", "w")
n = 2
maxn = 2**32
while n<=maxn :
    if IsPrimeByList(n, "2^32prime.txt"):
        f.write(str(n)+"\n")
    n += 1
f.close()

这次计算花费的时间远远超出我的预计，总用时大约100小时。输出文件2^32prime.txt有2.2G大，共203280868行。把它放在了百度云里：

http://pan.baidu.com/s/1bvIX5O

原本以为是Python效率比较低的原因，但用如下C程序读取并打印1224479671(当时只计算到这里)内的质数，就花了十多分钟。

#include<stdio.h>
int main(){
    char filename[] = "2^32prime.txt";
    FILE *fp;
    char StrLine[22];
    if( (fp = fopen(filename, "r")) == NULL){
        printf("error!\n");
        return -1;
    }
    while(!feof(fp)){
        fgets(StrLine, 22,fp);
        printf("%s\n", StrLine);
    }
    fclose(fp);
    return 0;
}

再联系电脑不到40%的CPU使用率于不到10%的内存使用率，看来我应该学习如何充分利用硬件资源了。
等我学好了再来计算2^64内的质数吧。