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》