Miller-Rabin 与 Pollard-Rho 学习笔记

跟着大佬脚步学习:D

Miller-Rabin

题目

给你一个数字 ((zle10^{18}))，判断它是不是质数。

公式

费马小定理：（ (p) 为质数）

[a^{p-1}equiv1pmod p ]

通俗地讲是 (a^{p-1}mod p=1)，但是逆定理不一定成立（大概率成立），即若 (a^{p-1}equiv1pmod p) 成立， (p) 也不一定是质数，不过如果不成立， (p) 一定不是质数。

二次探测定理：若 (p) 为质数（ (2) 除外）且 (a^2equiv1pmod p)，则 (aequiv 1) 或 ((p-1)pmod p) 。

算法

特判 (z) 是不是偶数，直接再见。

现在 (z) 是奇数了， ((z-1)) 为偶数，所以令 (z-1=2^q+p)，然后多次进行下面算法流程：

• 随机选择一个 (ain[1,z-1])，若 (a^{z-1} otequiv1pmod p) ，则它是合数，再见。（费马小定理）
• 否则枚举所有的 (bin [0,q-1])，如果 (Large a^{2^b}mod z ot=1) 但是 (Large a^{2^{b-1} imes p}mod z ot=1) 或 (Large(z-1)) 那么它是合数，再见。 （二次探测定理）

重复若干遍后如果还是没有再见，那么它就是素数了。

时间复杂度是 (O(klog^2n)) 的， (k) 为重复次数。但是如果数据小，效率会低于 (O(sqrt n)) 的暴力，所以酌情选择。

代码

注意要用龟速乘，怕溢出。（用了 LH 大奆的 (O(1)) '鬼'速乘）

Drand() 是小常数生成随机数（其实 rand() 更快，但是一般我们都是用 rand()*rand()+rand() 之类的以提高范围和随机性，这就会变很慢）。

#include<bits/stdc++.h>
#define rep(i,x,y) for(int i=x,WALL=y;i<=WALL;++i)
#define lon long long
using namespace std;
lon rad=1;

int rd(){
   int shu=0;char ch=getchar();
   while( !isdigit(ch) )ch=getchar();
   while( isdigit(ch) )shu=(shu<<1)+(shu<<3)+(ch^48),ch=getchar();
   return shu;
}

lon Drand(){
	unsigned lon x=rad;
	x^=x<<13,x^=x>>17,x^=x<<5;
	rad=x/4;return rad;
}

lon Dmul(lon p,lon q,lon mo) {
    lon shul=p*(q>>32ll)%mo*(1ll<<25)%mo;
    lon shur=p*(q & ( (1ll<<32)-1 ) )%mo;
    return (shul+shur)%mo;
}

lon Dpow(lon p,lon q,lon mo){
	lon tot=1;
	while(q){
		if(q&1)tot=Dmul(tot,p,mo)%mo;
		p=Dmul(p,p,mo)%mo,q=q>>1;
	}
	return tot;
}

bool rabin(lon z){
	if(z<=1)return 0;
	if(z==2)return 1;
	if(!(z&1))return 0;
	lon q=0,p=z-1;
	while(p&1)p=p/2,q++;
	int chktime=10;
	while(chktime--){
		lon x=Drand()%(z-1)+1;
		if(Dpow(x,z-1,z)^1)return 0;
		lon las=1;
		rep(i,0,q-1){
			lon tmp=Dpow(x,(1ll<<i)*p,z);
			if(tmp==1&&las^1&&las^(z-1))return 0;
			las=tmp;
		}
		if(las^1&&las^(z-1))return 0;
	}
	return 1;
}

int main(){
	int T=rd();
	while(T--){
		lon n=rd();
		puts(rabin(n)?"Prime":"notPrime");
	}
	return 0;
}

Pollard-Rho

贴个代码，待填

Luogu模板题，求出一个数的最大质因数

#include<bits/stdc++.h>
#define rep(i,x,y) for(int i=x,WALL=y;i<=WALL;++i)
#define lon long long
#define abs(x) ((x)<0?-(x):(x))
using namespace std;
lon rad=1,ans;

lon rd(){
   lon shu=0;char ch=getchar();
   while( !isdigit(ch) )ch=getchar();
   while( isdigit(ch) )shu=(shu<<1)+(shu<<3)+(ch^48),ch=getchar();
   return shu;
}

lon Drand(){
	unsigned lon x=rad;
	x^=x<<13,x^=x>>17,x^=x<<5;
	rad=x/4;return rad;
}

lon Dmul(lon p,lon q,lon mo){
	p=p%mo,q=q%mo;
    return (p*q-(lon)((long double)p/mo*q)*mo+mo)%mo;
}

lon Dpow(lon p,lon q,lon mo){
	lon tot=1;
	while(q){
		if(q&1)tot=Dmul(tot,p,mo)%mo;
		p=Dmul(p,p,mo)%mo,q=q>>1;
	}
	return tot;
}

bool mrabin(lon z){
	if(z<=1)return 0;
	if(z==2)return 1;
	if(!(z&1))return 0;
	lon q=0,p=z-1;//z=2^q+p
	while(!(p&1))p=p/2,q++;
	int chktime=10;
	while(chktime--){
		lon x=Drand()%(z-1)+1;
		if(Dpow(x,z-1,z)^1)return 0;
		lon las=1;
		rep(i,0,q-1){
			lon tmp=Dpow(x,(1ll<<i)*p,z);
			if(tmp==1&&las^1&&las^(z-1))return 0;
			las=tmp;
		}
		if(las^1&&las^(z-1))return 0;
	}
	return 1;
}

lon Dgcd(lon p,lon q){
	lon tmp;
	while(q)tmp=p,p=q,q=tmp%q;
	return p;
}

lon prho(lon z){
	int len=1;
	lon now=0,tot=1,lead,c=Drand()%(z-1)+1;
	while(1){
		lead=now;
		rep(i,1,len){
			now=(Dmul(now,now,z)+c)%z;
			tot=Dmul(tot,abs(lead-now),z);
			if( !(i%127) ){
				lon gcd=Dgcd(tot,z);
				if(gcd>1)return gcd;	
			}
		}
		lon gcd=Dgcd(tot,z);
		if(gcd>1)return gcd;
		len=len<<1;
	}
}

void work(lon z){
	if(z<=ans||z<2)return;
	if( mrabin(z) ){ans=max(ans,z);return;}
	lon x=z;
	while(x==z)x=prho(z);
	while( !(z%x) )z=z/x;
	work(x),work(z);
}

int main(){
	int T=rd();
	while(T--){
		lon n=rd();
		if( mrabin(n) )puts("Prime");
		else ans=0,work(n),printf("%lld
",ans); 
	}
	return 0;
}