Poj Solution 2447

http://poj.org/problem?id=2447

#include<stdlib.h>
#include<stdio.h>
#include<string.h>
#include<algorithm>
using namespace std;
typedef unsigned __int64 u64;
#define MAX 30
#define MAXN 5
u64 len, dig, limit;
u64 factor[MAXN];
u64 mod(u64 a, u64 b, u64 n){
        if(!a)return 0;
        else return ( ((a&dig)*b)%n + (mod(a>>len,b,n)<<len)%n )%n;
}
u64 by(u64 a, u64 b, u64 n){
            u64 p;
            p = 8, len = 61;
        while(p<n){
                    p<<=4;
                    len -=4;
        }
            dig = ((limit/p)<<1) - 1; 
        return mod(a,b,n);
}
u64 random(){//产生随机数 
            u64 a;
            a = rand();
            a *= rand();
            a *= rand();
            a *= rand();
        return a;
}
u64 square_multiply(u64 x, u64 c, u64 n){//(x^c)%n快速算法 
            u64 z=1;
        while(c){
                if(c%2==1)z = by(z,x,n);
                    x = by(x,x,n);
                    c=(c>>1);
        }
        return z;
}
bool Miller_Rabin(u64 n){//Miller_Rabin素数测试 
    if(n<2)return false;
        if(n==2)return true;
        if(!(n&1))return false;
            u64 k = 0, i, j, m, a;
            m = n - 1;
        while(m%2==0)m=(m>>1),k++;
        for(i=0;i<MAX;i++){
                    a = square_multiply(random()%(n-1)+1, m, n);//平方乘
                if(a==1)continue;
                for(j=0;j<k;j++){
                        if(a==n-1)break;
                            a = by(a,a,n);
                }
                if(j<k)continue;
                return false ;
        }
        return true;
}
u64 gcd(u64 a,u64 b){
        if(a==0) return b;
        else return gcd(b%a,a);
}
u64 f(u64 x, u64 n ){
        return (by(x,x,n)+1)%n;
}
u64 Pollard(u64 n){
                if(n<=2)return 0;
                if(!(n&1))return 2; 
            u64 i, p, x,xx;
        for(i=0;i<MAX;i++){                                 
                    x = random()%n; //或者直接用 x = i
                    xx = f(x,n);
                    p = gcd((xx+n-x)%n , n);
                while(p==1){
                            x = f(x,n);
                            xx = f(f(xx,n),n);
                            p = gcd((xx+n-x)%n,n)%n;
                }
                if(p)return p;
        }
        return 0;
}
u64 prime(u64 a){
        if(Miller_Rabin(a))return a;
            u64 t = Pollard(a);
        if(t!=0)
        {return prime(t);}
           
}
u64 Euclid(u64 a,u64 b,__int64 &x,__int64 &y)
{
 if(b==0)
 { 
     x=1,y=0;
  return a;
 } 
 u64 t,d;
 if(b!=0)
 d=Euclid(b,a%b,x,y);
 t=x;
 x=y;
 if(b!=0)
 y=t-a/b*y;
 return d;
}
int main(){
            u64 c,e,n,i,j,m,t,n0,T,ans,l;
        __int64 T0,x,y,d;
            limit = 1;
            limit = limit<<63; 
        while(scanf("%I64u%I64u%I64u",&c,&e,&n)!=EOF){
              m=0;n0=n;
               while(n%2==0){n/=2;factor[m++]=2;}
               while(n>1){
                        if(Miller_Rabin(n))break;
                            t = prime(n);
                            factor[m++] = t;
                        if(t!=0)
                            n/=t;
                }
               if(n>1)factor[m++]=n;
               //for(l=0;l<m;l++)printf("%I64un",factor[l]);
                   T0=(__int64)(factor[0]-1)*(factor[1]-1);
          T=(u64)T0;
               Euclid(e,T,x,y);
                   d=x;
      //printf("%I64dn",d);
      //while(d<0)d+=T0;
                   d=(d%T0+T0)%T0;
               //d=(__int64)d;
     // printf("%I64d %I64dn",d,T0);
                   ans=square_multiply(c,(u64)d,n0);
               printf("%I64un",ans);
                   
}
}