Miller_Rabin及Pollard_Rho算法模板

Miller_Rabin算法学习
Pollar_Rho算法学习
Pollard_Rho的代码实现上不同写法会有很大的时间差异，下面这份代码的常数似乎比较优秀，可以学习学习。
Code：

#include<cstdio>
#include<algorithm>
#define LL long long
using namespace std;
int cnt=0;
LL fac[10005],N;
int base[7]={2,3,7,61,24251,19260817};
inline LL mul(LL a,LL b,LL p){LL r=a*b-(LL)((long double)a/p*b+0.5)*p;return r<0?r+p:r;}
inline LL ksm(LL a,LL b,LL p){
	LL s=1;
	for(;b;b>>=1,a=mul(a,a,p)) if(b&1) s=mul(s,a,p);
	return s;
}
LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}
bool Miller_Rabin(LL n){
	if(n<2) return 0;
	for(int i=0;i<=5;i++) 
		if(base[i]==n) return 1; 
		else if(!(n%base[i])) return 0;
	LL d=n-1; int k=0;
	while(!(d&1)) d>>=1,k++;
	for(int s=0;s<=5;s++){
		LL x=ksm(base[s],d,n),y;
		for(int i=1;i<=k&&x!=1;i++,x=y)
			if((y=mul(x,x,n))==1&&x!=1&&x!=n-1) return 0;
		if(x!=1) return 0;
	}
	return 1;
}
const int M = 127;
LL Rho(LL n){
	LL c=rand()%n+1;
	LL x=rand()%n+1,y=x,q=1,d=1;
	for(int k=1;;k<<=1,y=x,q=1){
		for(int i=1;i<=k&&d==1;i++){
			x=(mul(x,x,n)+c)%n;
			q=mul(q,abs(x-y),n);
			if(!(i&M)) d=gcd(q,n);
		}
		if(d>1||(d=gcd(q,n))>1) break;
	}
	return d==n?Rho(n):d;
}
void Pollar(LL n){
	if(n==1) return;
	if(Miller_Rabin(n)) {fac[++cnt]=n;return;}
	LL p=Rho(n);
	while(n%p==0) n/=p;
	Pollar(p),Pollar(n);
}
int main()
{
	scanf("%lld",&N);
	Pollar(N);
	sort(fac+1,fac+1+cnt);
	for(int i=1;i<=cnt;i++) if(fac[i]!=fac[i-1]) printf("%lld\n",fac[i]);
}
/*
bool Miller_Rabin(LL n){
	if(n<2) return 0;
	for(int i=0;i<=5;i++) if(base[i]==n) return 1; else if(n%base[i]==0) return 0;
	LL d=n-1; int k=0;
	while(!(d&1)) d>>=1,k++;
	for(int s=0;s<=5;s++){
		LL x=Pow(base[s],d,n),y;
		for(int i=1;i<=k&&x!=1;i++,x=y)
			if((y=mul(x,x,n))==1&&x!=n-1) return 0;
		if(x!=1) return 0;
	}
	return 1;
}
LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}
LL Rho(LL n,LL c){
	LL x=rand()%n+1,y=x,q=1,d=1;
	for(int k=1;;k<<=1,q=1,y=x){
		for(int i=1;i<=k;i++){
			x=(mul(x,x,n)+c)%n;
			q=mul(q,abs(x-y),n);
			if(!(i&127)&&(d=gcd(q,n))>1) break;
		}
		if(d>1||(d=gcd(q,n))>1) return d==n?Rho(n,c+1):d;
	}
}
void Pollar(LL n){
	if(n==1) return;
	if(Miller_Rabin(n)) {mx=max(mx,n);return;}
	LL p=Rho(n,3);
	while(n%p==0) n/=p;
	Pollar(p),Pollar(n);
}
*/