xiaoguozi's Blog
Pay it forword - 我并不觉的自豪,我所尝试的事情都失败了······习惯原本生活的人不容易改变,就算现状很糟,他们也很难改变,在过程中,他们还是放弃了······他们一放弃,大家就都是输家······让爱传出去,很困难,也无法预料,人们需要更细心的观察别人,要随时注意才能保护别人,因为他们未必知道自己要什么·····

第一章 素数判定法现状

现在,确定性素数判定法已经有很多种,常用的有试除法、威廉斯方法、艾德利曼和鲁梅利法。它们的适用范围各不相同,威廉斯方法比较适合10^2010^50之间的数,艾德利曼和鲁梅利法适合大于10^50的数,对于32位机器数,由于都小于10^10,所以一般都用试除法来判定。

也许有人会问:“你为什么没有提马宁德拉.阿格拉瓦法呢?不是有人说它是目前最快的素数判定法吗?” 其实这是一个很大的误解,阿格拉瓦法虽然是log(n)的多项式级算法,但目前只有理论上的意义,根本无法实用,因为它的时间复杂度是O(log(n)^12),这个多项式的次数太高了。就拿最慢的试除法跟它来比吧,试除法的时间复杂度为O(n^(1/2)*log(n)^2),当n = 16时,log(n)^12 = 16777216,而n^(1/2)*log(n)^2 = 64,你看相差有多么大!如果要让两者速度相当,即log(n)^12 = n^(1/2)*log(n)^2,得出n = 10^43.1214,此时需要进行的运算次数为log(n)^12 = 10^25.873(注意:本文中log()函数缺省以2为底),这样的运算次数在一台主频3GHz的计算机上运行也要10^8.89707年才能运行完,看来我们这辈子是别指望看到阿格拉瓦法比试除法快的这一天啦!

除了这些确定性素数判定法外,还有基于概率的非确定性素数判定法,最常用的就是米勒-拉宾法。

对于32位机器数(四则运算均为常数时间完成),试除法的时间复杂度是O(n^(1/2)),而米勒-拉宾法的时间复杂度只有O(log(n))。所以后者要比前者快得多,但是由于米勒-拉宾法的非确定性,往往我们在需要确定解时仍然要依靠速度较慢的试除法。那是否可以通过扩展米勒-拉宾法,来找到一种更快的确定性素数判定法呢?结论是肯定的,本文就带你一起寻找这样一种方法。

第二章 2-伪素数表的生成

既然要扩展米勒-拉宾法,那首先我们应该知道为什么米勒-拉宾法是个非确定性素数判定法?答案很简单,由于伪素数的存在。由于米勒-拉宾法使用费尔马小定理的逆命题进行判断,而该逆命题对极少数合数并不成立,从而产生误判,这些使费尔马小定理的逆命题不成立的合数就是伪素数。为了研究伪素数,我们首先需要生成伪素数表,原理很简单,就是先用筛法得出一定范围内的所有素数,然后逐一判定该范围内所有合数是否使以2为底数的费尔马小定理的逆命题不成立,从而得出该范围内的2-伪素数表。我的程序运行了100分钟,得出了32位机器数范围内的2-伪素数表,如下:

341

561

645

1105

1387

1729

1905

2047

2465

2701

...

...

...

4286813749

4288664869

4289470021

4289641621

4289884201

4289906089

4293088801

4293329041

4294868509

4294901761

(共10403个,由于篇幅所限,中间部分省略。)

第三章 寻找2-伪素数的最小可判定底数

对于2-伪素数表的每一个伪素数,寻找最小的可以判定它们是合数的底数,我把这个底数称之为最小可判定底数。特别地,对于绝对伪素数,它的最小质因子即是它的最小可判定底数。由于已经证明了绝对伪素数至少有三个质因子,所以这个最小质因子一定不大于n^(1/3)。下面就是我找到的最小可判定底数列表:

341     3

561     3

645     3

1105    5

1387    3

1729    7

1905    3

2047    3

2465    5

2701    5

...

...

...

4286813749      3

4288664869      3

4289470021      5

4289641621      3

4289884201      3

4289906089      3

4293088801      3

4293329041      3

4294868509      7

4294901761      3

