Pollard_Rho大数分解模板题 pku-2191

题意：给你一个数n，  定义m=2^k-1,   {k|1<=k<=n}，并且 k为素数;  当m为合数时，求分解为质因数,输出格式如下:47 * 178481 = 8388607 = ( 2 ^ 23 ) - 1  

分析:要分解m，首先要判断m是否为合数，直接用米勒拉宾判断，但是后面的大合数分解，一开始用了试除法，直接给超时。所以，有更加快速的方法。（现学的）。使用Pollard_Rho大数分解算法。  

Pollard_Rho大数分解时间复杂度为n^1/4 

ac代码：

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

typedef long long llt;

int const Repeat = 10;
const int N = 105;
int prime[N];
bool vis[N], is_prime[N];
//利用二进制计算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 Prime()
{
    int cnt = 0;
    for (int i = 2; i <= N; ++i)
    {
        if (!vis[i])
        {
            prime[cnt++] = i;    is_prime[i] = 1;
        }
        for (int j = 0; j < cnt&&i*prime[j] <= N; ++j)
        {
            vis[i*prime[j]] = 1;
            if (i%prime[j] == 0)break;
        }
    }
    return cnt;
}

llt poww(llt x, llt n)
{
    llt res = 1;
    for (; n;x=x*x, n>>=1)
    if (n & 1)res = res*x;
    return res;
}
int main()
{
    int k = Prime();
    llt n;
    while (scanf("%lld", &n) != EOF)
    {
        for (int i = 0; prime[i] <= n; ++i)
        {
            llt pnum = poww(2, prime[i]) - 1;
            if (!Miller_Rabin(pnum, 5))
            {
                FCnt = 0;
                find(pnum);
                sort(Fac, Fac + FCnt);
                printf("%lld", Fac[0]);
                for (int i = 1; i < FCnt; ++i)
                    printf(" * %lld", Fac[i]);
                printf(" = %lld = ( 2 ^ %d ) - 1
", pnum, prime[i]);
            }
        }
    }
}