Pollard-Rho算法

当需要分解大质数的因子时

在绝大多数问题中，我们所涉及的到的需要找到一个因子的步骤，我们都可以使用试除法来解决

我们很容易想到，对于一个正整数 $N$ ，它的因数是关于 $\sqrt N$ 对称分布的，所以我们只需要在 $[1,\sqrt{N}]$ 的范围内扫一遍，复杂度是 $O(\sqrt{N})$ 的，在 $N$ 较小时，这是比较优秀的，但是当 $N$ 达到了 $1e18$ 级别时，这个方法的时间复杂度是我们无法接受的

对于这样的数，我们有一个可能更为优秀的思路: 随机猜一个因数。那么你的程序时间复杂度下限无疑是 $O(1)$ 的，显然的，对于 $N\ge 1e18$ ，我们得到正确答案的概率是 $\frac{1}{1000000000000000000}$ 。当然，我们如果在 $[1,\sqrt{N}]$ 里进行，概率会更高一些

那么，有没有办法来提高我们的成功率呢

通过悖论来提高成功率

我们来考虑这样的一种情况：在 $[1,1000]$ 的范围内找到一个数，如114 。那么随机选取一个数的成功率显然为 $\frac{1}{1000}$

如果我们选取两个数，让它们的差值为 114 ，即选出 $i,j$ ， $[i-j]=114$ 很显然，这样的概率是 $\frac{1}{500}$ 。那我们如果选多个数并两两取差呢？这即是一个经典的悖论

生日悖论，即有 $k$ 个人，他们生日互不相同的概率是多少。设有 $i$ 个人，它们两生日互不相同的概率为 $P_i$ 那么，可以得出 $P_i=\frac{365}{365}\cdot\frac{364}{365}\cdot\frac{363}{365}\cdots\frac{365-i}{365}\cdot\frac{365-i+1}{365}$ 那么有 $k$ 个人，他们之中存在至少两个生日相同的概率即为 $1-\prod_{i=1}^k\frac{365-i+1}{365}$ 可以大概算一下，当 $k\ge23$ 时，生日相同的概率高于 $50%$ ，而当 $k$ 的取值为 60 时 ，这个概率约为 $0.9999$

由于和人们的日常生活经验并不相吻合，所以这被称之为生日悖论

随机算法的优化

生日悖论的实质即是通过随机采样的方法，满足预想情况的组合比个体更多，来实现较高的成功率

那么对于查找因数，我们该如何选取来提高成功率呢？ 一个很简单的思路即是对于任意正整数 $k$ ,其与 $N$ 的最大公因数一定是 $N$ 的一个因子，即 $\forall k \in N^* , \gcd(k,N)|N$ ，只要我们选取恰当的 $k$ 使得 $\gcd(k,n)>1$ 即可得到 $N$ 的一个因子 $\gcd(k,n)$

于是我们选取了一组数 $x_1,x_2,x_3,x_4,\dots,x_k$ ，若存在 $\gcd([x_i-x_j],N)>1$ ，则称 $\gcd([x_i-x_j],N)$ 是 $N$ 的因子。在某篇论文中有过证明，需要选取的数的个数大约是 $O(N^{\frac{1}{4}})$ ，但是如果我们还需要两两枚举，时间复杂度又变回了 $O(N^{\frac{1}{2}})$

Pollard Rho和它的伪随机数列

我们考虑构造一个伪随机数列，然后取相邻两项来求 $\gcd$。为了生成一个高效的数列，Pollard 设计了这样一个函数：$f(x)=(x^2+c)\mod N$ 其中 $c$ 是一个随机的常数

我们随机取一个 $x_1$ 然后令 $x_2=f(x_1),x_3=f(x_2),\dots,x_i=f(x_{i-1})$ 然后相邻两项做差求 $\gcd$

不过之所以叫做伪随机数列，是因为这个函数是存在循环的，不难理解，这个数列的值域为 $[0,N-1]$ ，至多 $N$ 次就会完成一个循环

基于Floyd优化的 Pollard Rho

为了判断环的存在，我们引入了 Floyd 算法。我们可以是使用一个简单的思想，即龟兔赛跑 ，我们使 t 每次跳一格， r 每次跳两格。假设环长为 len ，当 $x_r=x_t$ 时，因为 $r=2t$ ，所以它们之间的路径差值正好就是 $len$ 的整数倍

