64位以内Rabin-Miller 强伪素数测试和Pollard rho 因数分解解析

在求解POJ1811题Prime Test中应用到的两个重要算法是Rabin-Miller强伪素数测试和Pollard r因数分解算法。前者可以在 的时间内以很高的成功概率判断一个整数是否是素数。后者可以在最优 的时间内完成合数的因数分解。这两种算法相对于试除法都显得比较复杂。本文试图对这两者进行简单的阐述，说明它们在32位计算机上限制在64位以内的条件下的实现中的细节。下文提到的所有字母均表示整数。

 

一、Rabin-Miller强伪素数测试

Rabin-Miller强伪素数测试的基本思想来源于如下的Fermat小定理：

如果p是一个素数，则对任意a有 。特别的，如果p不能整除a，则还有 。

利用Fermat小定理可以得到一个测试合数的有力算法：对 ，选择 ，计算 ，若结果不等于1则n是合数。若结果等于1则n可能是素数，并被称为一个以a为基的弱可能素数（weak probable prime base a，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和 。因此 对模n的平方根 也只可能是1和 。如果 本身还是一个偶数，我们可以再取一次平方根……将这些写成一个算法：

（Rabin-Miller强伪素数测试）记 ，其中d是奇数而s非负。如果 ，或者对某个 有 ，则认为n通过测试，并称之为一个以a为基的强可能素数（strong probable prime base a，a-SPRP）。若n是合数，则又称之为一个以a为基的强伪素数（strong pseudoprime）。

这个测试同时被冠以Miller的名字是因为Miller提出并证明了如下测试：如果扩展黎曼猜想（extended Riemann hypothesis）成立，那么对于所有满足 的基a，n都是a-SPRP，则n是素数。

尽管Monier和Rabin在1980年证明了这个测试的错误概率（即合数通过测试的概率）不超过 ，单个测试相对来说还是相当弱的（Pomerance、Selfridge和Wagstaff, Jr.证明了对任意 都存在无穷多个a-SPRP）。但由于不存在“强Carmichael数”（任何合数n都存在一个基a试之不是a-SPRP），我们可以组合多个测试来产生有力的测试，以至于对足够小的n可以用来证明其是否素数。

取前k个素数为基，并用 来表示以前k个素数为基的强伪素数，Riesel在1994年给出下表：

考虑到64位二进制数能表示的范围，只需取前9个素数为基，则对小于 的所有大于1的整数测试都是正确的；对大于或等于 并小于 的整数测试错误的概率不超过 。

Rabin-Miller强伪素数测试本身的形式稍有一些复杂，在实现时可以下面的简单形式代替：

对 ，如果 则认为n通过测试。

代替的理由可简单证明如下：

仍然记 ，其中d是奇数而s非负。若n是素数，由 可以推出 或 。若为前者，显然取 即可使n通过测试。若为后者，则继续取平方根，直到对某个 有 ，或 。无论 还是 ，n都通过测试。

function powermod(a, s, n)

{

   p := 1

   b := a

   while s > 0

   {

      if (s & 1) == 1 then p := p * b % n

      b := b * b % n

      s := s >> 1

   }

   return p

}

Rabin-Miller强伪素数测试的核心是幂取模（即计算 ）。计算幂取模有以下的 算法（以Sprache伪代码语言描述）：

function mul64to128(a, b)

{

   {ah, al} := {a >> 32, a & 0xffffffff}

   {bh, bl} := {b >> 32, b & 0xffffffff}

这个算法在32位计算机上实现的难点在于指令集和绝大部分编程语言的编译器都只提供了32位相乘结果为64位的整数乘法，浮点运算由于精度的问题不能应用于这里的乘法。唯一解决办法是模仿一些编译器内建的64位整数乘法来实现两个无符号64位相乘结果为128位的乘法。这个乘法可以将两个乘数分别分割成两个32位数来实现。为方便乘法之后的取模运算，运算结果应当用连续的128个二进制位来表示。以下是其伪代码：

   rl := al * bl

   c := al * bh

   rh := c >> 32

   c := c << 32

   rl := rl + c

   if rl < c then rh++   // 处理进位

   c := ah * bl

   rh := rh + (c >> 32)

   c := c << 32

   rl := rl + c

   if rl < c then rh++   // 处理进位

   rh := rh + ah * bh

   return {rh, rl}

}

乘法之后的取模运算可以用浮点运算快速完成。具体做法是乘积的高64位和低64位分别先对除数取模，然后再利用浮点单元合并取模。这里的浮点运算要求浮点单元以最高精度运算，计算前应先将浮点单元控制字中的精度控制位设置为64位精度。为保证精度，应当用80位浮点数实现此运算。伪代码如下：

function mod64(rh, rl, n)

{

   rh := rh % n

   rl := rl % n

   r := fmodl(rh * 2 ** 64, n)

   r := r + rl

   if r < n then r := r - n    // 处理进位

   r := fmodl(r, n)

   return r

}

至此，Rabin-Miller强伪素数测试的核心已经全部实现。

 

二、Pollard r因数分解算法

Pollard r因数分解算法又称为Pollard Monte Carlo因数分解算法。它的核心思想是：选取一个随机数a。再选一个随机数b。检查 是否大于1。若大于1， 就是n的一个因子。若不大于1，再选取随机数c，检查 和 。如此继续，直到找到n的一个非平凡因子。

它的主要实现方法是从某个初值 出发，通过一个适当的多项式 进行迭代，产生一个伪随机迭代序列 直到其对模n出现循环。多项式 的选择直接影响着迭代序列的长度。经典的选择是选取 ，其中 。不选择0和 的原因是避免当序列中某一项 时后续各项全部为1的情况。迭代序列出现循环的期望时间和期望长度都与 成正比。设n有一个非平凡因子p。由前面可知，迭代序列关于模p出现循环的期望时间和期望长度与 成正比。当循环出现时，设 ，记 ，则d是p的倍数。当 时，d就是n的一个非平凡因子。

