大数质因解：浅谈Miller-Rabin和Pollard-Rho算法
2017-07-19 08:54 Amphetamine:能发一下代码吗？
应我那位谜一样好友的邀约，我打算好好看一看Miller-Rabin和Pollard-Rho算法。很奇怪，各种地方有很多代码描述详细过程，但我仍旧很懵。也许是我太弱了，不能从那些“鱼龙混杂”的代码中找出本质上的共性。那么，我们现在来讨论一下吧。
首先，大整数分解现在仍然是个世界级的难题，在“公共密钥”的研究上有着重要的作用。
！！先判断质数！！
试除法：原始的根号算法
额。不想说了。正经一点。
Miller-Rabin：判断一个“p”是否是质数
首先，快速幂是极为必要的。
YYR说，这种东西是要靠RP的。但是，99.999%应该还是足以解决一些问题（或是大多数）了。
a^(p-1)≡1(mod p)，gcd(a,p)=1
费马小定理是极强的质数性质。只要找到一个a，与那个“p”不满足上述条件，就可以跳合数了？答案是肯定的。
但是，如果找不到呢？就是质数？呵呵，很可惜啊。因为存在Carmichael数：无论取多少个a，有一个不满足，算我输。
比如：561 = 11*51就是一个Carmichael数。
为了99.999%，YYR给了一种相对靠谱的方法。
【PART 1：必要的第一次粗筛】
先把2,3,5,7,11,13拿来筛一遍。
这样理论上可以筛掉：1-(1-1/2)*(1-1/3)*(1-1/5)*(1-1/7)*(1-1/11)*(1-1/13)=1-19.18%=80.82%的数。
当然，如果再加上17,19,23,29，理论上可以只会留下15.79%的数，也就是筛去了84.21%的数。
肯定，我们会发现，这个筛率会趋于收敛。
但是，让我们专门卡素数，怎么可能只用这种方法！！！
【PART 2：费马小定理与“二次探测定理”】
在进入这一步之前，我们还不能确定这个“p”的身份。
我们此时会先把p-1搞一搞，不停地除以2。直到发现p-1=2^k*t，而t是奇数。
这样有何用意？且让我们分析。有一个定理：
若a^2≡1(mod p)，则a≡±1(mod p)
YYR说到这里时，我说这很显然，他瞪大了眼睛。因为若p|(a^2-1)，就有p|(a+1)(a-1)。这是经典的欧几里得引理，伟大的欧几里得就是以此证明了唯一分解定理。我们有(a+1)%p==0或(a-1)%p==0。即a≡±1(mod p)。
然后，我们可以去随机地选取一些a。首先用快速幂算出a^t%p，然后把这个数自乘。得到了a^2t%p，然后再自乘。得到了a^4t%p，然后再自乘……我们可以发现，a^t%p自乘k次以后，就可以得到a^(2^k*t)%p，最后就是a^(p-1)%p了。
在这样的一条路上，我们可以充分验证上述“定理”。前面k-1个数中只要出现了1，而它的“父亲”（它的“父亲”自乘得到了它）不是±1，那么就可以跳合数了。最后，第k个数即a^(p-1)%p若不为1，也可以跳合数。
为什么这样效率很高？事实证明，若“p”通过一次测试，不是素数的概率为25%。若“p”通过x次测试，则不是素数的概率为1/(4^x)。当x=5时，素数误判的概率为1/(4^5)=1/1024=0.10%，已经直逼99.90%了。
而假如选取了4个a为2,3,5,7，则在2.5*10^13以内唯一一个失误的数为3215031751。
所以说，选几个a比较好？实践证明，这种东西可以悠着来。
【小结】
1.读入一个"p"。
2.用2,3,5,7,11,13粗筛，发现整除直接跳合数。
3.把p-1搞一搞，不停地除以2。得到奇数t，和一个数k。p-1=2^k*t。
4.随机地找一个数a，算出x=a^t。
5.last=x,x=x*x%p,若x==1且last!=1且last!=p-1，则跳合数。该步骤重复k次。
6.现在x=p-1，若x!=1，则跳合数。
7.按心情回到第4步，99.99%与99.9999%在神犇眼里有着本质区别。
！！然后再乱搞！！
FERMAT：平方差
如果是偶数，就把它釜底抽薪，变成奇数。
对于一个奇数N，若不为质数，则设N=c*d，明显c、d均为奇数。
不妨设c>d，令a=(c+d)/2,b=(c-d)/2，很明显，N=a*a-b*b。这就是了。
我们枚举一个a，有a*a-N=b*b，a>=sqrt(N),若a*a-N是一个完全平方数，就可以递归了。
但是，很明显。这样也是枚举，效率很低，得不偿失。
单纯的POLLARD RHO：生日悖论和判圈
精彩现在继续。有一篇文章很好：点击这里
单纯的Pollard Rho算法非常有趣，容易理解。虽然不是目前最快的算法，但它要比试除法快上多个量级。实现它的思想同样可以用于其他地方。
该算法最早于1975年由John M. Pollard提出，而Richard Brent于1980年提出了改进版本。Wikipedia上说得特别清楚。
上一个FERMAT算法中，我们试图分解N=c*d。那么，现在我们会继续这一思路。假设N只有c和d两个互异质因子，然后在[2,N]中随机出一个数，概率是2/（N-1），很小。
但是，随机选取23个人，他们中存在相同生日的概率大于50%。
这就是生日悖论。
但为什么是对的？如果是对的，然后呢？
【生日悖论的正确性】
我们设现在随机选取n个人(n<=365)，希望能够算出“存在相同生日”的概率。请注意，这个事件是指，我们能够从中找出至少两个生日一样的人。
设该概率为P，很显然它与“不存在相同生日”的概率Q之和为1。而对于Q，我们可以发现，如果人一个一个加入。对于第2个人，要让他与第1个人不同，概率为1-1/365；对于第3个人，要让他与第1个人和第2个人不同，概率为1-2/365；对于第4个人，要让他与第1个人、第2个人和第3个人不同，概率为1-3/365；对于第5个人，要让他与第1个到第4个人不同，概率为1-4/365；……；最后，对于第n个人，要让他与第1个到第n-1个人不同，概率为1-(n-1)/365。最后，即可得Q=(1-1/365)*(1-2/365)*(1-3/365)……(1-(n-1)/365)。
中间有一个不等式：1+x<=ex(x≠0)，我不想证了。
然后，就有Q=e(-1/365)*e(-2/365)*……*e(-(n-1)/365)=e(-n(n-1)/(2*365))
于是我们就可以算出，当n>=23时，他们中存在相同生日的概率大于50%。
其实，如果说一年有N天，那么只要n>=sqrt(N)，他们中存在相同生日的概率就会大于50%。（详见《算导》黑书）
【然后呢？】
现在，我们随机从[2,N]中选出n个数，当好刷出c或d的概率很小，但是如果它们两两作差呢？
使用类似的思路，我们可以猜（当然也可以用完全相同的证法），知道n~sqrt(N)时，概率约为50%。其实仔细想想也很有道理，因为sqrt(N)个数任意组队，有N/2种方案。
虽然很有道理，但是不可能直接两两作差，因为概率还是很低。但是，你可能已经想到了。
如果是gcd(|a-b|,N)呢？于是，事情变得很好看了。
因为N=c*d，c和d都很稀有，但是c和d的倍数便是一波大军。
所以，一个简单的策略如下：
- 在区间[2,N-1]中随即选取n个数，x1,x2, … … , xn
- 判断是否存在gcd(|xi-xj| ,N) >1, 若存在，gcd(|xi-xj| ,N) 是N的一个因子 (c 或 d)
这是很可靠的。如果N有多个因子，也能用类似于FERMAT的方法递归解决。
但是我们算一算。但是很早就出现了一个问题，我们需要选取大约N1/4个数，这个数量太大了，以至于我们并不能将其存放在内存中。
然后，怎么去造这n个数呢？
PollardRho说，很简单。
【流水线】
我们发现，n个数如果先造出来，再去判断实在太傻。于是PollardRho告诉我们，2个数就够了。
我们可以用不断地生成伪随机数。然后像流水线一样判断。
有一个随机函数特别棒，在此举例：f(x)=(x^2+a)%N
这个a已经被你给定了，可以为1当然也能是rand()，只是中途一般不作变化。
代码如下：
1 int find_factorplus(int N) { 2 a = 2; 3 for（ int i= 1; i <= 1000000; i++ ) { 4 b = f(a); 5 p = GCD( abs( b - a ) , N); 6 if( p > 1 ) return p;//Found factor: p 7 a = b; 8 } 9 return 0;//Failed. :-( 10 }
似乎很玄学，但是实际效果确实很棒。但不好的是，伪随机数有着玄学般的循环节。
【floyd来判圈】
对于循环节。你一定会觉得，用个vis数组就好了。似乎很可以，只要你存得下，不揪心就好。
刘汝佳先生说：想象一下，假设有两个小孩子在一个“可以无限向前跑”的跑道上赛跑，同时出发。但其中一个小孩的速度是另一个的两倍。如果跑道是直的，跑得快的小孩永远在前面；但如果跑道有环，则跑得快的小孩将“追上”跑得慢的小孩。
于是就这样了。这叫做“floyd判圈”算法。
1 int find_factorplus(int N) { 2 a = 2; 3 b = a; 4 do { 5 a = f(a);//a runs once 6 b = f(f(b));//b runs twice as fast 7 p = GCD( abs( b - a ) , N); 8 if( p > 1 ) return p;//Found factor: p 9 } while( b != a ); 10 return 0;//Failed. :-( 11 }
这样，我们就可以把退出条件温和化，只要发现有环，那就只有退出了。而不是暴力地把i从1 for 到 1,000,000。
如果算法失败了，我们只需要找到一个新的伪随机函数f(x)或是一个新的a就好了。不过请放心，大多数时候你并不会失败。
最后说一下，代码中a的初值是2，在实际生活中，你并不需要那么讲究，rand()一个也是不错的选择。
“最后”的POLLARD RHO：当与Miller-Rabin发生反应
我们可以发现pollard rho直到现在都还没有与Miller-Rabin有任何联系，但马上就不是了。
对于pollard rho，它可以在Θ(sqrt(p))的时间复杂度内找到N的一个小因子p，这一点我们曾论证过。可见，如果N的因子很多、因子值很小的整数N来说，效率是很优异的。
但是，如果反过来呢？如果说N是大整数，恰好因子很少、因子值很大？
例如，N=2*p，p为质数。你立马发现，N有一个因子2，然后你试图去解决p。然后，这个很优秀的算法成了根号算法。而且直到最后，你都很难判断这个p是否真的不可约。
但是，一旦拥有Miller-Rabin，一切便都已解决。
我们现在可以分析一下复杂度。N的质因子中，超过sqrt(N)的有且仅有一个。这样，即使运气极差，也能有相当的保障。
！！最后总结一下！！
斯堪福说，总结是好习惯。
对于Miller Rabin，我们需要一个快速幂，一个快速乘。先用2,3,5,7,11,13粗筛一遍，再将p的2抽尽，然后随机地选取一些数进行二次探测与费马小定理检验。
对于Pollard Rho，我们需要一个伪随机函数f，一个常数a，一个gcd，一个abs。使用floyd判圈算法。找到一个因子后递归解决，中间判断是否是质数。如果是，做记录。
当我们在做大数质因子分解时，质因子记录完毕后，我们常常会发现这是无序的。这就需要进行一下排序，然后离散化处理出每个质因子出现的次数。这样就解决了。就真的解决了。
致谢
yy儒学长（NOI AG THU）
Amphetamine
IDY002
BIGBALLON
参考
刘汝佳《算法竞赛入门经典》《算法竞赛入门经典训练指南》
林厚从、王宏《信息学奥赛值数学一本通》
丁尧尧idy002《信息学竞赛中的数学知识小结》
BIGBALLON《A Quick Tutorial on Pollard's Rho Algorithm》
