Miller-Rabin 素性测试和 Pollard-Rho 算法 学习笔记

Miller-Rabin 素性测试

用于判断一个数是否是质数。

首先有费马小定理：若 \(\gcd(a,n)=1\) 且 \(n\) 是质数，有 \(a^{n-1}\equiv 1\pmod{n}\)。

假如费马小定理的逆命题成立，只需要随机 \(a\) 使得 \(\gcd(a,n)=1\)，判断是否有 \(a^{n-1}\equiv 1\pmod{n}\)。

然而逆命题有小概率不成立。有这么一种数，满足对于任意 \(\gcd(a,n)=1\) 都有 \(a^{n-1}\equiv 1\pmod{n}\)，但是它是合数。

这需要引入二次探测定理：若 \(n\) 是奇质数，则 \(x^2\equiv 1\pmod{n}\) 的解为 \(x\equiv\pm1\pmod{n}\)。这容易证明，移项得到 \(n\mid(x+1)(x-1)\)，由于 \(n\) 是质数，只能让 \(n\mid x+1\) 或 \(n\mid x-1\)。

结合二次探测定理和费马小定理的逆命题。将 \(n-1\) 分解为 \(2^tu\)，构造数列 \(a^u,a^{2u},a^{4u},\cdots,a^{2^tu}\)，其中每一个都是前一个的平方，而且最后一个模 \(n\) 余 \(1\)，否则就不满足费马小定理了。如果这一串数中出现了不是 \(\pm1\) 的数平方得到 \(1\)，就可以判定这个数为合数。这说明，这串数要么全为 \(1\)，要么在出现了 \(-1\) 后全为 \(1\)。

细节：当 \(a\) 为 \(n\) 的倍数，跳过这一轮，因为此时费马小定理不可用。当 \(a\) 与 \(n\) 不互质，费马小定理也不可用，但此时 \(n\) 一定为合数，而此时 \(a^{n-1}\bmod n\) 必不为 \(1\)，所以不用担心。

复杂度是优秀的 \(O(k\log n)\)，\(k\) 为测试轮数。

在 OI 中，如果是 int 范围内，\(a\) 选 \(2,7,61\) 可以保证正确。如果是 unsigned long long 范围内，\(a\) 可以选 \(2,325,9375,28178,450775,9780504,1795265022\)。选前十二个质数也可以判断 \(2^{64}\) 以内的质数。这样 \(k\) 就是常数了。

long long qmul(long long a,long long b,long long n){
  return ((unsigned long long)a*b-(unsigned long long)((long double)a/n*b)*n+n)%n;
}
long long qpow(long long a,long long b,long long n,long long ans=1){
  for(a%=n;b;b>>=1)b&1&&(ans=qmul(ans,a,n)),a=qmul(a,a,n);
  return ans;
}
const long long A[7]={2,325,9375,28178,450775,9780504,1795265022};
bool isprime(long long n){
  if(n<=2||~n&1)return n==2;
  int t=__builtin_ctzll(n-1);
  long long u=(n-1)>>t,a,v;
  for(int i=0,s;i<7;i++){
    a=A[i],v=qpow(a,u,n);
    if(v<=1)continue;
    for(s=0;s<t;s++,v=qmul(v,v,n))if(v==n-1)break;
    if(s==t)return 0;
  }
  return 1;
}

Pollard-Rho 算法

概述

用于求 \(n\) 的任意一个非平凡因数（即除了 \(1,n\) 以外的任意一个因数）。

首先当然需要 Miller-Rabin 判断这个数是否为质数。

Pollard-Rho 是一个随机化算法，运用到了生日悖论：在 \([0,n-1]\) 中选一些数，只用选比较少的数就会出现两个数相同。

设选了 \(k\) 个数，则选的数互不相同的概率为 \(\prod_{i=0}^{k-1}\frac{n-i}n\)。发现这个函数衰减很快，解得当 \(k\) 为 \(\sqrt n\) 级别时，这个值就已经小于 \(\frac 1 2\)。那么出现两个数相同时 \(k\) 的期望也是 \(O(\sqrt n)\)。

假如随两个数 \(x_1,x_2(0\leq x_1,x_2<n)\)，且 \(\gcd(|x_1-x_2|,n)>1\)，就可以通过 \(\gcd(|x_1-x_2|,n)\) 得到一个非平凡因数。

构造伪随机函数 \(f(x)=x^2+c\)，其中 \(c\) 是随机常数。在生成的数列中挑选这两个数（模 \(n\) 意义下），可以证明期望 \(O(n^{\frac 1 4})\) 就能找到满足条件的数。

令 \(m\) 为 \(n\) 的最小非平凡因数，则 \(m\leq\sqrt n\)。在模 \(m\) 意义下，期望 \(O(\sqrt m)\) 项也就是 \(O(n^{\frac 1 4})\) 项数列就会进入循环节，因此取值只有约 \(O(n^{\frac 1 4})\) 种。可以在期望 \(O(n^{\frac 1 4})\) 的时间内找到两个数在模 \(m\) 意义下相等而在模 \(n\) 意义下不相等，此时 \(m\mid\gcd(|x_1-x_2|,n)\)，满足要求。

构造出来的函数需要有循环节，而且环长短，进入循环快。而 \(f(x)=x^2+c\) 就比较好。

Floyd 判环

搞两个指针指向选择的两个数。其中，一个指针每次迭代一步，另一个每次迭代两步。那么当两个指针指向的数相同，说明遇到了环。形成环时就要重新选择 \(c\) 再次分解。

Brent 判环

考虑倍增环的长度。每次，其中一个指针 \(t\) 向前跳 \(2^i\) 步。如果相遇则找到环，否则让另一个指针 \(s\) 移动到 \(t\) 的位置，增大 \(i\) 继续尝试。

这样做有一个好处：如果 \(t\) 经过的位置中，存在 \(\gcd(|x_s-x_t|,n)>1\)，那么把这些 \(|x_s-x_t|\) 乘起来，得到的 \(\gcd\) 还是大于一，那么每轮只用做一次 \(\gcd\)。

然而，当 \(t\) 移动的距离较大时，长时间不做 \(\gcd\) 可能造成浪费。设一个阈值 \(B\)，\(t\) 每移动 \(B\) 格也强制做一次 \(\gcd\)。

一般取 \(B=127\)，复杂度可以近似看作 \(O(n^{\frac 1 4})\)。

long long F(long long x,long long c,long long n){
  return (qmul(x,x,n)+c)%n;
}
long long pollard_rho(long long n){
  long long c=rand()%(n-1)+1,s=0,t=0,v=1;
  for(int g=1;;g<<=1,s=t){
    for(int i=1;i<=g;i++){
      t=F(t,c,n),v=qmul(v,abs(t-s),n);
      if(!v)return n;
      if(i%127==0){
        long long d=gcd(v,n);
        if(d>1)return d;
      }
    }
    long long d=gcd(v,n);
    if(d>1)return d;
  }
}

如果要分解质因数，还需要递归地分解，复杂度多一只 \(\log\)。

[[数学]]