《数据结构与算法分析：C语言描述》复习——第十章“算法设计技巧”——质数检验

浏览数：21 / 时间：2015年06月08日

2014.07.07 16:46

简介：

　　对于比较小的正整数n，我们习惯用逐个整除的方法检验n是否为质数。这种算法的复杂度是O(n^0.5)。对于int范围内的整数（最大是2147483647），开方以后不到五万，对于单次计算几乎是一瞬间完成，因此可以接受。但如果n是一个大数，比如10^100，这种算法的耗时就是天文数字了。这时需要开辟新算法来判断质数。

描述：

　　先来说说那么大的质数到底有何用处。当两个一百位的质数p与q相乘，能够过得到一个大概两百位的合数s。如果我不告诉你s有两个质因数p和q，你很可能一辈子也懒得管s到底是质数还是合数。在RSA加密算法中，s=pXq，即使公开了s，一般人也很难得出p和q（只要s够长，就能确保外人暂时无法分解出p和q）。而一旦得出了p和q，就能够解密信息。因此p和q好比两把钥匙，只是它们拼在一起时那条缝好比一条大峡谷（那么大），一个人的眼睛根本无法看清全貌。

　　因此，为了确保选取的两个大数p和q是质数，我们就得设计算法来检测它们。质数的英文是“prime number”，所以指数检测叫“primality testing”。

　　本节的关键知识点“随机算法”，也就是包含有随机数的算法。随机数稍后会在算法中出现。

　　首先介绍费马小定理：如果p是一个质数，并且整数a和p互质，那么a^(p-1)≡1(mod p)。

　　如果把条件限制得更严格，如果p是质数，那么介于0和p之间的整数a肯定是和p互质的。

　　仔细分析这条定理，有如下两点值得注意：

　　　　1. 如果p不是质数，并且a和p互质，那么上述公式也可能满足（此定理是充分而不必要的）

　　　　2. 如果a和p互质，但上述公式不满足，则p一定不是质数（逆否命题与原命题等价）

　　第一条值得深入思考，第二条则是基本的逻辑规律。

　　考虑一个数341，341=11 * 31，所有341是一个合数。2^340≡1(mod 341)，符合上述公式；而3^340≡56(mod 341)，不符合上述公式。

　　这样有可能被误判为质数的合数，被称为Carmichael Number，中文叫伪质数也行。

　　因此可以看出选定不同的a值可能得出不同的结论。因此当你选定的a得到的余数为1时，代表p很可能为质数；若余数不为1，则p肯定不是质数。

　　所以这种检测算法对于不同的(a,p)组合是有一定错误率的，对于某个固定的(a,p)组合，结果则是确定的。如何降低错误率呢？随机选择a值，多测试几次即可。

　　你愿意随机测多次，就能将错误率控制在自己期望的范围内。

　　至此算法的正确性和可行性已经有依据了。那么效率如何能胜任大数计算呢？

　　请观察这个公式：a^(p-1)≡1(mod p)。其中使用的取模运算，幂运算对于大数规模的实现还是比较基础的。而且幂运算通过快速幂的思路可以做到对数级别。

　　例如p是一个100位数，那么log(p)的大致范围只需要O(100)的时间。如此一来，时间开销就从天文数字降到了普通计算机也可接受的级别。

　　已经有不少例子证明了快速幂的强大威力：将连续的乘法（有时不一定是乘法）转化为幂运算，然后用O(log(n))的快速幂算法解决原本需要O(n)时间才能解决的线性问题。

　　这个算法，又为快速幂增加了一个例子。

　　以下的实现没有写大数的各种操作，因为实现篇幅可能会长很多。

实现：

 1 // Primality testing using Fermat‘s Lesser Theorem.
 2 #include <cstdlib>
 3 #include <ctime>
 4 #include <iostream>
 5 using namespace std;
 6 
 7 int primalityTest(int a, int i, int n)
 8 {
 9     // a is randomly chosen between (1, n) to test the primality of n.
10     int x, y;
11     
12     if (i == 0) {
13         return 1;
14     }
15     
16     x = primalityTest(a, i / 2, n);
17     if (x == 0) {
18         // x == 0 means n is composite.
19         // recursively return false
20         return 0;
21     }
22     
23     y = (x * x) % n;
24     if (i % 2 != 0) {
25         y = (y * a) % n;
26     }
27     
28     return y;
29 }
30 
31 bool isPrime(const int n)
32 {
33     if (n < 2) {
34         return false;
35     } else if (n < 4) {
36         return true;
37     }
38     
39     const int NUM_OF_TRIALS = 10;
40     int i = 0;
41     while (i < NUM_OF_TRIALS) {
42         if (primalityTest(2 + rand() % (n - 3), n - 1, n) != 1) {
43             return false;
44         }
45         ++i;
46     }
47     
48     return true;
49 }
50 
51 int main()
52 {
53     srand((unsigned int)time(nullptr));
54     int n;
55     
56     while (cin>> n && n > 0) {
57         cout << (isPrime(n) ? "Yes": "No") << endl;
58     }
59     
60     return 0;
61 }