Miller_Rabbin&&Pollard_Rho 学习笔记

时间:2023-03-09 05:49:11
Miller_Rabbin&&Pollard_Rho 学习笔记

Miller_Rabbin&&Pollard_Rho 学习笔记

Intro

首先我们考虑这样一个问题

给定一个正整数 \(p \ (p<=1e8)\),请判断它是不是质数

妈妈我会试除法!

于是,我们枚举 \(\sqrt p\) 以内的所有数,就可以非常轻松地得到正确答案。

贴一小段代码

bool check(int n)
{
if(n==1) return false;
for(int i=2;i*i<=n;++i)
if(!(n%i)) return false;
return true;
}

\(Extra\):给定 \(p\ (p<=1e5)\) 个正整数 \(p\ (p<=1e8)\) ,请判断它是不是质数

妈妈我会埃氏筛!

所谓筛法,就是利用了一条最最基本的原理:素数的倍数一定是合数,于是,我们便可以在一定数据范围之内来解决多次询问一个数是否为素数的问题

埃氏筛便是基于这个最简单的思想

时间复杂度约为\(O(nlog_{2}n)\)

再贴一小段代码

void prime()
{
memset(isprime,0,sizeof isprime);
for(int i=2;i<=n;i++)
if(isprime[i]==0)
for(int j=i+i;j<=n;j+=i)
isprime[j]=1;
}

\(Extra^{2}\):给定 \(p\ (p<=1e7)\) 个正整数 \(p \ (p<=1e8)\),请判断它们是不是质数

妈妈我有信仰,我能写出常数极为优秀的埃氏筛

这显然是不现实的,那可以怎样解决这个问题呢?

爸爸我会欧拉筛!

我们会发现,在用埃氏筛解决问题的时候,我们仍然进行了一些不必要的操作,比如\(30\),在程序运行中我们用了\(2,3,5\)三个素数来证明其为合数,这显然是十分不合常理的。

所以我们考虑对其进行优化。

我们考虑,对于每一对\((i,prime[j])\), 显然$i*prime[j] $为合数,我们将其筛去

当 $ i \ mod \ prime[j] =0 $ 时就说明剩余的情况,我们会在以后考虑到,因为每一个合数都可以分解为一个合数(或一个质数)和一个质数的乘积,这样就可以避免重复考虑

时间复杂度约为 \(O(n)\)

再贴一小段代码

void prime()
{
memset(isprime,0,sizeof isprime);
for(int i=2;i<=n;i++)
{
if(isprime[i]==0)
++cnt,prime[t]=i;
for(int j=1;j<=t&&i*prime[j]<=n;++j)
{
isprime[i*prime[j]]=1;
if(!(i%prime[j])) break;
}
}
}

\(Extra^{3}\):给定 \(p\ (p<=1e4)\) 个正整数 \(p\ (p<=1e16)\),请判断它们是不是质数

好像没辙了

Miller-Rabbin 算法

所谓 \(\texttt{Miller-Rabbin}\) 素性测试算法,就是在极短的的时间内,通过一种玄学科学的判断方法,使得判断出错的概率尽可能小。也就是说,这是一个正确率可控,但不稳定的算法

但为什么要学它呢?

至少,在OI界内够用嘛

因为,我们可以通过人为的操作,将其错误率降到我们可以接受的范围内。


前置芝士:欧拉定理

即当 \(p\) 是质数时

\[a^{\varphi (p)}\equiv 1 (mod \ p)
\]

由 \(\varphi (p)=p-1\) 得

\[a^{p-1}\equiv 1 (mod \ p)
\]

于是,既然当 \(p\) 是质数时,\(a^{p-1}\equiv 1 (mod \ p)\) 成立

那么,当 \(a^{p-1}\equiv 1 (mod \ p)\) 成立时,\(p\) 是否为质数呢?

显然不是!!!

那么,我们是否能够可以尝试一种方法,将一个数进行多次测试,这就是这种测试方法的一个初级想法。

但是,是否存在一些合数,能够通过所有的测试呢?

答案是肯定的。这样的书被人们称为\(\texttt{Carmichael}\)数,如 \(561,1105,1729\) ,有兴趣可以点击此链接

所以说,我们有没有一种办法,能够使得正确率进一步提高呢?

答案依然是肯定的。


二次探测

如果\(p\)为质数,且

\[x^2\equiv 1\pmod{p}
\]

那么 $ x=1 \ or \ p-1 $。

