数论 Pollard_Rho算法

Pollard_Rho算法解决大数的质因数分解。又是一个玄学算法..

我们的任务是对一个数字n进行质因数分解。可以发现，n的因数将会对称的分布在[1,sqrt(n)],和[sqrt(n),n]两个区间中，我们只需对前者扫一遍，即可求出所有的质因数，复杂度为\(0(\sqrt(n))\)

但是，假如数据范围是longlong范围呢？1e18呢？试除法似乎有点捉襟见肘…

emm，很明显，完整的扫一遍是不可能的，尝试随机一下。假设数字为在1e18范围内充斥着两个质因数，我们随机来一发，会有多大的几率中奖？的确，我们有一定的概率会得到一个因数，但概率显而易见不在我们可以接受的范围内。

生日悖论，指如果一个房间里有23个或23个以上的人，那么至少有两个人的生日相同的概率要大于50%。这就意味着在一个典型的标准小学班级(30人)中，存在两人生日相同的可能性更高。对于60或者更多的人，这种概率要大于99%。

这个生日悖论，条件显而易见，证明方法显而易见，只是因为与常识有所违背，所以称之为"悖论"。

这个悖论告诉我们的精华在于：在[1,n]内的数字中随机抽取数字，使之与指定数字p相同，抽取\(\sqrt(n)\)次便有了50%的概率；抽取次数越高，概率越可观。即，在范围内抽出数字使之符合期望是一件可以实现的事情。

这样子看来，虽然随机一遍成功的概率微乎其微，但通过多次随机，其期望步骤在可接受范围内。

但是，眼前还有很多问题没有解决

[1] 怎么随机出合适的数字？直接随机？

[2] 怎么合理的判断是否是因数？

对于上述问题算法，Pollard_Rho算法提出了一个切实可行的方法。

引出一个函数如下:

\[f(n)=(n^2+c) \quad mod\ x
\]

c是一个随机的常数，范围在[1,x-1]；x是当前进行质因数分解的数字。

从0开始，代入函数，并将结果再次代入函数，会顺序地得到一个数列a[..]。由于函数的特性，在一定范围内可以保证数列的随机性。数列中相邻两个数字做差，就可以得到一个随机数字n。求d=gcd(n,x)，倘若d大于1，d便是x的一个因数。

由于函数的某种玄学特性，这个做法相比普通随机可以节省大量时间。具体证明就..鸽了..╰(￣ω￣ｏ)

此时，我们得到了x的一个因数d，但这显然是不够的。我们需要将d和x/d继续分解，直到全部分解为质因数。

你知道为什么这个算法叫\(Pollard \quad \rho\)吗？（罗马字母\(\rho\)发音rho）因为通过这个函数得到数列可能会存在循环节，倘若出现了循环节，我们将反反复复进行大量无用尝试，而无法得到结果。因此，我们需要一个有效的判环操作。

通过生日悖论可知，这个环的期望长度在\(\sqrt n\)级别，但是我还是不会证明

floyd判环的思想在于设两个指针，让两个指针同时以不同的速率进行运动，这样的话倘若两个指针进入了环，期望只需要\(\sqrt n\)步，两个指针便会再次重逢，此时退出循环。枚举不同的c，并再次尝试。

既然我们已经有了两个指针了，就不再拘泥于必须要相邻的两个数字。两个指针的结果作差取绝对值即可。

至此，Pollard_Rho算法的基础理论，流程框架已经完整地展现出来了。

[1] 先用Millar_Rabin判断是否为质数，如果是质数说明已经到头；否则进行接下来的操作。

[2] 用函数生成数列，得到随机数n，在gcd(n,x)>1时，得到一个因数gcd(n,x)

[3] 倘若因为不合适的c导致了数列出现了循环节，利用floyd判环，及时退出即可

[42] 将x拆分为x和x/gcd(n,x)，递归下去

求解一次gcd的复杂度是\(o(log\ n)\)，大量求解gcd必然会降低效率。因此我们要慎重使用gcd

倘若使用频率过低，可能会出现出现了答案却不能及时退出；使用频率过高，效率就会显著受到影响

那么，就倍增路径！就是这么直接( ¯(∞)¯ )

每\(2^i\)步判断一次gcd，在此之前，递乘每一次的结果。由于某种特殊神秘的性质，对结果取模并不会对gcd产生影响。因此可以放心地取模

