[note] 素数判定与分解质因数

在某些毒瘤的数论题中，可能出现对 \(10^{18}\) 的范围内的数质因数分解的情况。这时，可以使用 Fermat 和 Miller-Rabin 算法进行素性判定，Pollard-Rho 算法寻找非平凡因子，两者结合以快速质因数分解。

Fermat 素性检验

Fermat 素性检验是最基本的概率性素性检验。该方法依靠费马小定理进行：

费马小定理
设 \(p\) 是素数，则对于任意的 \(a\) 满足 \(p\nmid a\)，都有 \(a^{p-1}\equiv 1 \pmod p\)

证明

首先证明对于 \(i=1,2,\dots,a-1\)，都有 \(ia\bmod p\) 互不相等。考虑反证，假设 \(ia\bmod p=ja\bmod p\)，那么 \(ia\equiv ja\pmod{p}\iff(i-j)a\equiv 0\pmod{p}\)，然而 \((i-j)\) 和 \(a\) 均不是 \(p\) 的倍数，所以矛盾。

那么 \(ia\bmod p\) 是一个 \(1\sim p-1\) 的排列，则

\[\prod_{i=1}^{p-1}i=\prod_{i=1}^{p-1}ia\bmod p\equiv\prod_{i=1}^{p-1}ia=a^{p-1}\prod_{i=1}^{p-1}i\pmod{p} \]

这说明

\[(a^{p-1}-1)\prod_{i=1}^{p-1}i\equiv0\pmod{p} \]

由于 \(\prod_{i=1}^{p-1}i\) 不是 \(p\) 的倍数，则 \(a^{p-1}-1\) 是 \(p\) 的倍数，也就是 \(a^{p-1}\equiv1\pmod{p}\)

因此，我们猜想费马小定理的逆命题：如果对于任意的 \(a\) 满足 \(p\nmid a\) 有 \(a^{p-1}\equiv1\pmod p\)，那么 \(p\) 为质数。我们在 \([2,p-1]\) 内随机几个 \(a\)，这样 \(p\nmid a\) 一定成立，如果存在 \(a\) 不满足 \(a^{p-1}\equiv1\pmod{p}\)，\(p\) 就是合数，否则就是质数。

然而费马小定理的逆命题是不成立的。事实上，对于某些合数，就算 \([2,p-1]\) 内的所有数都检验一遍都满足 \(a^{p-1}\equiv1\pmod{p}\)，比如 Carmichael 数，遇到这种数就原地歇菜了。因此，Fermat 素性检验是不够强的，下文将讲解更强的方法。

参考代码

