大整数质合判断

如何快速判断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为类型基数划分整数,有:

整数 性质 整数 性质 整数 性质 整数 性质
24k 2,3的倍数 24k+6 2,3的倍数 24k+12 2,3的倍数 24k+18 2,3的倍数
24k+1 可能为质数 24k+7 可能为质数 24k+13 可能为质数 24k+19 可能为质数
24k+2 2的倍数 24k+8 2的倍数 24k+14 2的倍数 24k+20 2的倍数
24k+3 3的倍数 24k+9 3的倍数 24k+15 3的倍数 24k+21 3的倍数
24k+4 2的倍数 24k+10 2的倍数 24k+16 2的倍数 24k+22 2的倍数
24k+5 可能为质数 24k+11 可能为质数 24k+17 可能为质数 24k+23 可能为质数

可能为质数的有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内的质数吧。

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注

8 + 16 =