Loading [MathJax]/jax/element/mml/optable/MathOperators.js

浅谈 Miller-Robbin 与 Pollard Rho


undefinedundefined

前言

MillerRobbinPollardRho 虽然都是随机算法,不过用起来是真的爽。

MillerRabin 算法是一种高效的质数判断方法。虽然是一种不确定的质数判断法,但是在选择多种底数的情况下,正确率是可以接受的。

PollardRho 是一个非常玄学的方式,用于在 O(n1/4) 的期望时间复杂度内计算合数n的某个非平凡因子。

事实上算法导论给出的是 O(p)pn 的某个最小因子,满足 pnp 互质。

但是这些都是期望,未必符合实际。但事实上 PollardRho 算法在实际环境中运行的相当不错。

注:以上摘自洛谷。

Miller-Robbin

前置芝士

1. 费马小定理

  • 内容:若 φ(p)=p1,p>1,则apa(modp)ap11(modp),(a<p)

  • 证明:戳这里

2. 二次探测定理

  • 内容:如果 φ(p)=p1,p>1,p>X ,且X21(modp),那么X=1 or p1

  • 证明:

X^2\equiv 1\pmod{p}

\therefore p|X^2-1

\therefore p|(X+1)(X-1)

\because p是大于X的质数

\therefore p=X+1\ or\ p\equiv X-1\pmod{p},即X=1\ or\ p-1

算法原理

由费马小定理,我们可以有一个大胆的想法:满足 a^{p-1}\equiv 1\pmod{p} 的数字 p 是一个质数。

可惜,这样的猜想是错误的,可以举出大量反例,如:2^{340}\equiv 1\pmod{341},然鹅 341=31*11

所以,我们可以取不同的 a 多验证几次,不过,\forall a<561,\,a^{560}\equiv 1\pmod{561},然鹅 561=51*11

这时,二次探测就有很大的用途了。结合费马小定理,正确率就相当高了。

这里推荐几个 a_i 的值: 2,3,5,7,11,61,13,17。用了这几个 a_i,就连那个被称为强伪素数的 46856248255981 都能被除去。

主要步骤

  • .将 p-1 提取出所有 2 的因数,假设有s 个。设剩下的部分为 d(这里提取所有2的因数,是为了下面应用二次探测定理) 。

  • 枚举一个底数 a_i 并计算 x=a_i^{d}\pmod p

  • y=x^{2}\pmod p,如果没有出现过 p-1,那么 p 未通过二次探测,不是质数。

  • 否则,若底数已经足够,则跳出;否则回到第二步。

简易代码

#define ll long long
ll p,a[]={2,3,5,7,11,61,13,17};
inline ll mul(ll a,ll b,ll mod)
{
ll ans=0;
for(ll i=b;i;i>>=1)
{
if(i&1) ans=(ans+a)%p;
a=(a<<1)%p;
}
return ans%p;
}
inline ll Pow(ll a,ll b,ll p)
{
ll ans=1;
for(ll i=b;i;i>>=1)
{
if(i&1) ans=mul(ans,a,p);
a=mul(a,a,p);
}
return ans%p;
}
bool Miller_Robbin(ll p)
{
if(p==2) return true;
if(p==1 || !(p%2)) return false;
ll d=p-1;int s=0;
while(!(d&1)) d>>=1,++s;
for(int i=0;i<=8 && a[i]<p;++i)
{
ll x=Pow(a[i],d,p),y=0;
for(int j=0;j<s;++j)
{
y=mul(x,x,p);
if(y==1 && x!=1 && x!=(p-1)) return false;
x=y;
}
if(y!=1) return false;
}
return true;
}

Pollard Rho

大致流程

  • 先判断当前数是否是素数(这里就可应用 Miller-Robbin ),如果是则直接返回

  • 如果不是素数的话,试图找到当前数的一个因子(可以不是质因子)

  • 递归对该因子和约去这个因子的另一个因子进行分解

如何找因子

一个一个试肯定是不行的。而这个算法的发明者采取了一种清奇的思路。(即采取随机化算法)

  • 我们假设要找的因子为p

  • 随机取一个 x、y,不断调整 x ,具体的办法通常是 x=x*x+c(c是随机的,也可以自己定)。

  • p=gcd(y-x,n) ,若p \in \left(1,n\right) ,则找到了一个因子,就返回。

  • 如果出现 x=y 的循环,就说明出现了循环,并不断在这个环上生成以前生成过一次的数,所以我们必须写点东西来判环:我们可以用倍增的思想,让y记住x的位置,然后x再跑当前跑过次数的一倍的次数。这样不断让y记住x的位置,x再往下跑,因为倍增所以当x跑到y时,已经跑完一个圈

  • 同时最开始设定两个执行次数i=1、k=2(即倍增的时候用) ,每次取 gcd++i ;如果 i==k ,则令 y=x ,并将 k 翻倍。

