Miller-Rabin 素性测试 与 Pollard Rho 大整数分解

(\)

Miller-Rabin 素性测试

---

考虑如何检验一个数字是否为素数。

经典的试除法复杂度 (O(sqrt N)) 适用于询问 (Nle 10^{16}) 的时候。

如果我们要把询问范围加到 (10^{18}) ，再多组询问呢？

Miller 和 Rabin 建立了Miller-Rabin 质数测试算法。

(\)

Fermat 测试

首先我们知道费马小定理：

[a^{p-1}equiv 1pmod p ]

当且仅当 (p) 为素数时成立。

逆命题是否成立呢？也就是说，如果对于任何的 (ale p) ，都有满足费马小定理，(p) 能否确定为素数？

答案是否定的，因为人们发现了 Carmichael数 ，最小的仅 (561) 。

(\)

二次探测定理

Miller-Rabin 素性测试其实是对上一方法的加强。

首先我们特判掉 (2) ，这样素数就只剩下奇素数。

如果 (p) 是奇素数，那么方程 (a^2equiv1pmod p) 的解只有两个：(a=1) 或 (a=p-1) 。

证明：若非模意义下，方程的解为 (pm 1) ，而模意义下数域为 ([0,p-1]) 所以 (-1) 在模意义下下为 (p-1) 。

(\)

Miller-Rabin 素性测试

对于一个待检验的数 (p) ，我们首先特判掉 (ple 2) 及 (p) 为偶数的情况。

然后选取随机数 (x) 去检验 (p) 。

首先若 (x,p) 不满足费马小定理，显然 (p) 不为素数。

若满足，我们继续讨论：

若当前质数是偶数，则我们使用二次探测定理继续检验：

根据二次探测定理，我们可以计算 (x^{frac{p-1}2}) 在模 (p) 意义下的值。

若得到的值不为 (1) 或 (p-1) ，则 (p) 不为素数。

若得到的值为 (p-1) ，或当前指数为奇数，则无法继续探测，递归结束。

若的到的值为 (1) ，且指数为偶数，则递归下去。

代码写出来大概是这样，注意 (10^{18}) 级别模数下的快速幂乘法会炸，需要使用手写快速乘。

const ll prm[12]={2,3,5,7,11,13,17,19,23,29,31,37};

inline bool MR(ll x){
  if(x==2) return 1;
  if(x<=1||x%2==0) return 0;
  for(R ll i=0,nowt,res;i<12;++i){
    if(prm[i]==x) return 1;
    nowt=x-1;
    if(qpow(prm[i],nowt,x)!=1ll) return 0;
    while(!(nowt&1)){
      nowt>>=1;
      res=qpow(prm[i],nowt,x);
      if(res!=1ll&&res!=x-1) return 0;
      if(res==x-1) break;
    }
  }
  return 1;
}

(\)

正确性及复杂度证明

已经证明，使用前 (12) 个素数去测试（就是上面给出来的），在 (10^{18}) 范围内所有数字的检测结果都是正确的。

博主太菜就不证明啦（但是结论是对的）

复杂度：直接对着代码算，单次检验复杂度（对一个 (10^{18}) 级别的数检验）： (12 imes log(p) imes log(p)^2approx2 imes 10^6)

其中后两个 (log) 是快速幂和快速乘。

要是快速乘太慢可以考虑网上的黑科技（我不确定正确性）

inline ll mul(ll x,ll y,ll p){
  ll tmp=(x*y-(ll)((long double)x/p*y+1e-8)*p);
  return tmp<0?tmp+p:tmp;
}

(\)

Pollard Rho 大整数分解

---

考虑如何求出一个数的标准分解。

经典的试除法复杂度还是 (O(sqrt N)) ，处理不了 (10^{18}) 级别的询问。

于是有了 Pollard Rho 算法。

(\)

生日悖论

此处的悖论意味反常识。

如果随机选 (30) 个人，他们之中至少有两人的生日是同一天的概率有多大？

我们可能会觉得很小，事实上 (23) 人时这个概率已经超过了 (50\%) 。

当选择 (100) 人时，这个概率已经达到了 (99.9999692751072\%) 。

(\)

分解思想

简单，暴力的思路：随机化！

我们随机出一个小于 (p) 的数 (x) ，若 (x|p) 则证明我们找到了一个 (p) 的约数，然后递归分解 (x) 和 (p/x) 。

注意找到 (1) 和 (p) 时不算做找到因子 ！

当我们发现当前待分解的数是质数的时候，就停止递归。

但是检验复杂度太大怎么办？Miller-Rabin！问题解决。

还有一个问题，随机化的数字显然不够科学。

此时提出了一个思路：枚举两个数，然后取他们差的绝对值。

这样枚举枚举到约数的概率还是不大。

又提出了基于这个的思路：求 (gcd(p,|a-b|)) 。

