# 【模板】SPOJ FACT0 大数分解 miller-rabin & pollard-rho

http://www.spoj.com/problems/FACT0/en/

miller-rabin & pollard-rho模板

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int xjb=10;
ll mmul(ll a, ll b, ll m){
ll d=((long double)a/m*b+1e-8);
ll r=a*b-d*m;
return r<0?r+m:r;
}
ll mpow(ll a, ll b, ll m){ll r=1;for(;b;b>>=1,a=mmul(a,a,m))if(b&1)r=mmul(r,a,m);return r;}
ll gcd(ll a, ll b){return a?gcd(b%a,a):b;}
int prime(ll n){
if(n==1) return 0;
if(n==2||n==3||n==5) return 1;
if(!(n&1)||(n%3==0)||(n%5==0)) return 0;
ll m=n-1; int k=0;
while(!(m&1)) m>>=1, k++;
for(int tt=0; tt<xjb; ++tt){
ll x=mpow(rand()%(n-2)+2,m,n), y=x;
for(int i=0; i<k; ++i){
x=mmul(x,x,n);
if(x==1&&y!=1&&y!=n-1) return 0;
y=x;
}
if(x!=1) return 0;
}
return 1;
}
ll f[105]; int M;
ll rho(ll n, ll c){
ll x=rand()%n, y=x, t=1;
for(int i=1, k=2; t==1; ++i){
x=(mmul(x,x,n)+c)%n;
t=gcd(x>y?x-y:y-x, n);
if(i==k) y=x, k<<=1;
}
return t;
}
void work(ll n){
if(n==1) return;
if(prime(n)){f[M++]=n; return;}
ll t=n;
while(t==n) t=rho(n, rand()%5+1);
work(t); work(n/t);
}
int main(){
srand(19260817);
ll n;
while(scanf("%lld", &n), n){
if(n==1){puts(""); continue;}
M=0;
work(n);
sort(f, f+M);
for(int i=0, c=1; i<M; ++i){
if(f[i]!=f[i+1]) printf("%lld^%d ", f[i], c), c=1;
else c++;
}
puts("");
}
return 0;
}

posted @ 2017-03-20 15:09  will7101  阅读(348)  评论(0编辑  收藏  举报