Pollard-Rho方法分解整数因子

对于大整数而言，要找出其因子是比较困难的，事实上仍然没有很好的解决方法。不过在ACM里面有几种“标准的”方法可以使用，例如Pollard-Rho方法。
假设题目是求$n$的因子，Pollard-Rho方法的基本思想是任意给定整数$y$和$x$，求$gcd(y-x,n)$。注意，所有运算都是在模$n$的条件下进行的。如果$gcd\neq1且gcd\neq{n}$，则$gcd$肯定是$n$的一个因子。无论有没有找到因子，均选取另外一对数$y$和$x$继续。注意到在模$n$的意义下，数对是有限的。可以期待这个过程能够输出将n分解的所需的足够的因子。
一般做法是生成一个数列$\lbrace{x_n}\rbrace$

$$x_{i+1}=(x_{i}^2-c){\%}n$$
$c$取一个随机数，但不要取值0或者2。$y$也无需另外取值，直接取序列 $\lbrace{x_1},x_2,x_4,x_8...\rbrace$即可。
注意到如果找到 $p$为$n$的因子，那么 $n/p$也是 $n$的因子。且$p$和 $n/p$各自的因子也一定是 $n$的因子，而且$p$和 $n/p$一定包含了 n <script type="math/tex" id="MathJax-Element-1014">n</script>的所有可能的因子。因此可以进行递归调用。

POJ1811要求判断给定数是否为质数，如果是合数则输出最小的质因子。先用Miller_Rabin方法判质数，然后对合数递归分解因子求出最小即可。

#include <cstdio>
#include <algorithm>
using namespace std;

typedef long long llt;

int const Repeat = 10;

//利用二进制计算a*b%mod
llt multiMod(llt a,llt b,llt mod){
    llt ret = 0LL;
    a %= mod;
    while( b ){
        if ( b & 1LL ) ret = ( ret + a ) % mod, --b;
        b >>= 1LL;
        a = ( a + a ) % mod;
    }
    return ret;
}

//计算a^b%mod
llt powerMod(llt a,llt b,llt mod){
    llt ret = 1LL;
    a %= mod;
    while( b ){
        if ( b & 1LL ) ret = multiMod(ret,a,mod),--b;
        b >>= 1LL;
        a = multiMod(a,a,mod);
    }
    return ret;
}

//Miller-Rabin测试,测试n是否为素数
bool Miller_Rabin(llt n,int repeat){
    if ( 2LL == n || 3LL == n ) return true;
    if ( !( n & 1LL ) ) return false;

    //将n分解为2^s*d
    llt d = n - 1LL;
    int s = 0;
    while( !( d & 1LL ) ) ++s, d>>=1LL;

    //srand((unsigned)time(0));
    for(int i=0;i<repeat;++i){//重复repeat次
        llt a = rand() % ( n - 3 ) + 2;//取一个随机数,[2,n-1)
        llt x = powerMod(a,d,n);
        llt y = 0LL;
        for(int j=0;j<s;++j){
            y = multiMod(x,x,n);
            if ( 1LL == y && 1LL != x && n-1LL != x ) return false;
            x = y;
        }
        if ( 1LL != y ) return false;
    }
    return true;
}

llt Fac[100];//质因数分解结果（刚返回时是无序的）
int FCnt;//质因数的个数。数组小标从0开始

llt gcd(llt a,llt b){
    if (0L == a || 0L == b) return 1;
    if ( a < 0 ) a = -a;
    if ( b < 0 ) b = -b;
    while(b){
        llt t = a % b;
        a = b;
        b = t;
    }
    return a;
}
llt Pollard_Rho(llt n,llt c){
    llt i = 1, k = 2;
    llt x = rand() % n;
    llt y = x;
    while(1){
        ++i;
        x = ( multiMod(x,x,n) + c ) % n;
        llt d = gcd(y-x,n);
        if (d!=1LL && d!=n) return d;
        if (y == x) return n;
        if ( i == k ) y = x, k <<= 1;
    }
}

void find(llt n){
    if ( 4LL == n ){
        Fac[0] = Fac[1] = 2LL;
        FCnt = 2;
        return;
    }
    if ( Miller_Rabin(n,Repeat) ){
        Fac[FCnt++] = n;
        return;
    }

    llt p;
    while( ( p = Pollard_Rho(n,rand()%(n-3)+3) ) == n );

    find(p);
    find(n/p);
}

int main(){
    int kase;
    scanf("%d",&kase);
    while(kase--){
        llt n;
        scanf("%I64d",&n);

        FCnt=0;
        find(n);


        if(Fac[0] == n) printf("Prime\n");
        else printf("%I64d\n",*min_element(Fac,Fac+FCnt));
    }
    return 0;
}