米勒来宾素数测试与Pollard rho算法
米勒来宾素数测试
此处只考虑如何判断奇素数
根据费马小定理：
当\(p\)为（奇）素数，\(\gcd(a,p)=1\)时，\(a^{p-1}\equiv 1\pmod p\)
就可以写一个用快速幂判断一个数是不是奇素数
#include<bits/stdc++.h>
using namespace std;
const int plist[]= {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97};
//素数表
const int L=25;
int fastpow(int a,int k,int mod) {
int base=1;
while(k) {
if(k&1) base=base*1ll*a%mod;
a=a*1ll*a%mod;
k>>=1;
}
return base;
}
bool isprime(int p) {
if(p<=100) {//在表内就直接判断
for(int i=0; i<L; ++i) {
if(plist[i]==p)return true;
}
return false;
}
for(int i=0; i<L; ++i) {//否则看看是否符合费马小定理
if(fastpow(plist[i],p-1,p)!=1)return false;
}
return true;
}
int n;
int main() {
scanf("%d",&n);
printf(isprime(n)?"YES":"NO");
return 0;
}
但这个效率低，而且错误的概率很大。
我们可以发现，当\(p\)为奇素数时，\(p-1\)为偶数，就可以写成\(2^r\times m\)的形式（\(m\)为奇数）。
当\(a^m\equiv 1\pmod p\)时，一定通过测试。
否则当\(p\)为奇素数时，必定有一个数\(k(0\leq k<r)\)使\(a^{2^k\times m}\equiv p-1\pmod p\)，这样再进行平方（指数乘二）后才能是\(1\)。
所以我们可以用\(k\)从\(0\)枚举到\(r-1\)，若有一个\(k\)使\(a^{2^k\times m}\equiv p-1\pmod p\)则通过测试，否则\(p\)不为素数。
同时计算\(a^{2^k\times m}\)不用快速幂，直接用\((a^{2^{k-1}\times m})^2\)计算即可
const int plist[]= {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97};
//素数表
const int L=25;
inline ll mul(const ll&a,const ll&b,const ll&mod) { //可能爆ll的计算a*b%mod
ll r=((long double)a*b/mod+0.5); //注：用这个不要用太大的模数
r=a*b-r*mod;
return r<0?r+mod:r;
}//具体为什么可以这样，搜索“骆可强 快速乘”
inline ll fastpow(ll a,ll k,const ll&mod) { //快速幂
ll ans=1;
while(k) {
if(k&1)
ans=mul(ans,a,mod);
a=mul(a,a,mod);
k>>=1;
}
return ans;
}
bool mr(const ll&n,const ll&a) { //用a来hack掉n，是合数就返回true
ll t=n-1;
const ll p=t;
while(!(t&1))t>>=1; //除到t为奇数
ll buf=fastpow(a,t,n);
if(buf==1||buf==p)return 0; //通过测试
while((t<<=1)<p) { //和上面1到r-1是一样的
buf=mul(buf,buf,n);
if (buf==p)return 0; //通过测试
}
return 1;
}
bool ptest(ll n) {
if(n<2ll)
return 0;
if(!(n&1))return n==2; //偶数的情况
for(int i=0; i<L; ++i) {//小范围的判断整除
if(n==plist[i])return true;
if(!(n%plist[i]))return false;
}
return n<=10000ll?1:!mr(n,2)&&!mr(n,61)&&!mr(n,233);//用这三个素数就够了
}
时间复杂度\(O(\log n)\)
Pollard rho质因数分解
小结论1（生日悖论）
随便找\(23\)个人,他们之中有任意两人生日相同的概率已经大于\(50\%\)
证明：
不相同的概率：\(\dfrac{365}{365}\times\dfrac{364}{365}\times\cdots\times\dfrac{365-23+1}{365}\)
相同的概率即是\(1\)减去上式，手工算/程序算可以算出是大于\(50\%\)的。
扩展：
随机选出\([1,n]\)中的整数，期望写到\(\sqrt{n}\)个整数时出现重复。
算法思路
先通过一个方法找到一个因子\(s\)，之后递归\(s\)和\(\dfrac{n}{s}\)，递归到质数时更新答案即可。
现在要考虑怎么找一个因子。
假设合数\(n\)有一个因子\(q(1<q<n)\)，则我们一定可以找到两个数\(x,y\)使\(x\equiv y\pmod q\)且\(x\not\equiv y\pmod n\)
也就是说，我们可以用\(\gcd(|x-y|,n)\)来求出一个对应的因子\(q\)。
问题便转换成求一个\((x,y)\)。
考虑一个数列\(f_0=1,f_i=F(f_{i-1})\mod n\)
其中\(F\)某种单变量映射,结果和只和参数相关，例如\(F(x)=x^2+c\)（\(c\)可以随机取一个）
这个函数还算随机（即算个伪随机数），所以根据上面的小结论，大约到\(\sqrt{n}\)项就可能有环。
每个\(f_i\)的值即可表示为一个类似于\(\rho(rho)\)的东西：（假设\(q(1<q<n)\)为\(n\)的一个因数）
第一张是模\(q\)意义下的，其中编号为\(i\)个节点表示\(f_i\mod q\)的值
第二张是模\(n\)意义下的，其中编号为\(i\)个节点表示\(f_i\mod n\)的值
可以发现\(f_{11}\equiv f_{4}\pmod q,f_{11}\not\equiv f_4\pmod n\)，也就是出现环了，我们就可以通过\(\gcd(|f_{11}-f_{4}|,n)\)求出\(q\)。
所以当且仅当\(f_{11}\equiv f_4\pmod n\)时匹配失败。
这个情况特别小，因为\(q\)的环一般远远比\(n\)的环要小，所以更换另一个\(c\)再匹配的复杂度可以忽略不计。
接下来就是要考虑判断环。
追及法求环
我们可以用两个变量\(x1,x2\)，初始时\(x1=1，x2=2\)，每次\(x1\)往后跳一个，\(x2\)往后跳两个（即\(v_{x2}=2v_{x1}\)）
也就是看\(f_{x1},f_{x2}\)就是看\(f_i\)和\(f_{2i}\)。
那么如果\(q\)有环,前者一定会被后者追上（即多跑整数个环）。
如果\(f_i=f_{2i}=0\)，则匹配失败。
否则如果\(\gcd(n,|f_{2i}-f_{i}|)>1\)则匹配成功。
根据小结论，这个是期望\(\sqrt{q}\)的。
若失败，则重置\(F(x)=x^2+c\)的\(c\)即可。
小优化
由于\(q\leq\sqrt{n}\)，得\(\sqrt{q}\leq n^{0.25}\)，所以我们要计算\(n^{0.25}\)次\(\gcd\)
也就是说，时间复杂度为\(O(n^{0.25}\log n)\)。
我们发现当\(\gcd(a,n)>1,\gcd(b,n)\geq1\)时，\(\gcd(ab\mod n,n)>1\)
我们可以一次存下\(z\)要测试的个值，同时求出它们的积\(\mod n\)，然后求与\(n\)的\(\gcd\)
大于\(1\)时，这\(z\)个值一定有一个使与\(n\)的\(\gcd\)大于\(1\)，找到它即可。
否则重新再找\(z\)个值，不断重复上面的操作。
一般\(z\)取\(128\)较合适，这样我们就优化掉了一个\(\log\)
时间复杂度\(O(n^{0.25})\)
#include <bits/stdc++.h>
using namespace std;//以下大段为米勒来宾素数测试
const int plist[]= {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97};
const int L=25;
typedef long long ll;
inline ll read() {
bool f=1;
ll x=0;
int c=getchar();
while(c<=47||c>=58)f&=c!=45,c=getchar();
while(c>=48&&c<=57) {
x=(x<<3)+(x<<1)+(c&15);
c=getchar();
}
return f?x:-x;
}
inline ll mul(const ll&a,const ll&b,const ll&mod) {
ll r=((long double)a*b/mod+0.5);
r=a*b-r*mod;
return r<0?r+mod:r;
}
inline ll fastpow(ll a,ll k,const ll&mod) {
ll ans=1;
while(k) {
if(k&1)
ans=mul(ans,a,mod);
a=mul(a,a,mod);
k>>=1;
}
return ans;
}
bool mr(const ll&n,const ll&a) {
ll t=n-1;
const ll p=t;
while(!(t&1))t>>=1;
ll buf=fastpow(a,t,n);
if(buf==1||buf==p)return 0;
while((t<<=1)<p) {
buf=mul(buf,buf,n);
if (buf==p)return 0;
}
return 1;
}
bool ptest(ll n) {
if(n<2ll)
return 0;
if(!(n&1))return n==2;
for(int i=0; i<L; ++i) {
if(n==plist[i])return true;
if(!(n%plist[i]))return false;
}
return n<=10000ll?1:!mr(n,2)&&!mr(n,61)&&!mr(n,233);
}//上面都是米勒来宾素数测试
inline ll gcd(ll a,ll b) {//gcd
a=abs(a),b=abs(b);
if(!a)return b;
for(; b; swap(a,b))
a%=b;
return a;
}
const int lim=128;//z
ll sav[150];//存z个值的数组
int tot;
ll prho(const ll&n,const ll&c) {//追及法求环
ll x1=c+1,x2=mul(x1,x1,n)+c,buf=1;//先取第一个fi
tot=0;
if(x1>=n)x1-=n;
if(x2>=n)x2-=n;
while(1) {
buf=mul(x1-x2,buf,n);//乘上要测试的数的值
++tot;
sav[tot]=x1-x2;//存下要测试的数
if(tot==lim) {//够z个数了
if(gcd(buf,n)>1)break;//满足要求
tot=0;//重新开始存
}
x1=mul(x1,x1,n)+c;
x2=mul(x2,x2,n)+c;
x2=mul(x2,x2,n)+c;//往后跳
}
for(int i=1; i<=tot; ++i) {
buf=gcd(sav[i],n);
if(buf>1)return buf;//返回值
}
return n;//不成功
}
void solve(const ll&n,ll&ans) {
if(n<=ans)return ;//剪枝
if(ptest(n)) {//是素数
ans=max(ans,n);
return ;
}
const ll p=n-1;
ll sav=prho(n,rand()%p+1);
while(sav==n)sav=prho(n,rand()%p+1);//没成功就换一个c值
solve(sav,ans);//递归
solve(n/sav,ans);
}
ll pollard_rho(ll n) {
ll ans=0;
if(!(n&1)) {
ans=2;
while(!(n&1))n>>=1;
}
const int p=sqrt(min(10000ll,n));
for(int i=3; i<=p; i+=2)
if(!(n%i)) {
ans=i;
n/=i;
while(!(n%i))n/=i;
}//上面是先把小范围的素数判了
solve(n,ans);
return ans;
}
int main() {
srand(233);
for(int T=read(); T; --T) {
ll n=read();
ll ans=pollard_rho(n);
if(ans==n)printf("Prime\n");
else printf("%lld\n",ans);
}
return 0;
}
