poj 2191 Mersenne Composite Numbers

首先利用miller_rabin测试是否为素数；

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

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