Miller-Rabin素数测试 与 Pollard-Rho算法
参考文献
\(\text{Miller-Rabin}\) 素数测试
引入
当我们需要判断一个数是否是素数，我们往往会想到一个非常暴力的 \(O(\sqrt n)\) 算法：
- 枚举 \(i=2...\sqrt n\)，逐一判断 \(n\) 是否是 \(i\) 的倍数。如果存在 \(i\) 满足 \(i|n\)，那么 \(n\) 是合数，否则 \(n\) 是素数
bool isprime(ll n)
{
for(ll i=2;i*i<=n;i++)
if(n%i==0) return false;
return true;
}
当 \(n\) 达到 \(10^{18}\) 级别，或者其他情况下时间不允许时，这种方法就会被淘汰。
值得高兴的是，我们还有另一种方法，即—— \(\text{Miller-Rabin}\) 素数测试。
费马小定理
费马小定理是一个著名的定理。
- 费马小定理：对于 \(a\) 与素数 \(p\)，若满足 \(1\le a<p\)，那么有 \(a_{p-1}\equiv 1\pmod p\)。
可惜的是，费马小定理并不能帮我们判断 \(p\) 是否是素数，因为有些合数 \(p\)，也可能满足上式。
考虑其逆反命题，若满足 \(a_{p-1}\not\equiv 1\pmod p\)，那么 \(p\) 是合数，因此我们可以选取若干个 \(a\) 带入式子，判断 \(p\) 是否是合数。
可惜这样做的正确率还是很低，我们需要更加高科技的产物来帮助。
二次探测定理
- 二次探测定理：若 \(p\) 是奇素数，那么同余方程 \(x^2\equiv 1\pmod p\) 的解为 \(x_1\equiv -1\equiv p-1\pmod p,\space x_2\equiv 1\pmod p\)。
这是一个很强大的定理，我们尝试用它来判断 \(p\) 是否为奇素数。
那么如果一正整数 \(c\) 为偶数，并且 \(a^c\equiv 1\pmod p\)，尝试判断 \(a^{\frac c2}\) 模 \(p\) 意义下是否为 \(\pm 1\)。
如果 \(a^{\frac c2}\) 在模 \(p\) 意义下不是 \(1\) 或 \(-1\)，那么 \(p\) 一定不是素数。
否则若是 \(1\)，我们可以令 \(c=\frac c2\)，继续判断；否则，只能结束判断。
算法流程
上面的简易流程给了我们启示。\(c\) 从大到小需要不断计算快速幂，考虑从小到大。
-
先多次随机选取数字 \(a\) ，然后执行以下流程。
-
分解 \(n-1=u\times 2^t\)。
-
如果 \(a^u\equiv 1\pmod n\)，结束本轮测试。
-
令 \(c=u\)。
-
判断在枚举过程中是否出现过 \(a^c\equiv -1\pmod p\) 的情况，若存在则本轮测试成功；否则失败，\(n\) 不是素数。
如果一开始 \(a^u\) 模 \(n\) 是 \(1\)，我们就无法通过二次探测定理来测试。否则，我们在枚举 \(c\) 的过程中，一定是“其中出现过一个 \(a^c\) ，满足模 \(p\) 意义下为 \(-1\)，之后所有数模 \(p\) 意义下都为 \(1\)”。
注意到“随机选取 \(a\)”很吃 \(\text{RP}\)，这里给出通常选取的数。
-
对于 \([1,2^{32})\) 内的数，只需要选取 \(a=2,7,61\) 即可。
-
对于 \([1,2^{64})\) 内的数，只需要选取 \(a=2,325,9375,28178,450775,9780504,1795265022\) 即可。
-
（\(\text{OI}\) 中最常用的选取方法）对于 \([1,2^{78})\) 内的数，只需要选取 \(a=2,3,5,7,11,13,17,19,23,27,29,31,37\) 即可，也就是前 \(12\) 个素数。
理论上时间复杂度为 \(O(k\log ^2n)\)，其中 \(k\) 为选取的素数个数。
const ll Pr[]={0,2,325,9375,28178,450775,9780504,1795265022ll};
ll mul(ll a,ll b,ll p)
{
return (__int128)a*b%p;
}
ll power(ll a,ll b,ll p)
{
ll ans=1;
while(b)
{
if(b&1) ans=mul(ans,a,p);
a=mul(a,a,p);
b>>=1;
}
return ans;
}
bool Miller_Rabin(ll n)
{
if(n<3||n%2==0) return (n==2);
ll u=n-1, t=0;
while(u%2==0) u>>=1, ++t;
for(ll i=1;i<=7;i++)
{
ll tmp=power(Pr[i],u,n);
if(tmp==1) continue;
ll fl=0;
for(ll j=1;j<=t;j++)
{
if(tmp==n-1)
{
fl=1;
break;
}
tmp=mul(tmp,tmp,n);
}
if(!fl) return false;
}
return true;
}
\(\text{Pollard-Rho}\) 算法
引入
- 一个问题：求 \(n\) 的一个非平凡因子，非平凡因子指一个 \(x\) 满足 \(x\not=1,x\not=n,x|n\)。\(n\le 10^{18}\)。
我们完全可以用 \(O(\sqrt n)\) 的时间内解决这个问题。但当 \(n\) 过大时，这种方法显然超时。
生日悖论
- 生日悖论：有 \(k\) 个人，若 \(k\ge 23\)，那么这些人中至少有两人同一天生日的概率 \(\ge 50\%\)。
这个结论似乎与我们的数学直觉很矛盾，考虑证明。设 \(A\) 为这 \(k\) 个人中不存在两人同一天生日的概率，那么
当我们带入 \(k=23\) 时，惊奇的发现 \(P(A)\approx 0.492703,\space 1-P(A)>50\%\)，即至少有两人在同一天生日的概率 \(\ge 50\%\)。
这启示我们，满足答案的组合往往优于单个。
伪随机数列
注意到 \(\gcd(d,n)\) 一定是 \(n\) 的因数，我们可以尝试找一些 \(d\)，满足 \(1<\gcd(d,n)<n\)。
如果 \(n\) 不是素数，设 \(p\) 为 \(n\) 的最小质因子我们可以随机生成 \(\sqrt[4]n\) 个不同数的数列，根据生日悖论，我们期望存在两个数列中的数 \(a,b\)，满足 \(|a-b|\pmod p\)，那么我们就可以令 \(d=|a-b|\)。
\(\text{Pollard}\) 构造了这样一个数列：\(a_1=c\)，对于 \(i(i>1)\)，有 \(a_i=(a_{i-1}^2+c)\mod n\)，其中 \(c\) 是一个随机生成的常量。
- 引理：存在位置 \(i\) 以及一个正整数 \(s\)，满足对于 \(j(j\ge i)\)，有 \(a_{j+s}=a_j\)。即，这个数列在某个位置开始存在循环节。
割巢原理就能证。
如果把这些数放在坐标系上，这些循环节就是一个环，这就非常像希腊字母 \(\rho\)，这也是 \(\text{Pollard-Rho}\) 名字的由来。
根据生日悖论，我们最多走 \(O(\sqrt p)\) 步就能到达环，环长期望 \(O(\sqrt p)\) 也就是 \(O(\sqrt[4]n)\)。
\(\text{Floyd}\) 判环
考虑维护两个位置 \(i,j\)，逐一判断 \(\gcd(|a_i-a_j|,n)\) 是否满足条件。我们可以令 \(i\) 每次 \(+1\)，\(j\) 每次 \(+2\)，形成距离差。只要足够随机，这个方法就是对的。
当 \(a_i=a_j\) 时，即 \(i,j\) 在环上重合于一点，我们直接返回 \(n\)。
倍增优化
我们频繁调用 \(\gcd\) 函数，这可能使程序非常慢。因为 \(a_i<n\)，所以对于所有 \(i,j\) 都有 \(|a_i-a_j|<n\)，\(\gcd\) 有个性质，就是如果 \(\gcd(a,b)>1,\space c>1\)，那么 \(\gcd(ac,b)>1\)。进一步的，\(\gcd(ac,b)=\gcd(ac\mod b,b)\)。
因此，我们可以尝试维护一段区间的的 \(|a_i-a_j|\) 的乘积，最后把乘积与 \(n\) 进行 \(\gcd\) 判断。
设每过 \(2^w-1\) 个数进行一次判断，通常取 \(w=7\) 最优。
ll f(ll x,ll c,ll n)
{
return (mul(x,x,n)+c)%n;
}
ll gcd(ll a,ll b)
{
if(!b) return a;
return gcd(b,a%b);
}
ll Pollard_Rho(ll n)
{
ll s=0, t=0, c=1ll*rand()%(n-1)+1;
ll tmp=1;
for(ll i=1;;i<<=1,s=t,tmp=1)
{
for(ll j=1;j<=i;j++)
{
t=f(t,c,n);
tmp=mul(tmp,abs(t-s),n);
if(j%127==0)
{
ll d=gcd(n,tmp);
if(d>1) return d;
}
}
ll d=gcd(tmp,n);
if(d>1) return d;
}
}
