POJ 2191 Mersenne Composite Numbers 整数分解

题意：分解2^K次方以内的梅森复合数。

#include<cstdio>
#include<cstring>
#include<ctime>
#include<cstdlib>
#include<cmath>
#include<algorithm>
using namespace std;
#define lint __int64

lint ans[65];
int p[65], cnt;

lint gcd ( lint a, lint b )
{
    while ( b != 0 )
    {
        lint c = a % b;
        a = b;
        b = c;
    }
    return a;
}

/*
lint mod_mult ( lint a, lint b, lint n )
{
    lint s=0;
    a = a % n;
    for ( ; b; b >>= 1 )
    {
        if ( b & 1 ) s += a;
        if( s >= n) s -= n;
        a = a << 1;
        if ( a >= n ) a -= n;
    }
    return s;
}
*/


lint mod_mult ( lint a, lint b, lint n )
{
    lint ret = 0;
    for ( ; b; b >>= 1, a = (a+a) % n )
        if ( b & 1 )
            ret = (ret + a) % n;
    return ret;
}

lint mod_exp ( lint a, lint b, lint n )
{
    lint ret = 1;
    for ( ; b; b >>= 1, a = mod_mult(a,a,n) ) // a = (a*a)%n;
        if ( b & 1 )
            ret = mod_mult(ret,a,n); // ret = (ret*a)%n;
    return ret;
}

//witness只能检测奇数
//判断奇数n是否是合数，若是则返回true,不是返回false
bool witness ( lint a, lint n )
{
    int i, t = 0;
    lint m = n-1, x, y;
    while ( m % 2 == 0 ) { m >>= 1; t++; }
    x = mod_exp ( a, m, n );

    for ( i = 1; i <= t; i++ )
    {
        y = mod_exp(x,2,n);
        if ( y==1 && x!=1 && x!=n-1 )
            return true;
        x = y;
    }
    if ( y != 1 ) return true;
    return false;
}


bool miller_rabin1 ( lint n, int times = 10 )
{
    if ( n==1 || (n!=2&&!(n%2)) || (n!=3&&!(n%3)) || (n!=5&&!(n%5)) || (n!=7&&!(n%7)) )
        return false;
    if ( n == 2 ) return true; //由于witness只检测奇数，偶素数2需要特殊处理

    while ( times-- )
    {
        lint a = rand()%(n-1)+1;
        if ( witness(a, n) ) return false;
    }
    return true;
}

bool miller_rabin2 ( lint n, int times = 10 )
{
    if ( n==1 || (n!=2&&!(n%2)) || (n!=3&&!(n%3)) || (n!=5&&!(n%5)) || (n!=7&&!(n%7)) )
        return false;

    while ( times-- )
    {
        lint m = mod_exp ( rand()%(n-1) + 1, n-1, n );
        if ( m != 1 ) return 0;
    }
    return true;
}

lint rho ( lint n )
{
    lint i, k, x, y, d;
    i = 1; k = 2;
    srand ( time(NULL) );
    y = x = rand() % n;
    lint c =  rand() % (n-1) + 1;
    while ( true )
    {
        i++;
        x = (mod_mult(x,x,n)+c) % n;
        d = gcd ( y - x, n );
        if ( d > 1 && d < n ) return d;
        if ( x == y ) break;
        if ( i == k ) { y = x; k = k*2; }
    }
    return n;
}

void pollard ( lint n )
{
    if ( n == 1 ) return; //1不加入因式分解
    if ( miller_rabin1(n) ) { ans[cnt++] = n; return; }

    lint m = n;
    while ( m >= n )
        m = rho ( n );
    pollard ( m );
    pollard ( n/m );
}


void prime ()
{
    int pn = 0;
    bool f[100] = { 0 };
    for ( int i = 2; i < 100; i++ )
    {
        if ( f[i] ) continue;
        p[pn++] = i;
        for ( int j = 2; i * j < 100; j++ )
            f[i*j] = 1;
    }
}

int main()
{
    prime();
    int k,i,j;
    scanf("%d",&k);
    for ( i = 1; p[i] <= k; i++ )
    {
            lint n=((lint)1<<p[i])-1;
            if ( miller_rabin1(n) ) continue;
            memset ( ans, 0, sizeof(ans) );
            cnt = 0;
            pollard ( n );
            sort ( ans, ans + cnt );
            for ( j = 0; j < cnt-1; j++ )
                printf("%I64d * ", ans[j]);
            printf("%I64d = %I64d = ( 2 ^ %d ) - 1\n",ans[j],n,p[i]);
    }
    return 0;
}