完整代码

#include<cstdio>
#include<ctime>
#include<map>
#include<algorithm>
using namespace std;
#define rg register int
#define V inline void
#define I inline int
#define db double
#define B inline bool
#define ll long long
#define F1(i,a,b) for(rg i=a;i<=b;++i)
#define F3(i,a,b,c) for(rg i=a;i<=b;i+=c)
#define ed putchar('\n')
#define bl putchar(' ')
template<typename TP>V read(TP &x)
{
TP f=1;x=0;register char c=getchar();
for(;c>'9'||c<'0';c=getchar()) if(c=='-') f=-1;
for(;c>='0'&&c<='9';c=getchar()) x=(x<<3)+(x<<1)+(c^48);
x*=f;
}
template<typename TP>V print(TP x)
{
if(x<0) x=-x,putchar('-');
if(x>9) print(x/10);
putchar(x%10+'0');
}
int T;
ll n,ans;
ll a[]={2,3,5,7,11,61,13,17,24251};
template<typename TP>inline ll Gcd(TP a,TP b) {return !b?a:Gcd(b,a%b);}
template<typename TP>inline ll mul(TP a,TP b,TP p)
{
ll ans=0;
for(TP i=b;i;i>>=1)
{
if(i&1) ans=(ans+a)%p;
a=(a<<1)%p;
}
return ans%p;
}
template<typename TP>inline ll Pow(TP a,TP b,TP p)
{
ll ans=1;
for(TP i=b;i;i>>=1)
{
if(i&1) ans=mul(ans,a,p);
a=mul(a,a,p);
}
return ans%p;
}
B Miller_Robbin(ll n)
{
if(n<2) return false;
if(n==2) return true;
if(n%2==0) return false;
ll d=n-1;int s=0;
while(!(d&1)) d>>=1,++s;
for(rg i=0;i<=8 && a[i]<n;++i)
{
ll x=Pow(a[i],d,n),y=0;
F1(j,0,s-1)
{
y=mul(x,x,n);
if(y==1 && x!=1 && x!=(n-1)) return false;
x=y;
}
if(y!=1) return false;
}
return true;
}
inline ll Pollard_Rho(ll n)
{
ll x,y,c,i,k;
while(true)
{
ll x=rand()%(n-2)+1;
ll y=rand()%(n-2)+1;
ll c=rand()%(n-2)+1;
i=0,k=1;
while(++i)
{
x=(mul(x,x,n)+c)%n;
if(x==y) break;
ll d=Gcd(abs(y-x),n);
if(d>1) return d;
if(i==k) {y=x;k<<=1;}
}
}
}
V Find(ll n)
{
if(n==1) return;
if(Miller_Robbin(n)) {ans=max(ans,n);return;}
ll p=n;
while(n<=p) p=Pollard_Rho(p);
Find(p);Find(n/p);
}
int main()
{
read(T);srand(time(0));
while(T--)
{
ans=0;
read(n);Find(n);
if(ans==n) puts("Prime");
else print(ans),ed;
}
return 0;
}

emmmm \cdots

这数据也太毒瘤了吧!!

看来要疯狂卡常了

优化1(不如叫做卡常?)

蛋定的分析一波,我们发现除了 Pollard-RhoO(n^{1/4}) 的期望时间复杂度外, gcd 和龟速乘都是 O(\log N) 的。

虽然这种复杂度已经很优秀了,可对于本题的数据(T≤3501≤n≤10^{18}),还是太 \cdots

所以我们要果断摒弃这种很 low 的龟速乘,改用一种暴力溢出的快速乘:

  • 简单原理: a \times b \ \ mod \ \ p=a \times b−\left \lfloor \frac{a \times b}{p} \right \rfloor

  • long double 来处理这个 \left \lfloor \frac{a \times b}{p} \right \rfloor

  • 然后处理一下浮点误差就可以了。

  • 模数较大时可能会出锅。

  • 不过出锅概率很小 \cdots

如下

inline ll mul(ll a,ll b,ll mod)
{
a%=mod,b%=mod;
ll c=(long double)a*b/mod;
ll ret=a*b-c*mod;
if(ret<0) ret+=mod;
else if(ret>=mod) ret-=mod;
return ret;
}

实践证明,战果辉煌:6pts -> 94pts !!!

优化2(正解)

虽然关于龟速乘的 O(\log n) 的恶劣影响得到了一定遏制,不过,我还是好想 AC 啊!

通过办法1打表 \cdots

正确 AC 姿势如下:

