64位以内Rabin-Miller 强伪素数测试和Pollard rho 因数分解算法的实现

64位以内Rabin-Miller强伪素数测试和Pollard因数分解算法的实现在求解POJ1811题PrimeTest中应用到的两个重要算法是Rabin-Miller强伪素数测试和Pollard因数分解算法。前者可以在nOlog的时间内以很高的成功概率判断一个整数是否是素数。后者可以在最优4nO的时间内完成合数的因数分解。这两种算法相对于试除法都显得比较复杂。本文试图对这两者进行简单的阐述，说明它们在32位计算机上限制在64位以内的条件下的实现中的细节。下文提到的所有字母均表示整数。一、Rabin-Miller强伪素数测试Rabin-Miller强伪素数测试的基本思想来源于如下的Fermat小定理：如果p是一个素数，则对任意a有paapmod。特别的，如果p不能整除a，则还有papmod11。利用Fermat小定理可以得到一个测试合数的有力算法：对1n，选择1a，计算nanmod1，若结果不等于1则n是合数。若结果等于1则n可能是素数，并被称为一个以a为基的弱可能素数（weakprobableprimebasea，a-PRP）;若n是合数，则又被称为一个以a为基的伪素数（pseudoprime）。这个算法的成功率是相当高的。在小于25,000,000,000的1,091,987,405个素数中，一共只用21,853个以2为基的伪素数。但不幸的是，Alford、Granville和Pomerance在1994年证明了存在无穷多个被称为Carmichael数的整数对于任意与其互素的整数a算法的计算结果都是1。最小的五个Carmichael数是561、1,105、1,729、2,465和2,801。考虑素数的这样一个性质：若n是素数，则1对模n的平方根只可能是1和1。因此1na对模n的平方根21na也只可能是1和1。如果21n本身还是一个偶数，我们可以再取一次平方根……将这些写成一个算法：（Rabin-Miller强伪素数测试）记dns21，其中d是奇数而s非负。如果nadmod1，或者对某个sr0有nadrmod12，则认为n通过测试，并称之为一个以a为基的强可能素数（strongprobableprimebasea，a-SPRP）。若n是合数，则又称之为一个以a为基的强伪素数（strongpseudoprime）。这个测试同时被冠以Miller的名字是因为Miller提出并证明了如下测试：如果扩展黎曼猜想（extendedRiemannhypothesis）成立，那么对于所有满足2log21na的基a，n都是a-SPRP，则n是素数。尽管Monier和Rabin在1980年证明了这个测试的错误概率（即合数通过测试的概率）不超过41，单个测试相对来说还是相当弱的（Pomerance、Selfridge和Wagstaff,Jr.证明了对任意1a都存在无穷多个a-SPRP）。但由于不存在“强Carmichael数”（任何合数n都存在一个基a试之不是a-SPRP），我们可以组合多个测试来产生有力的测试，以至于对足够小的n可以用来证明其是否素数。取前k个素数为基，并用k来表示以前k个素数为基的强伪素数，Riesel在1994年给出下表：321,972,326,370,024,942,526,193,897,56401,002,541,205,143,073,566,360,553,1041,689,705,135,316,234,41321,728,071,550,341321,728,071,550,341383,660,749,474,3747,898,302,152,2751,031,215,3001,326,25653,373,1047,21110987654321考虑到64位二进制数能表示的范围，只需取前9个素数为基，则对小于8的所有大于1的整数测试都是正确的；对大于或等于8并小于642的整数测试错误的概率不超过1821。Rabin-Miller强伪素数测试本身的形式稍有一些复杂，在实现时可以下面的functionpowermod(a,s,n){p:=1b:=awhiles0{if(s&1)==1thenp:=p*b%nb:=b*b%ns:=s1}returnp}functionmul64to128(a,b){{ah,al}:={a32,a&0xffffffff}{bh,bl}:={b32,b&0xffffffff}简单形式代替：对1n，如果nanmod11则认为n通过测试。代替的理由可简单证明如下：仍然记dns21，其中d是奇数而s非负。若n是素数，由nanmod11可以推出naansdmod12112或nadsmod112。若为前者，显然取1sr即可使n通过测试。若为后者，则继续取平方根，直到对某个sr1有nadrmod12，或nadmod12。无论nadmod1还是nadmod1，n都通过测试。Rabin-Miller强伪素数测试的核心是幂取模（即计算nasmod）。计算幂取模有以下的nOlog算法（以Sprache伪代码语言描述）：这个算法在32位计算机上实现的难点在于指令集和绝大部分编程语言的编译器都只提供了32位相乘结果为64位的整数乘法，浮点运算由于精度的问题不能应用于这里的乘法。唯一解决办法是模仿一些编译器内建的64位整数乘法来实现两个无符号64位相乘结果为128位的乘法。这个乘法可以将两个乘数分别分割成两个32位数来实现。为方便乘法之后的取模运算，运算结果应当用连续的128个二进制位来表示。以下是其伪代码：rl:=al*blc:=al*bhrh:=c32c:=c32rl:=rl+cifrlcthenrh++//处理进位c:=ah*blrh:=rh+(c32)c:=c32rl:=rl+cifrlcthenrh++//处理进位rh:=rh+ah*bhreturn{rh,rl}}functionmod64(rh,rl,n){rh:=rh%nrl:=rl%nr:=fmodl(rh*2**64,n)r:=r+rlifrnthenr:=r-n//处理进位r:=fmodl(r,n)returnr}乘法之后的取模运算可以用浮点运算快速完成。具体做法是乘积的高64位和低64位分别先对除数取模，然后再利用浮点单元合并取模。这里的浮点运算要求浮点单元以最高精度运算，计算前应先将浮点单元控制字中的精度控制位设置为64位精度。为保证精度，应当用80位浮点数实现此运算。伪代码如下：至此，Rabin-Miller强伪素数测试的核心已经全部实现。二、Pollard因数分解算法Pollard因数分解算法又称为PollardMonteCarlo因数分解算法。它的核心思想是：选取一个随机数a。再选一个随机数b。检查nba,gcd是否大于1。functionpollard_floyd(n){x:=x0y:=x0d:=1whiled==1{x:=f(x)y:=f(f(y))若大于1，nba,gcd就是n的一个因子。若不大于1，再选取随机数c，检查nbc,gcd和nac,gcd。如此继续，直到找到n的一个非平凡因子。它的主要实现方法是从某个初值nx10出发，通过一个适当的多项式xf进行迭代，产生一个伪随机迭代序列,,,,,,111321xffxfxxxx直到其对模n出现循环。多项式xf的选择直接影响着迭代序列的长度。经典的选择是选取axxf2，其中namod2,0。不选择0和2的原因是避免当序列中某一项nxkmod1时后续各项全部为1的情况。迭代序列出现循环的期望时间和期望长度都与n成正比。设n有一个非平凡因子p。由前面可知，迭代序列关于模p出现循环的期望时间和期望长度与p成正比。当循环出现时，设pxxkjmod，记nxxdkj,gcd，则d是p的倍数。当nd时，d就是n的一个非平凡因子。但是在分解成功之前，p是未知的，也就无从直接判断循环是否已经出现。仍然设迭代序列关于模p出现循环，pxxkjmod并设kj。由取模运算的性质可知nxfxfkjmod，即pxxkjmod11。故对任意0c，总有pxxckcjmod。记循环部分长度为t，若kt，则kt2，那么nxxkkmod2。因此，只要对迭代序列中的偶数项kx和对应的2kx计算nxxdkk,gcd2。只要nd，d就是n的一个非平凡因子。以上就是Pollard原来建议使用的Floyd循环判断算法。由此得到Pollard因数分解算法的第一个伪代码：d:=gcd(x–y,n)if1dANDdnthenreturndifd==nthenreturnFAILURE}}functionpollard_brent(n){x:=x0y:=x0k:=0i:=1d:=1whiled==1{k:=k+1x:=f(x)d:=gcd(x–y,n)if1dANDdnthenreturndifd==nthenreturnFAILUREifk==ithen{y:=xi:=i1}}}Floyd循环判断算法对储存整个迭代序列的空间要求很高，一般实现时都使用时间换空间的办法，同时计算kx和kx2来进行判断（如上面的伪代码）。Brent提出了另外一种效率更高的循环判断算法：每步只计算kx，当k是2的方幂时，令kxy；每一步都拿当前的kx和y计算nyxdk,gcd。伪代码如下：由于循环出现的时空期望都是pO，Pollard因数分解算法的总体时间复杂度也就是pO。在最坏情况下，其时间复杂度可能接近nO，但在一般条件下，时间复杂度都可以认为是4nO。参考资料：ChrisCaldwell“Fermat,probable-primalityandpseudoprimes.”ChrisCaldwell“Strongprobable-primalityandapracticaltest.”EricW.Weisstein“Brent’sFactorizationMethod.”EricW.Weisstein“PollardRhoFactorizationMethod.”EricW.Weisstein“Rabin-MillerStrongPseudoprimeTest.”EricW.Weisstein“StrongPseudoprime.”UnknownAuthor“Pollard’s”