Miller-Rabin算法与Pollard's Rho算法总结

Miller-Rabin素数测试

---

费马小定理：

当$p$是素数时，对于$a\in [1,p-1]$，满足

$$a^{p-1}\equiv 1 (mod \ p)$$
但是其逆命题是假命题，也就是说对于 $a\in [1,p-1]$，我们测试所有的 $a^{p-1}$，就算全部 $=1$， $p$也不一定是素数（可能是伪素数），但出错的概率很小。于是可以选择若干个$a$来测试，再加上下文中的二次探测，正确率就十分接近100%了。

二次探测：

首先有一个显然的性质：$x^2\equiv 1(mod \ p)\Leftrightarrow x\equiv 1/x \equiv p-1(mod \ p)$。
我们设$p-1=2^s\cdot r$，序列$b_i=a^{2^i\cdot r}$；
那么$b_s=1$，而且$b_{0..s-1}$要么全部$=1$，要么存在$b_k=p-1$且$b_{k+1..s-1}=1$。
满足上述条件才可能是素数。

总结一下，Miller-Rabin测试就是选择若干个用来测试的$a$（一般选择前$10$个素数即可），进行二次探测，如果全部通过就可以认定其为素数了。

Pollard’s Rho大数质因数分解

---

假设要质因数分解$n$，我们使用一个伪随机数$f(seed)=(seed^2+c) \mod n$。先设两个变量$x=y=seed$，让$y=f(y)$迭代并检测$\gcd(|x-y|,n)$是否$>1$，迭代$k$次之后把$x$赋为$y$并把$k$乘$2$（倍长迭代），直到$\gcd>1$找到了一个因子，或者$x=y$（进入循环）则重新换一个$c$寻找。
该算法期望在迭代$O(\sqrt p)$次找到$n$的一个因子$p$，考虑到$n$的质因子中$>\sqrt n$的最多只有一个，所以复杂度大概是$O(n^{\frac{1}{4}})$的。

板子题：BZOJ3667

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
#define F(x) (mul(x,x,n)+c)%n
using namespace std;
const ll pri[]={2,3,7,61,10007,24251};
ll mx;
ll rd(ll x)
{
    return (((((((ll)rand()<<15)+rand())<<15)+rand())<<15)+rand())%x;
}
ll gcd(ll x,ll y)
{
    if(!y) return x;
    return gcd(y,x%y);
}
ll mul(ll a,ll b,ll mod)
{
    a=(a>=mod?a%mod:a);
    b=(b>=mod?b%mod:b);
    ll tmp=(a*b-(ll)((long double)a/mod*b+1e-2)*mod);
    return tmp<0?tmp+mod:tmp;
}
ll ksm(ll a,ll b,ll mod)
{
    ll re=1;
    for(a%=mod;b;b>>=1,a=mul(a,a,mod))
        if(b&1) re=mul(re,a,mod);
    return re;  
}
ll check(ll a,ll n,ll r,ll s)
{
    ll t=ksm(a,r,n),u=t;
    while(s--)
    {
        t=mul(t,t,n);
        if(t==1&&u!=1&&u!=n-1) return 0;
        u=t;
    }
    if(t!=1) return 0;
    return 1;
}
ll MR(ll n)
{
    ll r=n-1,s=0;
    while(!(r&1)) r>>=1,s++;
    for(int i=0;i<6;i++)
    {
        if(pri[i]==n) return 1;
        if(!check(pri[i],n,r,s)) return 0;
    }
    return 1;
}
ll rho(ll n)
{
    ll a,b,c,d=1,k=2;
    a=b=rd(n);c=rd(n);
    for(ll i=1;d==1;i++)
    {
        a=F(a);
        d=gcd(a>=b?a-b:b-a,n);
        if(i==k) b=a,k<<=1;
    }
    if(d==n) return rho(n);
    return d;
}
void solve(ll x)
{
    if(x<=1) return ;
    if(MR(x)) {mx=max(mx,x);return ;}
    ll t=rho(x);
    solve(t);solve(x/t);
}
int main()
{
    srand(666623333);
    int ca;
    scanf("%d",&ca);
    while(ca--)
    {
        ll x;
        scanf("%lld",&x);
        mx=0;
        solve(x);
        if(mx==x)puts("Prime");
        else printf("%lld\n",mx);
    }
    return 0;
}