inline int f(int x,int c,int n){return (x*x+c)%n;}
inline int Pollard_Rho(int n)
{
    int c=rand()%(n-1)+1;
    int t=f(0,c,n),r=f(f(0,c,n),c,n);
    while(t!=r)
    {
        int d=gcd(abs(r-t),n);
        if(d>1)return d;//找到一个合法因子
        t=f(t,c,n),r=f(f(r,c,n),c,n);
    }
    return n;//没有找到合法因子，重新选取 c
}

通过倍增优化 Pollard Rho

考虑到对于每一个 x，我们都要进行一遍 $O(\log{N})$ 的 $\gcd$ ，这是相当慢的，我们考虑每次选取一段连续区间相乘的差，因为方便和一些玄学原因，我们考虑每次倍增查找质因子，同时如果累计次数超过127次，我们就再统计一次

inline int f(int x,int c,int n){return (x*x+c)%n;}
inline int Pollard_Rho(int n)
{
    int l=0,r=0,c=rand()%(n-1)+1;
    int sum=1;
    for(register int g=1;;g<<=1,l=r,sum=1)
    {
        for(register int step=1;step<=g;++step)
        {
            r=f(r,c,n);
            sum=sum*abs(r-l)%n;
            if(!(step%127))
            {
                int d=gcd(sum,n);
                if(d>1)return d;
            }
        }
        int d=gcd(sum,n);
        if(d>1)return d;
    }
}

至此，我们的 Pollard-Rho 算法已经可以在一个可以接受的复杂度内分解大整数的因子

模板题和其余优化

模板

对于 Pollard-Rho 来说，如果对象是一个大质数，我们的时间复杂度将会异常糟糕。因此，我们一般先用 Miller-Robin素性测试 来判断该数是否为质数，然后再进行分解

#include<stdio.h>
#include<iostream>
#include<math.h>
#include<cstring>
#include<stdlib.h>
#include<time.h>
#include<algorithm>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
ll T;
ll n,ans;

inline ll mul(ull a, ull b, ull p)
{
    ll c=(ll)a*b-(ll)((ull)((long double)a*b/p)* p);
    return c<0?c+p:(c>p?c-p:c);
}

struct Miller_Robin
{
    int p[9]={2,3,7,13,17,61,2333,4567,24251};
    inline ll ksm(ll a,ll b,ll p){ll res=1;while(b){if(b&1)res=mul(res,a,p);b>>=1;a=mul(a,a,p);}return res;}
    inline bool check(ll x,ll p)
    {
        ll k=x-1;
        while(k)
        {
            ll c=ksm(p,k,x);
            if(c^1&&c^(x-1))return false;
            if(k&1||!(c^(x-1)))return true;
            k>>=1;
        }
        return true;
    }
    inline bool isprim(ll x)
    {
        if(x<2)return false;
        for(register int i=0;i<9;++i)
            if(!(x^p[i]))return true;
            else if(!check(x,p[i]))return false;
        return true;
    }
}MR;

struct Pollard_Rho
{
    ll gcd(ll x,ll y){return y?gcd(y,x%y):x;}
    inline ll f(ll x,ll c,ll n){return (mul(x,x,n)+c)%n;}
    inline ll getp(ll n)
    {
        ll c=rand()%(n-1)+1;
        ll l=0,r=0;
        ll sum=1;
        for(register int g=1;;g<<=1,sum=1,l=r)
        {
            for(register int step=1;step<=g;++step)
            {
                r=f(r,c,n);
                sum=mul(sum,abs(r-l),n);
                if(!(step%127))
                {
                    ll d=gcd(sum,n);
                    if(d>1)return d;
                }
            }
            ll d=gcd(sum,n);
            if(d>1)return d;
        }
    }
    void solve(ll x)
    {
        if(x<=ans||x<2)return;
        if(MR.isprim(x)){ans=x;return;}
        ll p=x;
        while(p>=x)
            p=getp(x);
        while(!(x%p))x/=p;
        solve(x),solve(p);
    }
}PR;

signed main(void)
{
    srand((ull)20041022);
    cin>>T;
    while(T--)
    {
        cin>>n;
        ans=0;
        PR.solve(n);
        if(ans==n||!ans)printf("Prime\n");
        else printf("%lld\n",ans);
    }
}