这个概率非常大了，因为 (p) 的约数中，每一个数在 (p) 的范围内倍数极多。

注意上一做法概率小是因为 差值不是 (p) 的约数 ，而不是因为 (gcd) 不为 (1) 。

这个做法大概需要枚举到 (sqrt[4]{n}) 级别的数就可以找到答案。

(\)

随机数

其实随机数怎么选取的问题没有解决。

首先是随机数生成器，rand() 函数其实并不可靠。

有一个函数表现特别好，是 (f(x)=(x^2+a)\%p) ，其中 (a) 是随意规定的一个量。

开始随意选取两个数 (x,y) ，之后每一次都让 (x=f(x),y=f(y)) 就可以用了。

还有一个问题是判环，我们发现一些情况下迭代几次就会出现以前的取值。

这时可以强上 (Hash) 表。出现环以后我们令 (a=a+1) 然后继续重新做就可以，复杂度不会受影响。

(\)

Floyd 判环

Floyd 在判环上提出了新的做法。

两个人在环上赛跑， (A) 是 (B) 速度的二倍，两人同起点同时出发，下一次 (AB) 相遇时他们两人一定都跑过一圈了。

利用的就是这个原理。我们让 (x,y) 初值相同，每次令 (x) 只做一次变化，令 (y) 做两次变化，两数再次相同时停止。

实测表现非常优秀。注意两数相同的时候求出 (d) 会等于 (n) ，退出循环。

ll rho(ll n,ll a){
  ll x=2,y=2,d=1;
  while(d==1){
    x=mul(x,x,n)+a;
    y=mul(y,y,n)+a; y=mul(y,y,n)+a;
    d=gcd(abs(x-y),n);
  }
  if(d==n) return rho(n,a+1);
  return d;
}

(\)

Brent 判环

令一种判环的方式。

Brent 每次只计算 (x_k) ,当 (k) 是 (2) 的幂次时，(y=x_k)，每次计算 (d=gcd(x_k−y,n)) 。

理论复杂度比上一方法优秀，实际测试表现一般。退出循环的原理相同。

ll rho(ll n,ll a){
  ll x=2,y=2,d=1,k=0,i=1;
  while(d==1){
    ++k;
    x=mul(x,x,n)+a;
    d=gcd(abs(x-y),n);
    if(k==i){y=x;i<<=1;}
  }
  if(d==n) return rho(n,a+1);
  return d;
}

(\)

一道例题

---

多组询问，判断给出数是否为素数，若不是则输出最大素因子。

提供一份模板 觉得自己常数比较优秀的可以去洛谷交一交

使用的是 Floyd 判环。BZOJ 9848Ms。

#include<cmath>
#include<cstdio>
#include<cctype>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#define R register
#define gc getchar
using namespace std;
typedef long long ll;
 
inline ll rd(){
  ll x=0; char c=gc();
  while(!isdigit(c)) c=gc();
  while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=gc();}
  return x;
}
 
const ll prm[12]={2,3,5,7,11,13,17,19,23,29,31,37};
 
inline ll mul(ll x,ll y,ll p){
  ll tmp=(x*y-(ll)((long double)x/p*y)*p);
  return tmp<0?tmp+p:tmp;
}
 
inline ll qpow(ll x,ll t,ll p){
  ll res=1;
  while(t){
    if(t&1) res=mul(res,x,p);
    x=mul(x,x,p); t>>=1;
  }
  return res;
}
 
inline bool MR(ll x){
  if(x==2) return 1;
  if(x<=1||x%2==0) return 0;
  for(R ll i=0,nowt,res;i<12;++i){
    if(prm[i]==x) return 1;
    nowt=x-1;
    if(qpow(prm[i],nowt,x)!=1ll) return 0;
    while(!(nowt&1)){
      nowt>>=1;
      res=qpow(prm[i],nowt,x);
      if(res!=1ll&&res!=x-1) return 0;
      if(res==x-1) break;
    }
  }
  return 1;
}
 
inline ll Abs(ll x){return x>0ll?x:-x;}
 
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
 
ll rho(ll n,ll a){
  ll x=2,y=2,d=1;
  while(d==1){
    x=mul(x,x,n)+a;
    y=mul(y,y,n)+a;
    y=mul(y,y,n)+a;
    d=gcd(Abs(x-y),n);
  }
  if(d==n) return rho(n,a+1);
  return d;
}
 
ll calc(ll n){
  if(n<=1) return 1;
  if(MR(n)) return n;
  ll t=rho(n,2);
  return max(calc(t),calc(n/t));
}
 
void print(ll x){
  if(!x) return;
  print(x/10);
  putchar('0'+x%10);
}
 
int main(){
  ll t=rd(),n,ans;
  while(t--){
    ans=calc(n=rd());
    if(ans==n) puts("Prime");
    else{print(ans);putchar('
');}
  }
  return 0;
}