- 我们发现在 Pollard-Rho 中如果长时间随机化而得不到结果,gcd带来的 O(\log n) 还是很伤肾的!!那有没有办法优化呢?答案是肯定的。
  • 在生成x的操作中,龟速乘所模的数就是n,而要求的就是n的某一个约数,即现在的模数并不是一个质数

  • 根据取模的性质:如果模数和被模的数都含有一个公约数,那么这次模运算的结果必然也会是这个公约数的倍数。所以如果我们将若干个(y−x) 相乘,因为模数是 n ,所以如果若干个(y−x)中有一个与n有公约数,最后的结果定然也会含有这个公约数

  • 所以可以多算几次(y−x)的乘积再来求gcd (一般连续算127次再求一次gcd)。

  • 不过\cdots,如果在不断尝试x的值时碰上一个环,就可能会还没算到127次就跳出这个环了,就无法得出答案;同时,可能127次计算之后,所有(y−x)的乘积都变成了n的倍数(即\prod_{i=1}^{127} {(y-x)} \equiv 0 \pmod{n}

  • 所以我们可以不仅在每计算127次之后求gcd、还要在倍增时(即判环时)求gcd,这样既维护了其正确性,又判了环!!

  • 完整AC代码:

#include<cstdio>
#include<ctime>
#include<algorithm>
using namespace std;
#define rg register int
#define V inline void
#define I inline int
#define db double
#define B inline bool
#define ll long long
#define F1(i,a,b) for(rg i=a;i<=b;++i)
#define F3(i,a,b,c) for(rg i=a;i<=b;i+=c)
#define ed putchar('\n')
#define bl putchar(' ')
template<typename TP>V read(TP &x)
{
TP f=1;x=0;register char c=getchar();
for(;c>'9'||c<'0';c=getchar()) if(c=='-') f=-1;
for(;c>='0'&&c<='9';c=getchar()) x=(x<<3)+(x<<1)+(c^48);
x*=f;
}
template<typename TP>V print(TP x)
{
if(x<0) x=-x,putchar('-');
if(x>9) print(x/10);
putchar(x%10+'0');
}
int T;
ll n,ans;
ll a[]={2,3,5,7,11,61,13,17};
template<typename TP>inline ll Gcd(TP a,TP b) {return !b?a:Gcd(b,a%b);}
inline ll mul(ll a,ll b,ll mod)
{
a%=mod,b%=mod;
ll c=(long double)a*b/mod;
ll ret=a*b-c*mod;
if(ret<0) ret+=mod;
else if(ret>=mod) ret-=mod;
return ret;
}
template<typename TP>inline ll Pow(TP a,TP b,TP p)
{
ll ans=1;
for(TP i=b;i;i>>=1)
{
if(i&1) ans=mul(ans,a,p);
a=mul(a,a,p);
}
return ans%p;
}
B Miller_Robbin(ll n)
{
if(n<2) return false;
if(n==2) return true;
if(n%2==0) return false;
ll d=n-1;int s=0;
while(!(d&1)) d>>=1,++s;
for(rg i=0;i<=8 && a[i]<n;++i)
{
ll x=Pow(a[i],d,n),y=0;
F1(j,0,s-1)
{
y=mul(x,x,n);
if(y==1 && x!=1 && x!=(n-1)) return false;
x=y;
}
if(y!=1) return false;
}
return true;
}
inline ll Pollard_Rho(ll n)
{
while(true)
{
ll x=rand()%(n-2)+1;
ll y=rand()%(n-2)+1;
ll c=rand()%(n-2)+1;
ll i=0,k=1,b=1;
while(++i)
{
x=(mul(x,x,n)+c)%n;
b=mul(b,abs(y-x),n);
if(x==y || !b) break;
if(!(i%127) || i==k)
{
ll d=Gcd(b,n);
if(d>1) return d;
if(i==k) y=x,k<<=1;
}
}
}
}
V Find(ll n)
{
if(n<=ans) return;
if(Miller_Robbin(n)) {ans=max(ans,n);return;}
ll p=Pollard_Rho(n);
while(n%p==0) n/=p;
Find(p),Find(n);
}
int main()
{
read(T);srand(time(0));
while(T--)
{
ans=0;
read(n);Find(n);
if(ans==n) puts("Prime");
else print(ans),ed;
}
return 0;
}
posted @   EinNiemand  阅读(1388)  评论(4编辑  收藏  举报
编辑推荐:
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· .NET Core 托管堆内存泄露/CPU异常的常见思路
· PostgreSQL 和 SQL Server 在统计信息维护中的关键差异
· C++代码改造为UTF-8编码问题的总结
阅读排行:
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· 清华大学推出第四讲使用 DeepSeek + DeepResearch 让科研像聊天一样简单!
· 实操Deepseek接入个人知识库
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库
点击右上角即可分享
微信分享提示