通过统计这个列表,我发现了一个规律,那就是所有的最小可判定底数都不大于n^(1/3),由前述可知,对于绝对伪素数,这个结论显然成立。而对于非绝对伪素数,虽然直观上觉得它应该比绝对伪素数好判定出来,但是我无法证明出它的最小可判定底数都不大于n^(1/3)。不过没关系,这个问题就作为一个猜想留给数学家来解决吧,更重要的是我已经通过实验证明了在32位机器数范围内这个结论成立。

我们还有没有更好的方法来进一步减小最小可判定底数的范围呢?有的!我们可以在计算平方数时进行二次检测,下面是进行了二次检测后重新计算的最小可判定底数列表:

341     2

561     2

645     2

1105    2

1387    2

1729    2

1905    2

2047    3

2465    2

2701    2

...

...

...

4286813749      2

4288664869      2

4289470021      2

4289641621      2

4289884201      2

4289906089      2

4293088801      2

4293329041      2

4294868509      2

4294901761      3

很显然,二次检测是有效果的,经过统计,我发现了新的规律,那就是经过二次检测后所有的最小可判定底数都不大于n^(1/6),真的是开了一个平方呀,哈哈!这个结论的数学证明仍然作为一个猜想留给数学家们吧。我把这两个猜想叫做费尔马小定理可判定上界猜想。而我已经完成了对32位机器数范围内的证明。

通过上面总结的规律,我们已经可以设计出一个对32位机器数进行素数判定的 O(n^(1/6)*log(n)) 的确定性方法。但是这还不够,我们还可以优化,因为此时的最小可判定底数列表去重后只剩下了5个数(都是素数):{2,3,5,7,11}。天哪,就是前5个素数,这也太容易记忆了吧。

不过在实现算法时,需要注意这些结论都是在2-伪素数表基础上得来的,也就是说不管如何对2的判定步骤必不可少,即使当2>n^(1/6)时。

还有一些优化可以使用,经过实验,当n>=7^6时,可以不进行n^(1/6)上界限制,而固定地用{2,5,7,11}去判定,也是100%正确的。这样就可以把判定次数降为4次以下,而每次判定只需要进行4log(n)次乘除法(把取余运算也看作除法),所以总的计算次数不会超过16log(n)。经过实验,最大的计算次数在n=4294967291时出现,为496次。


自己写了个随机化素数判定算法:
 1 #include <iostream>
 2 #include <cmath>
 3 #include <algorithm>
 4 
 5 using namespace std;
 6 typedef unsigned __int64 longlong_t;
 7 
 8 bool _IsLikePrime(longlong_t n, longlong_t base)
 9 {
10     longlong_t power = n-1;
11     longlong_t result = 1;
12     longlong_t x = result;
13     longlong_t bits = 0;
14     longlong_t power1 = power;
15     //统计二进制位数
16     while (power1 > 0)
17     {
18         power1 >>= 1;
19         bits++;
20     }
21     //从高位到低位依次处理power的二进制位
22     while(bits > 0)
23     {
24         bits--;
25         result = (x*x)%n;
26         //二次检测
27         if(result==1&&x!=1&&x!=n-1)
28             return false;
29 
30         if ((power&((longlong_t)1<<bits)) != 0)
31             result = (result*base)%n;
32 
33         x = result;
34     }
35     //费尔马小定理逆命题判定
36     return result == 1;
37 }
38 inline bool JudgePrime(int n,int s)
39 {
40     srand(20);
41     for(int i=0;i<s;i++){
42         int a=rand()%(n-1)+1;
43         if(!_IsLikePrime(n,a))
44             return false;
45     }
46     return true;
47 }
48 int main()
49 {
50     int n;
51     while(scanf("%d",&n)!=EOF){
52         int num=0;
53         int cnt=0;
54         for(int i=0;i<n;i++){
55             scanf("%d",&num);
56             if(JudgePrime(num,20))
57                 cnt++;
58         }
59         printf("%d\n",cnt);
60     }
61     return 0;
62 }
s越大,稳定性越好,但效率会低点。原作者有更好的判定算法,对作者表示感谢。。
posted on 2008-01-26 09:37 小果子 阅读(1195) 评论(0)  编辑 收藏 引用 所属分类: 学习笔记

只有注册用户登录后才能发表评论。
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理