BZOJ3667: Rabin-Miller算法 （Miller-Rabinpol_rho特技快速乘学习笔记）

学习了两个新算法，判素数和分解大合数。

Miller-Rabin(以下简称MR)算的时候大概就是这样子算。

那么怎么分解质因数呢？

这两篇博文应该已经说得比较清楚了。
http://www.cnblogs.com/thythy/p/5493624.html http://blog.csdn.net/maxichu/article/details/45459533

rho每次只找一个因数，要找n的因数，每次随机一个x出来，首先令y=x，然后令y不断做如下变化y=y*y+c(c是自己定的一个随机数)，然后用abs(x-y)去跟n求gcd，如果不为1，则找到一个，有一个步长扩二倍的优化（详见代码），不太明白为啥，估计是为了防止死循环。

找到一个因数之后递归求解（配合MR），找到所有质因数，期望复杂度$O(n^{1/4})$(口胡)

//QWsin
#include<cmath>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int a[]={
   2,3,5,7,11,13,17,19,23,29};ll maxs=0;ll mul(ll a,ll b,ll p){ll d=(long double)a/p*b+0.5;ll ret=a*b-d*p;if(ret<0) ret+=p;return ret;
}ll Pow(ll x,ll p,ll MOD)
{ll ret=1;for(;p;p>>=1,x=mul(x,x,MOD)) if(p&1) ret=mul(ret,x,MOD);return ret; 
}inline int MR(ll n)
{if(n==2) return 1;ll t=n-1;int cnt=0;while((t&1)==0) t>>=1,++cnt;for(int i=0;i<10;++i){if(a[i]==n) return 1;ll fac=Pow(a[i],t,n),nxt;for(int i=1;i<=cnt;++i){nxt=mul(fac,fac,n);if(nxt==1&&fac!=1&&fac!=n-1) return 0;fac=nxt;}if(fac!=1) return 0;}return 1;
}ll gcd(ll a,ll b){
   return b?gcd(b,a%b):a;}ll pol_rho(ll n,ll c)
{ll k=2,x=rand()%n,y=x,p=1;for(ll i=1;p==1;++i){y=(mul(y,y,n)+c)%n;p=x>y?x-y:y-x;p=gcd(p,n);if(i==k) x=y,k<<=1;}return p;
}inline void solve(ll n)
{
// printf("%lld\n",n);if(MR(n)){maxs=max(maxs,n);return ;}ll t=n;while(t==n||t==1) t=pol_rho(n,rand()%(n-1));solve(n/t);solve(t);
}int main()
{srand(20010827);int T;cin>>T;while(T--) {ll n;cin>>n;maxs=0;solve(n);if(maxs==n) puts("Prime");else printf("%lld\n",maxs);}return 0;
}