Pollard rho模板

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<cstdlib>
 4 #include<algorithm>
 5 #include<ctime>
 6 #include<cmath>
 7 #include<iostream>
 8 using namespace std;
 9 
10 #define LL long long
11 LL n;
12 #define maxs 80
13 LL fac[maxs],num[maxs],lf=0;
14 LL mul(LL a,LL b,LL p) 
15 {
16     LL tmp=a,ans=0;
17     while (b) {if (b&1) ans=(ans+tmp)%p;tmp=(tmp+tmp)%p;b>>=1;}
18     return ans;
19 }
20 LL random(LL x) {return (LL)((double)rand()/RAND_MAX*(x-1)+0.5);}
21 LL f(LL x,LL c,LL p) {return (mul(x,x,p)+c)%p;}
22 LL gcd(LL a,LL b) {if (a<b) return gcd(b,a);return b?gcd(b,a%b):a;}
23 LL play(LL x,LL c)
24 {
25     LL a=random(x-1)+1,b=f(a,c,x),i=1,k=2;
26     if (a==b) return x;
27     while (++i)
28     {
29         LL sig=gcd(fabs(a-b),x);
30         if (sig>1 && sig<x) return sig;
31         b=f(b,c,x);
32         if (a==b) return x;
33         if (i==k) a=b,k<<=1;
34     }
35 }
36 LL pow_mod(LL a,LL b,LL p)
37 {
38     LL tmp=a,ans=1;
39     while (b) {if (b&1) ans=mul(ans,tmp,p);tmp=mul(tmp,tmp,p);b>>=1;}
40     return ans;
41 }
42 bool is_prime(LL x)
43 {
44     LL num=x-1,k=0;while (!(num&1)) num>>=1,k++;
45     for (int i=1;i<=20;i++)
46     {
47         LL t=random(x-2)+1,jud=pow_mod(t,num,x),now=jud;
48         for (int j=1;j<=k;j++)
49         {
50             jud=mul(jud,jud,x);
51             if (jud==1 && now!=1 && now!=x-1) return 0;
52             now=jud;
53         }
54         if (jud!=1) return 0;
55     }
56     return 1;
57 }
58 void rho(LL x,LL c)
59 {
60     if (x==1) return;
61     if (is_prime(x)) {fac[++lf]=x;return;}
62     LL p=x,tmp=c;
63     while (p==x) p=play(x,tmp--);
64     rho(p,c);rho(x/p,c);
65 }
66 int main()
67 {
68     srand(time(NULL));
69     scanf("%lld",&n);
70     memset(fac,0,sizeof(fac));
71     rho(n,120);sort(fac+1,fac+1+lf);num[1]=1;LL nlf=1;
72     for (int i=2;i<=lf;i++)
73         if (fac[i-1]==fac[i]) num[nlf]++;
74         else num[++nlf]=1,fac[nlf]=fac[i];
75     lf=nlf;
76     for (int i=1;i<lf;i++) printf("%lld^%lld*",fac[i],num[i]);printf("%lld^%lld\n",fac[lf],num[lf]);
77     return 0;
78 }