POJ 2191 Mersenne Composite Numbers 整数分解

#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;
}