这个应该是显然吧,移项直接就可以得到。

然后,我们就可以利用这一个性质来完成 \(\texttt{Miller-Rabbin}\) 素性测试。


具体做法:

  • 将 \(p-1\) 提取出所有 \(2\) 的因数,假设有 \(num\) 个。设剩下的部分为 \(tot\)。

  • 枚举(或随机)一个底数 \(a\) 并计算 \(r=a^{tot}\pmod p\)。

  • 将 \(r\) 连续平方 \(num\) 次。如果没有出现过 \(p-1\) ,那么 \(p\) 未通过二次探测,不是质数。

  • 否则,若已经测试完毕,则跳出。否则回到第二步。

个人常用的模数:\(3,5,7,11,13,37,79,97\) (还是挺好记的吧)

进行一次二次探测的时间复杂度为为 \(O(\log_2n)\)。

所以测试一个数是否为素数的时间最坏为 \(O(T\log_2n)\),其中 \(T\) 为测试次数。

贴代码

typedef long long ll;
bool miller_rabbin(ll x,ll a,ll d)
{
if(x==1||(!(x&1))) return 0;
if(x==2) return 1;
while(!(d&1))
d>>=1;
ll t=ksm(a,d,x);
while(d!=x-1&&t!=1&&t!=x-1)
{
t=ksc(t,t,x);
d<<=1;
}
return t==x-1||(d&1)==1;
}
bool solve(ll n)
{
for(int i=1;i<=6;++i)
{
if(n==te[i]) return 1;
if(!miller(n,te[i],n-1)) return 0;
}
return 1;
}

Part II


Pollard-Rho


A Brief introduction

\(Pollard-Rho\)是一个基于随机算法的大整数分解算法,可以在较短时间内求出其因数。


考虑这样一个问题

从\(\left [2,n \right ]\)中随机一个数字,求其为n的约数的概率。

这个概率显然是很低的,不符合我们的期望

我们考虑优化这个算法

Birthday Trick (生日悖论)

我们稍微修改一下这个问题:

从[1,1000]中随机选取两个数  x 和 y,x−y=42 (x≠y) 的概率是多少?

能够粗略地算出,此时的概率大约为\(\frac{1}{500}\)。

如果我们不再坚持仅选取1个数并且这个数必须为42,而是能够选取2个数并且他们的差刚好等于42,那么,成功的概率被提高了。

如果我们在[1,1000]中选取 k 个数 x1,……,xk, 

取得的k个数中满足xi−xj=42的概率是多少呢?

大约在 k=30 的时候,成功的概率已经超过一半。 

我们有一个大区间[1,1000],但是仅仅生成约 30 个随机数就让我们成功的概率超过了一半,大约在 k=100 的时候,我们几乎可以确保我们是成功的。这是一个非常重要的观察, 它被称为生日悖论

以上摘自某国外论文


于是,我们就可以通过这样一个方法来大幅度的提高程序效率

为了使复杂度较稳定,我们可以使随机数不那么随机

根据\(f(x)=x^2+c\) 来随机生成,复杂度较优

但是这样生成可能有环出现,于是有两种方法可以来解决这个问题

i.floyd法

我们可以设置两个变量,像追击问题那样,让一个变量以另一个变量的二倍速来扩展,这样,当他们再一次相等时,至少跑过了一圈。

ii.倍增优化法

暂时略

然后关于本代码,由于为了通过Luogu上的毒瘤测试点,加了非常多的玄学优化,在此暂时不做叙述。

留坑,待填

贴代码

ll f(ll n,ll c,ll md)
{
return (ksc(n,n,md)+c)%md;
}
ll ans=0;
ll pollard_rho(ll n)
{
ll w=rand()%(n-1)+1;
ll a,b;
a=b=rand()%(n-1)+1;
ll i=0,j=1;
ll p=1;
ll z=1;
while(++i)
{
a=f(a,w,n);
if(b>a)p=b-a;
else p=a-b;
z=ksc(z,p,n);
if(a==b||(!z)) return n;
if(!(i%667)||i==j)
{
p=gcd(z,n);
if(p>1) return p;
if(i==j)
b=a,j<<=1;
}
}
return n;
}
void solve(ll n)//求一个数的最大质因数
{
if(n==1) return ;
if(n<ans) return ;
if(miller_rabbin(n))
{
ans=max(ans,n);
return;
}
ll t=n;
while(t==n)
t=pollard_rho(t);
solve(t);
solve(n/t);
}