long long qpow(long long a,long long b,long long mod)
{
	long long res=1;
	while(b)
	{
		if(b&1) res=res*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return res;
}
bool Fermat(long long p)
{
	for(int i=1;i<=8;i++)
	{
		int a=rand()%(n-2)+2;
		if(qpow(a,p-1)!=1) return false;
	}
	return true;
}

Miller-Rabin 素性检验

Miller-Rabin 是一种正确率更高的素性检验算法。首先引入一个定理：

二次探测定理
对于质数 \(p\)，\(x^2\equiv1\pmod p\) 的解为 \(x\equiv1\pmod{p}\) 或 \(x\equiv p-1\pmod{p}\)

证明

证明是简单的。\(x^2\equiv1\pmod{p}\iff x^2-1=(x+1)(x-1)\equiv0\pmod{p}\)，则 \(x-1\equiv0\pmod{p}\) 或 \(x+1\equiv0\pmod{p}\)，得证。

根据这个，我们得到定理：

对于奇质数 \(p\)，令 \(d\times2^r=p-1\)（\(a\) 为奇数），则对于任意的 \(a\le p\)，以下两个条件至少有一个成立：

1. \(a^d\equiv1\pmod{p}\)；
2. 存在 \(0\le i<r\)，满足 \(a^{d\times 2^i}\equiv p-1\pmod{p}\)。

证明

令 \(s_i=a^{d\times2^{i}}\)，由费马小定理，\(s_r=a^{d\times2^r}=a^{p-1}\equiv1\pmod{p}\)。设 \(s_i\) 为最小的满足 \(s_i\equiv1\pmod p\) 的 \(i\)，则分类讨论：

1. \(i=0\)：则 \(s_i=a^d\equiv1\pmod{p}\)；
2. \(i>0\)：则 \(s_i=s_{i-1}^2\equiv1\pmod{p}\)，由于 \(i\) 是最小的，则由二次探测定理得 \(s_{i-1}\equiv p-1\pmod{p}\)

猜想这个定理的反命题：随机几个 \(a\)，判断一下符不符合这个定理，不符合就是合数，否则是质数。

遗憾的是，随机 \(a\) 能正确通过 Miller-Rabin 的概率至多为 \(\frac{1}{4}\)。不过在 OI 中，使用 \(\{2,325,9375,28178,450775,9780504,1795265022\}\) 或 \(\{2,3,5,7,11,13,17,19,23,29,31,37\}\) 中的数作为 \(a\) 进行检测，能保证 long long 范围内的数 100% 正确判断。

实现细节上，注意 long long 和 __int128 的使用。

参考代码

long long pr[7]={2,325,9375,28178,450775,9780504,1795265022};
long long qpow(long long a,long long b,long long p)
{
    int128 res=1;
    while(b)
    {
        if(b&1) res=res*a%p;
        a=(int128)a*a%p;
        b>>=1;
    }
    return res;
}
bool Miller_Rabin(long long x,long long p)
{
    if(x<3||x%2==0) return x==2;
    long long d=x-1,c=0;
    while(d%2==0) d>>=1,c++;
    int128 t=qpow(p,d,x);
    if(t==1) return true;
    for(int i=0;i<c;i++)
    {
        if(t==x-1) return true;
        t=t*t%x;
    }
    return false;
}
bool isprime(long long x)
{
    if(x<V) return !p[x];
    if(x<=1) return false;
    for(int i=0;i<7;i++)
    {
        if(x==pr[i]) return true;
        if(x%pr[i]==0) return false;
        if(!Miller_Rabin(x,pr[i])) return false;
    }
    return true;
}

Pollard-Rho 算法

Pollard-Rho 算法是一种随机化算法，能在期望 \(O(\sqrt{p})\) 的时间复杂度内找到 \(n\) 的一个非平凡因子（即除了 \(1\) 或 \(n\) 以外的因子），其中 \(p\) 为 \(n\) 的最小质因数。

首先有一个生日悖论：

生日悖论
随机生成一个值域在 \([1,n]\) 的序列，第一个重复出现的数字前面期望有 \(\sqrt\frac{\pi n}{2}\) 个数。

因此，我们随机生成一个序列 \(x\)，期望在生成到 \(\sqrt{p}\) 个数的时候就出现 \(x_i\equiv x_j\pmod{p}\)，此时 \(\gcd(|x_i-x_j|,n)\) 即为所求。但由于 \(p\) 是不确定的，我们只能用 \(\gcd(|x_i-x_j|,n)>1\) 这个条件判断每一对 \((x_i,x_j)\)，复杂度退化到 \(O(\sqrt{n})\) 的了。我们应重新设计 \(x\) 的生成方式，使得不用对每一组 \((x_i,x_j)\) 进行判断。

设计函数 \(f(x)=(x^2+c)\bmod n\)（\(c\) 为随机的常数），令伪随机序列 \(x\) 满足 \(x_0=0\)，\(x_i=f(x_{i-1})\)，将推导过程画图，会长这样：

注意到 \(f\) 函数的取值为 \([0,n-1]\) 之间，因此必然会形成如图所示的 \(\rho\) 形。

为什么要这样设计呢？因为 \(f\) 函数有神秘性质：若 \(i\equiv j\pmod{p}\)，则 \(f(i)\equiv f(j)\pmod{p}\)，对 \(i\) 和 \(j\) 做 \(f\) 映射相当于两个点都在图中走了一步，两个点距离不变。因此，对于 \(i-j\) 相等的所有 \((x_i,x_j)\) 只需检验一遍就够。

实现上，维护两个指针，一快一慢，快的一次跳 \(2\) 步，慢的一次跳 \(1\) 步，这样就能检查所有距离的 \((x_i,x_j)\)，进入环后快的追慢的还会相遇。复杂度 \(O(n^{\frac{1}{4}}\log n)\)。

注意 \(\gcd(|x_i-x_j|,n)\) 可能为 \(n\)，或者啥都没得到，此时应该重选 \(c\) 再做一遍。一定记得先对 \(n\) 做素性检验，防止质数把算法搞炸。

参考代码

//代码来源：知乎 Pecco，侵删
long long Pollard_Rho(long long N)
{
    if (N == 4)
        return 2;
    if (is_prime(N))
        return N;
    while (1)
    {
        long long c = rand()%(n-1)+1; 
        auto f = [=](long long x) { return ((__int128)x * x + c) % N; };
        long long t = f(0), r = f(f(0));
        while (t != r)
        {
            long long d = gcd(abs(t - r), N);
            if (d > 1)
                return d;
            t = f(t), r = f(f(r));
        }
    }
}

参考文章

https://www.cnblogs.com/biyimouse/p/18924291

https://zhuanlan.zhihu.com/p/267884783

https://oi-wiki.org/math/number-theory/pollard-rho/

https://oi-wiki.org/math/number-theory/prime/