[数论模板] BZOJ 3667 Rabin-Miller算法

还是发一下吧 以后找版子容易些

多年前的代码丑

#include<cstdio>
#include<cstdlib>
#include<ctime>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;

inline char nc()
{
	static char buf[100000],*p1=buf,*p2=buf;
	if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }
	return *p1++;
}

inline void read(ll &x)
{
	char c=nc(),b=1;
	for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
	for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

inline void write(ll x,char c)
{
	if (x==0) { putchar('0'); putchar(c); return; }
	if (x<0) putchar('-'),x=-x;
	char s[55]={0},len=0;
	while (x) s[++len]=x%10+'0',x/=10;
	for (int i=len;i;i--) putchar(s[i]); putchar(c);
}

inline ll random(ll n)
{
	return (ll)((double)rand()/(RAND_MAX+1)*n);
}

inline ll gcd(ll a,ll b)
{
	return b!=0?gcd(b,a%b):a;
}

inline ll Mul(ll a,ll b,ll p)
{
	a%=p; b%=p;
	ll t=(a*b-(ll)((double)a/p*b+0.5)*p);
	return t<0?t+p:t;
}

inline ll Pow(ll a,ll b,ll p)
{
	ll t=1;
	for(;b;b>>=1,a=Mul(a,a,p))
		if(b&1)
			t=Mul(t,a,p);
	return t;
} 

inline bool Witness(ll a,ll p,ll s,ll t)
{
	ll x0,x1;
	x0=Pow(a,t,p);
	while (s--)
	{
		x1=Mul(x0,x0,p);
		if (x1==1 && x0!=1 && x0!=p-1)
			return true;
		x0=x1;
	}
	if (x0!=1) return true;
	return false;
}

inline bool Miller(ll p)
{
	if (p<2) return false;
	if (p==2 || p==3) return true;
	if (~p&1) return false;
	if(p%6!=1 && p%6!=5) return false; 
	ll s=0,t=p-1;
	while (~t&1) s++,t>>=1;
	ll a;
	for (int i=1;i<=5;i++)
	{
		a=random(p-1)+1;
		if (Witness(a,p,s,t))
			return false;
	}
	return true;
}

ll n,maximum=0LL;

ll rho(ll n,ll c)  
{  
    ll k=2,x=random(n),y=x,p=1;  
    for(ll i=1;p==1;i++)  
    {  
        x=(Mul(x,x,n)+c)%n;  
        p=abs(x-y);  
        p=gcd(n,p);  
        if(i==k) y=x,k<<=1;  
    }  
    return p;  
}  

void Calc(ll n)  
{  
    if(n==1) return;  
    if(Miller(n))  
    {  
        maximum=max(maximum,n);  
        return;  
    }  
    ll t=n;  
    while(t==n) t=rho(n,random(n-1)+1);  
    Calc(t); Calc(n/t);  
}

int main()
{
	ll Q;
	freopen("t.in","r",stdin);
	freopen("t.out","w",stdout);
	read(Q);
	while (Q--)
	{
		maximum=0LL;
		read(n);
		if (Miller(n))
			printf("Prime\n");
		else
		{
			Calc(n);
			write(maximum,'\n');
		}
	}
	return 0;
}

UPD 新模板get！orz hzq

#include<cstdio>
#include<cstdlib>
#include<ctime>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;

int prime[9]={2,3,5,7,11,13,17,19,23};
unsigned long long RR;
inline long long R(long long Mo){return (RR+=4179340454199820289ll)%Mo;}
inline long long mul(long long x,long long y,long long Mo){
  long long tmp=(x*y-(long long)((long double)x/Mo*y+0.5)*Mo);
  return tmp<0?tmp+Mo:tmp;
}
inline long long Pow(long long x,long long k,long long Mo){
  long long ans=1;
  for(;k;k>>=1,x=mul(x,x,Mo))if(k&1)ans=mul(ans,x,Mo);
  return ans;
}
inline bool MR(long long n){
  if(n<=1)return false;
  for(int i=0;i<9;i++)if(n==prime[i])return true;
  long long d=n-1;int tmp=0;
  while((d&1)==0)d>>=1,tmp++;
  for(int i=0;i<9;i++){
    long long x=Pow(prime[i],d,n),p=x;
    for(int j=1;j<=tmp;j++){
      x=mul(x,x,n);
      if((x==1)&&(p!=1)&&(p!=n-1))return false;
      p=x;
    }
    if(x!=1)return false;
  }
  return true;
}
inline long long f(long long x,long long c,long long Mo){return (mul(x,x,Mo)+c)%Mo;}
inline long long gcd(long long x,long long y){return x==0?y:gcd(y%x,x);}
inline long long check(long long c,long long n){
  long long x=R(n),y=f(x,c,n),p=n;
  while((x!=y)&&((p==n)||(p==1))){
    if(x>y)p=gcd(n,x-y);
    else p=gcd(n,y-x);
    x=f(x,c,n);y=f(f(y,c,n),c,n);
  }
  return p;
}

ll maxv;

inline void rho(long long n){
  if(n<=1)return;
  if(MR(n)){maxv=max(maxv,n);return;}
  while(true){
    long long tmp=check(R(n-1)+1,n);
    if(tmp!=n && tmp!=1){rho(tmp),rho(n/tmp);return;}
  }
}

ll n,Q;

#define read(x) scanf("%lld",&(x))
int main(){
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  read(Q);
  while (Q--){
    maxv=0;
    read(n);
    if (MR(n))
      printf("Prime\n");
    else{
      rho(n);
      printf("%lld\n",maxv);
    }
  }
  return 0;
}