数论之大素数检测与大整数分解(10^18)

大素数检测

大素数检测常用的方法为miller_rabin。网上讲这个方法的文章博客已经很多了，我在这里就只转载一篇别人的文章。

原文出处：http://www.cppblog.com/zoyi-zhang/archive/2008/09/23/62572.aspx

 1 大素数的检验
 2 费马小定理：a^(p-1) mod p = 1（p是素数＆＆a<p＆＆a>0）
 3 
 4  首先我们证明这样一个结论：如果p是一个素数的话，那么对任意一个小于p的正整数a，a, 2a, 3a, …, (p-1)a
 5 除以p的余数正好是一个1到p-1的
 6 排列。例如，5是素数，3, 6, 9, 12除以5的余数分别为3, 1, 4, 2，正好就是1到4这四个数。
 7 反证法，假如结论不成立的话，那么就是说有两个小于p的正整数m和n使得na和ma除以p的余数相同。不妨假设n>m，
 8 则p可以整除a(n-m)。但p是素
 9 数，那么a和n-m中至少有一个含有因子p。这显然是不可能的，因为a和n-m都比p小。
10 用同余式表述，我们证明了：
11 (p-1)! ≡ a * 2a * 3a * … * (p-1)a (mod p)
12     也即：
13 (p-1)! ≡ (p-1)! * a^(p-1) (mod p)
14     两边同时除以(p-1)!，就得到了我们的最终结论：
15 1 ≡ a^(p-1) (mod p)
16 
17 费马小定里的欧拉推广：a^φ(m) ≡ 1 (mod m)（其中φ(m)为与m互质的数的个数）
18 
19 证明与费马小定理类似
20 
21 但是费马小定理的逆命题并不正确，即，当满足a^(p-1) mod p = 1的数p不一定是素数，例如p=341，a=2，
22 而此时p=11＊13
23 
24 后来，人们又发现了561, 645, 1105等数都表明a=2时Fermat小定理的逆命题不成立。虽然这样的数不多，
25 但不能忽视它们的存在。于是，人们把所
26 有能整除2^(n-1)-1的合数n叫做伪素数(pseudoprime)，意思就是告诉人们这个素数是假的。
27 
28 不满足2^(n-1) mod n = 1的n一定不是素数；如果满足的话则多半是素数。这样，一个比试除法效率更高的
29 素性判断方法出现了：制作一张伪素数
30 表，记录某个范围内的所有伪素数，那么所有满足2^(n-1) mod n = 1且不在伪素数表中的n就是素数。之所以
31 这种方法更快，是因为我们可以使用
32 二分法快速计算2^(n-1) mod n 的值，这在计算机的帮助下变得非常容易；在计算机中也可以用二分查找有序
33 数列、Hash表开散列、构建Trie树等
34 方法使得查找伪素数表效率更高。
35 
36 有人自然会关心这样一个问题：伪素数的个数到底有多少？换句话说，如果我只计算2^(n-1) mod n的值，事先
37 不准备伪素数表，那么素性判断出
38 错的概率有多少？研究这个问题是很有价值的，毕竟我们是OIer，不可能背一个长度上千的常量数组带上考场。
39 统计表明，在前10亿个自然数中共
40 有50847534个素数，而满足2^(n-1) mod n = 1的合数n有5597个。这样算下来，算法出错的可能性约为0.00011。
41 这个概率太高了，如果想免去建
42 立伪素数表的工作，我们需要改进素性判断的算法。
43 
44 最简单的想法就是，我们刚才只考虑了a=2的情况。对于式子a^(n-1) mod n，取不同的a可能导致不同的结果。一
45 个合数可能在a=2时通过了测试，
46 但a=3时的计算结果却排除了素数的可能。于是，人们扩展了伪素数的定义，称满足 a^(n-1) mod n = 1的合数n叫
47 做以a为底的伪素数(pseudoprime to base a)
48 。前10亿个自然数中同时以2和3为底的伪素数只有1272个，这个数目不到刚才的1/4。这告诉我们如果同时验证a=2
49 和a=3两种情况，算法出错的概率
50 降到了0.000025。容易想到，选择用来测试的a越多，算法越准确。通常我们的做法是，随机选择
51 
52 若干个小于待测数的正整数作为底数a进行若干次测试，只要有一次没有通过测试就立即把这个数扔回合数的世界。
53 这就是Fermat素性测试。
54 
55 人们自然会想，如果考虑了所有小于n的底数 a，出错的概率是否就可以降到0呢？没想到的是，居然就有这样的合数，
56 它可以通过所有a的测试
57 （这个说法不准确，详见我在地核楼层的回复）。 Carmichael第一个发现这样极端的伪素数，他把它们称作Carmichael数。
58 你一定会以为这样的数一定很大。错。第一个Carmichael 数小得惊人，仅仅是一个三位数，561。前10亿个自然数中
59 Carmichael数也有600个之多。Carmichael数的存在说明，我们还需要继续加强素性判断的算法。
60 
61 Miller和Rabin两个人的工作让Fermat素性测试迈出了革命性的一步，建立了传说中的 Miller-Rabin素性测试算法。新的测
62 试基于下面的定理：如果p是素数，x是小于p的正整数，且x^2 mod p = 1，那么要么x=1，要么x=p-1。这是显然的，因为
63 x^2 mod p = 1相当于p能整除x^2-1，也即p能整除(x+1)(x-1)。由于p是素数，那么只可能是x-1能被p整除(此时x=1)或x+1能
64 被p整除(此时x =p-1)。
65 
66 这就是Miller-Rabin素性测试的方法。不断地提取指数n-1中的因子2，把n-1表示成d*2^r（其中d是一个奇数）。那么
67 我们需要计算的东西就变成了a的d*2^r次方除以n的余数。于是，a^(d * 2^(r-1))要么等于1，要么等于n-1。
68 如果a^(d * 2^(r-1))等于1，定理继续适用于a^(d * 2^(r-2))，这样不断开方开下去，直到对于某个i满足
69 a^(d * 2^i) mod n = n-1或者最后指数中的2用完了得到的a^d mod n=1或n-1。这样，Fermat小定理加强为如下形式：
70 尽可能提取因子 2，把n-1表示成d*2^r，如果n是一个素数，那么或者a^d mod n=1，或者存在某个i使得a^(d*2^i) mod n=n-1 ( 0<=i<r ) 
71 （注意i可以等于0，这就把a^d mod n=n-1的情况统一到后面去了）
72 Miller-Rabin 素性测试同样是不确定算法，我们把可以通过以a为底的Miller-Rabin测试的合数称作以a为底的强伪素
73 数(strong pseudoprime)。
74 第一个以2为底的强伪素数为2047。第一个以2和3为底的强伪素数则大到1 373 653。
75 
76  Miller- Rabin算法的代码也非常简单：计算d和r的值（可以用位运算加速），然后二分计算a^d mod n的值
77 ，最后把它平方r次。程序的代码比想像中的更简单，我写一份放在下边。虽然我已经转C了，但我相信还有很多
78 人看不懂C语言。我再写一次Pascal 吧。函数IsPrime返回对于特定的底数a，n是否是能通过测试。如果函数返回
79 False，那说明n不是素数；如果函数返回True，那么n极有可能是素数。注意这个代码的数据范围限制在longint，
80 你很可能需要把它们改成int64或高精度计算。
81     我们下面来演示一下上面的定理如何应用在Fermat素性测试上。前面说过341可以通过以2为底的Fermat测试，
82 因为 2^340 mod 341=1。如果341真是素数的话，那么2^170 mod 341只可能是1或340；当算得2^170 mod 341确实等于
83 1时，我们可以继续查看2^85除以341的结果。我们发现，2^85 mod 341=32，这一结果摘掉了341头上的素数皇冠，
84 面具后面真实的嘴脸显现了出来
85 
86 
87   对于大数的素性判断，目前Miller-Rabin算法应用最广泛。一般底数仍然是随机选取，但当待测数不太大时，
88 选择测试底数就有一些技巧了。比如，如果被测数小于4 759 123 141，那么只需要测试三个底数2, 7和61就足
89 够了。当然，你测试的越多，正确的范围肯定也越大。如果你每次都用前7个素数(2, 3, 5, 7, 11, 13和17)进行
90 测试，所有不超过341 550 071 728 320的数都是正确的。如果选用2, 3, 7, 61和24251作为底数，那么10^16内
91 唯一的强伪素数为46 856 248 255 981。这样的一些结论使得Miller-Rabin算法在OI中非常实用。通常认为&