这个，从洛谷上看到的。尝试了一下发现确实虽然不愿意承认有实际效果。（ac掉最后一个tle点）

方法是：每127步求一遍gcd。

有点难以解释..127除了是质数以外，似乎没有其它特殊的性质。就当作一个玄学优化技巧吧。

在洛谷上有Pollard_Rho的板子题，代码如下：

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;

const int MAX=1e6+5,TOP=1e4;
bool vis[MAX]; int prime[MAX];
ll max_factor;

int T; ll n;

ll read();
void oula();
ll qsm(ll a,ll k,ll m);
ll qsc(ll a,ll b,ll m);
bool millar_rabin(ll);
void pollard_rho(ll);
ll frac(ll);
ll func(ll,ll,ll);
ll gcd(ll a,ll b);

int main(){
    #ifndef ONLINE_JUDGE
    freopen("test.in","r",stdin);
    #endif

    oula();

    scanf("%d",&T);
    while(T--){
        n=read();
        max_factor=0;
        pollard_rho(n);
        if(max_factor==n) printf("Prime\n");
        else printf("%lld\n",max_factor);
    }

    return 0;
}

ll read(){
    char tmp=getchar(); ll sum=0; bool flag=false;
    while(tmp<'0'||tmp>'9'){
        if(tmp=='-') flag=true;
        tmp=getchar();
    }
    while(tmp>='0'&&tmp<='9'){
        sum=(sum<<1)+(sum<<3)+tmp-'0';
        tmp=getchar();
    }
    return flag?-sum:sum;
}

void oula(){
    for(int i=2;i<=TOP;++i){
        if(!vis[i]) prime[++prime[0]]=i;
        for(int j=1;j<=prime[0]&&prime[j]*i<=TOP;++j){
            vis[prime[j]*i]=true;
            if(i%prime[j]==0) break;
        }
    }
}

ll qsm(ll a,ll k,ll m){
    ll ans=1,base=a,mod=m;
    while(k){
        // if(k&1) ans=ans*base%mod;
        if(k&1) ans=qsc(ans,base,mod);
        // base=base*base%mod;
        base=qsc(base,base,mod);
        k>>=1;
    }
    return ans;
}

ll qsc(ll a,ll b,ll m){
    // ll ans=0,base=a,mod=m;
    // while(b){
    //  if(b&1) ans=(ans+base)%mod;
    //  base=(base+base)%mod;
    //  b>>=1;
    // }
    // return ans;

    ll z=(ld)a/m*b;
    ll res=(ull)a*b-(ull)z*m;
    return (res+m)%m;

    return ((ull)a*b-(ull)((ld)a/m*b*b))%m;
}

bool millar_rabin(ll x){
    if(x==2) return true;
    if(!(x&1)||x==0||x==1) return false;
    ll s=0,t=x-1;
    while(!(t&1)) s++,t>>=1;
    for(int p=1;p<=30&&prime[p]<x;++p){
        ll b=qsm((ll)prime[p],t,x),k;
        for(ll i=1;i<=s;++i){
            k=qsc(b,b,x);
            if(k==1&&b!=1&&b!=x-1) return false;
            b=k;
        }
        if(b!=1) return false;
    }
    return true;
}

void pollard_rho(ll x){
    if(x==1) {max_factor=1; return;}
    if(millar_rabin(x)) {max_factor=max(x,max_factor); return;}
    ll p;
    do{
        p=frac(x);
    }while(p>=x);
    while(x%p==0) x/=p;
    if(x!=1) pollard_rho(x); pollard_rho(p);
}

ll frac(ll x){
    ll s,t,c,d;
    c=rand()%(x-1)+1;
    s=0; t=func(0,c,x);
    ll goal=1,step=0,val=1;
    while(s!=t){
        step++;
        val=qsc(val,abs(s-t),x);
        if(step==goal){
            d=gcd(val,x);
            if(d>1) return d;
            goal<<=1ll;
        }
        s=func(s,c,x);
        t=func(func(t,c,x),c,x);
    }
    return x;
}

ll func(ll x,ll c,ll mod){
    return (qsc(x,x,mod)+c)%mod;
}

ll gcd(ll a,ll b){
    return b==0?a:gcd(b,a%b);
}

一个非常有用的pdf