function pollard_floyd(n)

{

   x := x0

   y := x0

   d := 1

   while d == 1

   {

      x := f(x)

      y := f(f(y))

但是在分解成功之前，p是未知的，也就无从直接判断循环是否已经出现。仍然设迭代序列关于模p出现循环， 并设 。由取模运算的性质可知 ，即 。故对任意 ，总有 。记循环部分长度为t，若 ，则 ，那么 。因此，只要对迭代序列中的偶数项 和对应的 计算 。只要 ，d就是n的一个非平凡因子。以上就是Pollard原来建议使用的Floyd循环判断算法。由此得到Pollard r因数分解算法的第一个伪代码：

function pollard_brent(n)

{

   x := x0

   y := x0

   k := 0

   i := 1

   d := 1

   while d == 1

   {

      k := k + 1

      x := f(x)

      d := gcd(x – y, n)

      if 1 < d AND d < n then return d

      if d == n then return FAILURE

      if k == i then

      {

        y := x

        i := i << 1

      }

   }

}

      d := gcd(x – y, n)

      if 1 < d AND d < n then return d

      if d == n then return FAILURE

   }

}

由于循环出现的时空期望都是 ，Pollard r因数分解算法的总体时间复杂度也就是 。在最坏情况下，其时间复杂度可能接近 ，但在一般条件下，时间复杂度都可以认为是 。

 

参考资料：

Chris Caldwell “Fermat, probable-primality and pseudoprimes.”

Chris Caldwell “Strong probable-primality and a practical test.”

Eric W. Weisstein “Brent’s Factorization Method.”

Eric W. Weisstein “Pollard Rho Factorization Method.”

Eric W. Weisstein “Rabin-Miller Strong Pseudoprime Test.”

Eric W. Weisstein “Strong Pseudoprime.”

Unknown Author “Pollard’s r.”

下面附一个米勒罗宾算法和布莱德算法的代码给大家当模板：

#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#include <functional>
const int MAX_COUNT = 6;
unsigned __int64 allFactor[65], nFactor;

unsigned __int64 fMultiModular( unsigned __int64 a, unsigned __int64 b, const unsigned __int64 n)
{
    unsigned __int64 back = 0, temp = a % n;
    while ( b > 0 )
    {
        if ( b & 0x1 )
        {
            if ( (back = back + temp) > n )
                back -= n;
        }
        if ( (temp <<= 1) > n )
            temp -= n;
        b >>= 1;
    }
    return back;
}

unsigned __int64 GCD(unsigned __int64 a, unsigned __int64 b)
{
    if(b==0)return a;
    else return GCD(b,a%b);
}
unsigned __int64 modular_exp(const unsigned __int64 &a, unsigned __int64 b, const unsigned __int64 &n)
{
    unsigned __int64 d(1), dTemp(a % n);
    while ( b > 0 )
    {
        if ( b & 0x1 )
            d = fMultiModular(d, dTemp, n);
        dTemp = fMultiModular(dTemp, dTemp, n);
        b >>= 1;
    }
    return d;
}
bool fWitNess(const unsigned __int64 &a, const unsigned __int64 &n, char t, const unsigned __int64 &u)
{
    unsigned __int64 x, y;
    x = modular_exp(a, u, n);
    while ( t-- )
    {
        y = fMultiModular(x, x, n);
        if ( y == 1 && x != 1 && x != n-1 )
            return true; //must not
        x = y;
    }
    return y != 1;
}
bool miller_rabin(const unsigned __int64 &n, int count)
{
    if ( n == 1 )
        return false;
    if ( n == 2 )
        return true;
    unsigned __int64 a, u;
    char t(0);
    for (u = n-1; !(u & 0x1); u>>=1)++t;
    for (int i = 1; i <= count; ++i)
    {
        a = rand() % (n-1) + 1;
        if ( fWitNess(a, n, t, u) )
            return false;
    }
    return true;
}
unsigned __int64 pollard_rho(const unsigned __int64 &c, const unsigned __int64 &num)
{
    int i(1), k(2);
    unsigned __int64 x = rand() % num;
    unsigned __int64 y = x, comDiv;
    do
    {
        ++i;
        if ( (x = fMultiModular(x, x, num) - c) < 0 )
            x += num;
        if ( x == y )
            break;
        comDiv = GCD((y-x+num)%num, num);
        if ( comDiv > 1 && comDiv < num )
            return comDiv;
        if ( i == k )
        {
            y = x;
            k <<= 1;
        }
    }
    while ( true );
    return num;
}
void fFindFactor(const unsigned __int64 &num)
{
    if ( miller_rabin(num, MAX_COUNT) == true )
    {
        allFactor[++nFactor] = num;
        return;
    }
    unsigned __int64 factor;
    do
    {
        factor = pollard_rho(rand()%(num-1) + 1, num);
    }
    while ( factor >= num );

    fFindFactor(factor);
    fFindFactor(num/factor);
}

int main()
{
    int t,i;
    unsigned __int64 test_num,min ,factor;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%I64u",&test_num);
        if(miller_rabin(test_num,5))
            printf("Prime\n");
        else
        {
            min = test_num;
            nFactor = 0;
            fFindFactor(test_num);
            for(i = 1; i <= nFactor ; i++)
            {
                if(allFactor[i] < min)
                    min = allFactor[i];
            }
            printf("%I64u\n",min);
        }
    }
}