您的位置:首页 > 其它

poj 2191 Mersenne Composite Numbers

2012-07-17 20:47 495 查看
首先利用miller_rabin测试是否为素数;

再利用pallord进行质因子分解;

View Code

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#define LL  unsigned long long
using namespace std;
LL p[10] = { 2,3,5,7,11,13,17,19,23,29 };
int cnt = 0; LL num[100];
LL Multi( LL a , LL b, LL n )
{
LL sum = 0;
while( b )
{
if( ( b&1 ) ) sum = ( sum + a )%n;
a <<= 1;
b >>= 1;
if( a >= n ) a %=n;
}
return sum;
}
LL Quick_mod( LL a, LL b, LL n )
{
LL sum = 1;
while( b )
{
if( ( b&1 ) ) sum = Multi( sum , a , n );
a = Multi( a , a ,n );
b >>= 1;
}
return sum;
}
bool Miller_rabin( LL n )
{
LL m = n -1,d,L;
int k = 0;
if( n ==2 ) return true;
if( n < 2 || !( n & 1 ) ) return false;
while( !( m &1 ) )
{
k ++ ; m >>=1;
}
for( int i = 0 ; i < 9 ; i ++ )
{
if( p[i] >= n ) return true;
d = Quick_mod( p[i] , m , n );
if( d==1 ) continue;
for( int j = 0 ; j < k ; j ++ )
{
L = Multi( d , d ,n );
if( L == 1 && d != 1 && d !=( n -1 ) )
return false;
d = L;
}
if( d != 1 ) return false;
}
return true;
}
LL Gcd( LL a ,LL b )
{
return b == 0 ? a : Gcd( b, a%b );
}
LL Pallord( LL n )
{
LL x,y,c = 0,d,i=1,k=2;
x = y = rand()%( n -1 )     + 1;
while( c ==0 || c == 2 )
{
c = abs(rand())%( n -1 ) + 1;
}
do
{
i ++;
d = Gcd( y + n - x ,n );
if( d > 1 && d < n )
return d;
if( i == k )
{
y = x ; k <<= 1;
}
x = (Multi( x , x , n ) + c)%n;
}while( x != y );
return n;
}
void Pallord_Min( LL n )
{
LL a,b,p;
if( n == 1 ) return ;
if( Miller_rabin( n ) )
{
num[cnt++] = n;
//        printf( "safas%I64u\n",n );
return ;
}
p = Pallord( n );
Pallord_Min( p );
Pallord_Min( n / p );
//    return a < b ? a:b;
}
int main( )
{
LL k = 2;
int m;
while( scanf( "%d",&m )==1 )
{
k = 2;
for( int i = 2 ; i <= m ; i ++ )
{
k *= 2;
if( Miller_rabin( (LL)i ) && Miller_rabin( k - 1 ) == false)
{
cnt = 0;
Pallord_Min( k -1 );
sort( num , num + cnt );
for( int j = 0; j < cnt-1; j ++ )
printf( "%I64u * ",num[j] );
printf( "%I64u = %I64u = ( 2 ^ %d ) - 1\n",num[cnt-1],k-1,i );
}
}
}
//    system( "pause" );
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: