poj 2429 GCD & LCM Inverse

一个对Miller_rabin与pallord的一个运用；

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<queue>
#include<set>
#include<map>
#include<cstring>
#include<vector>
#define LL unsigned long long
using namespace std;
LL p[10] = { 2,3,5,7,11,13,17,19,23,29 },num[100],ans,sq;
int cnt;
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;
}
LL 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 < 10; 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;
   LL i = 1,k=2,d;
   x = y = abs( 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 )
{
    if( n == 1 ) return ;
    if( Miller_rabin( n ) )
    {
       num[cnt++] = n;
       return ;    
    }    
    LL p = Pallord( n );
    Pallord_Min( p );
    Pallord_Min( n/p );
}
void Solve( int n ,LL val )
{
    if( n >= cnt )
    {
        if( val > ans && val <= sq )
            ans = val;
        return ;    
    }
    Solve( n + 1, val*num[n] );
    Solve( n + 1 , val );
}
bool cmp( LL a ,LL b )
{
    return a < b;    
}
int main(  )
{
    LL n,m;
    while( scanf( "%I64u %I64u",&n,&m )==2 )
    {
        if( n > m ) swap( n , m );
        cnt = 0;
        m /= n;
        Pallord_Min( m );
        sort( num , num + cnt ,cmp);
        int j = 0;
        for ( int i = 1; i < cnt; i++ )  
        {  
        while ( num[i-1] == num[i] && i < cnt )  
            num[j] *= num[i++];  
        if ( i < cnt ) num[++j] = num[i];  
        }  
        ans = 1;
        sq = (LL)sqrt(m*1.0);
        cnt = j + 1;
        Solve( 0 ,1 );    
        printf( "%I64u %I64u\n",ans*n,(m/ans) * n );
    }
    //system( "pause" );
    return 0;
}