Miller_Rabin和Pollard_Rho算法

Miller_Rabin和Pollard_Rho算法

SemiWaker

一直觉得随机化不稳定就没学。

Miller_Rabin

目的：判质数。

当然，不是随机一个数然后去模，随机一个数去模正确率太低了。

我们要用数学的方法。

费马小定理：
如果p是质数，，那么有：

$$a^{p-1}\equiv 1\pmod p$$
我们看它的逆命题：
如果对于所有 $1<a<p$，都满足上式，那么p是质数。

然而这个命题是错的，不过还是很靠谱的。

为了使它更加靠谱，要么我们背数表，要么我们加上别的东西去判断。

如果p是质数，那么$x^2\equiv 1\pmod p$只有两个解$x=1$或$x=p-1$
因为$x^2-1\equiv 0\pmod p$，即$(x+1)(x-1)\equiv 0\pmod p$。
显然只有两个解$x=1$或$x=p-1$。

好的，这个命题是对的，但是我们不能暴力验证。

我们将两个东西合起来验证。
首先对于任意一个a，有：$a^{p-1}\equiv 1\pmod p$
我们考虑把p-1分解，$p-1=d\times2^r$
d是奇数。
为什么我们要分解出所有2出来呢？为了构造平方。

$$a^{d\times2^r}\equiv 1\pmod p$$
我们把前面的开方：
$$(a^{d\times2^{r-1}})^2\equiv 1\pmod p$$
如果 $a^{d\times2^{r-1}}$不是1或p-1的话，p就不是质数。
如果 $a^{d\times2^{r-1}}$是1，我们可以继续开方，直到 $a^d$为止。
如果 $a^{d\times2^{r-1}}$是p-1，那么继续开方没有意义了。

实际操作中，从$a^d$开始一直平方下去比较好。
（只用一次快速幂。）
从$a^d$开始不停平方是怎样的呢？
则一列数是这样的：
* $a^d$
* $a^2d$
* $a^4d$
* $a^8d$
* ……
* $a^d\times 2^r$
从$a^d$开始，我们考虑第一个出现的1或者n-1。
出现了一个n-1，那么后面的都是1。
出现了一个1，而且不是$a^d$，说明有一个不是1也不是n-1的数的平方是1，那么p就不是质数。

所以实际上我们是这么做的：从$a^d$开始，一直平方直到平方了r次或者出现n-1或1。
然后我们可以判断：如果出现了n-1，那么是质数。
如果出现了1，当停止位置是$a^d$的时候是质数。
其他情况都不是质数。

直接上代码：

bool Miller_Rabin(LL x)
{
    if (x==2) return 1;
    if (x<=1 || x&1==0) return 0;
    for (int i=0;i<7 && P_TEST[i]<x;++i)
    {
        LL y=x-1;
        while ((y&1)==0) y>>=1;
        LL t=GetPow(P_TEST[i],y,x);
        while (y!=x-1 && t!=1 && t!=x-1)
        {
            t=QMMul(t,t,x);
            y<<=1;
        }
        if (!(t==x-1 || (y&1) && t==1)) return 0;
    }
    return 1;
}

其中GetPow是快速幂，QMMul是快速模乘（因为x可能是10^18，直接平方会爆）
我测试用的a是2\3\5\7\11\13\17这7个最小的质数，按照Matrix67的说法，$10^{14}$以内都没问题。再大就多用几个质数就好了。比直接随机稳定些。

---

Pollard_Rho

目的：分解质因数。

为了分解质因数，首先要找到一个因数，所以问题实际上就是找到一个因数。

直接随机一个因数显然是不行的。
优化一点，随机一个$\sqrt n$以内的因数，也是不行的。

为什么这两种算法都不行呢？
因为这两种算法实际上就是在一个区间内找一个目标点，一次找到的概率分别是$\frac 1n$和$\frac 1{\sqrt n}$
通过计算可以得出50%正确率需要的随机次数。
然而我不会算，我用程序暴力找了下，大概是区间大小的10倍左右。

这种算法找到的概率低的原因，是因为区间很大，而目标点只有一个。

我们考虑增加目标点。

怎么增加呢？
我们随便加一个点x进去，再随机一个点y，如果|x-y|=p，且p|n的话就刚好找到了一个因数。

在这种情况下，目标点变成了两个：x+p和x-p。

两个还是不够，我们再加：

随便加m个点进去，求任意两个点的差，判断是否是n的因数。

此时，目标点变成了2m个，随机到的概率就很大了。

大致可以这么算：

$$1-\prod_{i=1}^{m-1}\frac {n-2i}n$$
考虑不一定每次都加2n个点进来，简化一点：
$$1-\prod_{i=1}^{m-1}\frac {n-i}n$$
这样大概会在 $\sqrt n$左右得到50%的概率，然后会急剧上升，到 $\sqrt n$的几倍的时候，就可以有99.999%的概率能找到了。

然而放$\sqrt n$个点还不如直接暴力枚举1~$\sqrt n$。

我们考虑一次性加多一点目标点：
随便加n个点进去，求任意两个点的差，判断和n的gcd是否>1。
也就是说，如果我们找到了n的因数的倍数，那么也算。

n肯定有一个$\le\sqrt n$的因数，所以n的因数的倍数至少有$\sqrt n$个。

那么变成：

$$1-\prod_{i=1}^{m-1}\frac {n-i\sqrt n}n$$
大致相当于：
$$1-\prod_{i=1}^{m-1}\frac {\sqrt n-i}{\sqrt n}$$
那么需要放的点数为 $n^{\frac14}$，应该足够少了。

这个东西大概叫做生日悖论，也就是在一个很多人的集体里，同一天生日的概率会很高。两个人同一天生日，相当于他们两出生日期的差是365的倍数。比一个人刚好在某一天生日的概率要高。从维基找的图：

Pollard_Rho算法就是基于这么一个思想，找两个相差为n的因数的倍数的数。
虽然这么说，但是你不能暴力放$n^{\frac 14}$这么多个点然后两两枚举求gcd，这样又变成$O(\sqrt n)$的时间复杂度了。

Pollard_Rho算法实际上每次就随机两个点然后求差.

直接随机两个点求差实际上是很愚蠢的，那就和直接随机一个$\sqrt n$范围内的点差不多。

所以Pollard_Rho算法引入另外一个东西：循环。

假设我们有一个mod n意义下的随机数生成函数g(x)，由种子x可以得到确定的随机数g(x)。
（这个g(x)一般就用$g(x)=x^2+C$即可，C是一个你设的常数。）

从随便一个x开始，我们生成一个随机数序列：x g(x) g(g(x)) g(g(g(x))) ……

因为g(x)是mod n意义下的，所以这个序列肯定会出现循环。

设$n=p\times q$，那么长为p的循环很有可能会比长为n的循环先完成。

在找到循环的同时，我们也就能找到因数。

怎么找循环呢？实际上Pollard_Rho算法用的是Floyd找环法。

设两个变量为x和y，x和y一开始想等，x每次走一步，y每次走两步。
因为y比x快，所以y最终会追上x，然后找到循环。
（追尾的形状就像$\rho$这个字母一样，由此得名。）

算法伪代码如下:

x=y=rand()
d=1
while d=1
    x=g(x)
    y=g(g(y))
    if (x==y) break
    d=gcd(|x-y|,n)

if d>1 找到因数d
else 找不到因数

这个算法不是一定能找到的，要多试几次。

接下来还有问题：
计算次数太多。
每一层循环我们调用3次g(x)函数，这是不优的。

所有有以下优化：
y不动，设一个最大步数k，每当达到最大步数的时候，使y=x。

这样就能跳过一开始的不循环的一段。

最大循环次数怎么设呢？显然不能比环的大小小，也不能太大，直接设并不好。
所以我们不直接设最大循环次数，我们倍增它，一开始k=2，每达到最大步数时，让k翻倍。

这个优化主要优化了计算次数，实际循环次数是增加了的，但是实测加上优化会比较快。

通过这个Pollard_Rho算法，我们就能找到n的一个因数。
要分解n的质因数的话，实际上需要递归多次调用。

设每次找到n的因数p，然后我们递归地分解n/p 和 p。
边界条件是n是质数， 这个用Miller_Rabin判就可以了。

直接上代码：

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>

using namespace std;

typedef long long LL;
typedef unsigned long long ULL;

LL Fac[1000];
int Fn;
const int P_TEST[7]={2,3,5,7,11,13,17};
const int RhoLimit=10000;
const int RhoC=13;
LL GetPow(LL x,LL k,LL MOD);
LL gcd(LL x,LL y);
bool Miller_Rabin(LL x);
LL Pollard_Rho(LL n);
void Work(LL n);
LL Random(LL x);
LL QMMul(LL x,LL y,LL MOD);

int main()
{
    freopen("pollard_rho.in","r",stdin);
    freopen("pollard_rho.out","w",stdout);

    LL n;
    scanf("%I64d",&n);
    Work(n);
    sort(Fac+0,Fac+Fn);
    for (int i=0;i!=Fn;++i) printf("%I64d ",Fac[i]);printf("\n");

    return 0;
}
LL GetPow(LL x,LL k,LL MOD)
{
    LL ans=1;
    for (;k>0;k>>=1,x=QMMul(x,x,MOD))
        if (k&1) ans=QMMul(ans,x,MOD);
    return ans;
}
bool Miller_Rabin(LL x)
{
    if (x==2) return 1;
    if (x<=1 || x&1==0) return 0;
    for (int i=0;i<7 && P_TEST[i]<x;++i)
    {
        LL y=x-1;
        while ((y&1)==0) y>>=1;
        LL t=GetPow(P_TEST[i],y,x);
        while (y!=x-1 && t!=1 && t!=x-1)
        {
            t=QMMul(t,t,x);
            y<<=1;
        }
        if (!(t==x-1 || (y&1) && t==1)) return 0;
    }
    return 1;
}
LL Pollard_Rho(LL n)
{
    LL x=Random(n-2)+2; 
    LL y=x;
    int STEP=2;
    for (int i=1;;++i)
    {
        x=QMMul(x,x,n)+RhoC;
        if (x>=n) x-=n;
        if (x==y) return -1;
        LL d=gcd(abs(x-y),n);
        if (d>1) return d;
        if (i==STEP)
        {
            i=0;
            y=x;
            STEP<<=1;
        }
    }
}
void Work(LL n)
{
    if (Miller_Rabin(n)) {Fac[Fn++]=n;return;}
    LL p;
    for (int i=0;i!=RhoLimit;++i)
    {
        p=Pollard_Rho(n);
        if (p!=-1) break;
    }
    if (p==-1) return;
    Work(p);
    Work(n/p);
}
LL Random(LL x)
{
    unsigned long long p=rand()*rand();
    p*=rand()*rand();p+=rand();
    return p%x;
}
LL gcd(LL x,LL y)
{
    if (y==0) return x;
    return gcd(y,x%y);
}
LL QMMul(LL x,LL y,LL MOD)
{
    LL ans=0;
    for (;y>0;y>>=1)
    {
        if (y&1)
        {
            ans+=x;
            if (ans>=MOD) ans-=MOD;
        }
        x+=x;
        if (x>=MOD) x-=MOD;
    }
    return ans;
}

---

把这两个随机化算法学完之后，发现随机化算法其实还是很靠谱的。前提：要有数学证明。By SemiWaker