Pollard-Rho算法

prelogue

怎么感觉我这个人和随机化关系这么好。

鲤鱼我是从这篇博客中进行学习的。其中也有很多地方借鉴了大佬的博客，侵权删。

Pollard-Rho 算法

Pollard-Rho 算法是一种求非 1 非自身的因子的高效算法。

main body

我们求因子平常是用的复杂度为 \(O(\sqrt n)\) 的试除法，如果 \(n\) 这个数字很大，我们这个算法就很低效，需要一些卡常技巧来帮助我们卡过去这个题目，甚至是直接就过不去。

但是想从根本来解决问题，我们就得要换一种更加高效的算法 Pollard-Rho 算法。

这个算法是一种带有随机化思想的方法。有点类似于猴子排序（就是那个很 naive 很 trivial 的算法，你能说他效率低嘛，那只能说明你 rp 不好）

我们如果单纯的用这个随机数来进行匹配，效率极其低下，甚至可以说是 \(O(+\infty)\) 的。那我们就要想着怎么对它进行优化。

在这之前，我们要先引入一个知识点生日悖论：

生日悖论

生日悖论是说在一个房间里面有 23 个人，则他们生日相同的概率就会大与 \(\frac{1}{2}\)，反映出来就是这个式子 \(\prod _ {i=0} ^ {22} \frac{365 - i}{365} < \frac{1}{2}\)

这就告诉我们，在一个区间里面不断的生成随机数。很快就会得到重复的数字。期望大约是在根号级别（手玩一下就不难发现，大概就是根号级别的）。

应用到原题上面来，我们随即的最坏情况为 \(n = p^2\)，期望在生成了\(\sqrt n = n ^ {\frac{1}{4}}\) 个数后，可以出现两个在模 \(p\) 的情况下相同的数。那么记这两个数的差的绝对值为 \(d\)，就一定满足 \(d \equiv 0 (\bmod p)\) 那么就同时满足\(gcd(d, n) > 1\)。但是此时我们需要对这些数字一一验证，从整体上来说，时间复杂度并没有更优。

---

Pollard 对于这个问题提出了一个新的解决方案：一种特别的伪随机数生成器用来生成 \([0, N - 1]\) 之间的伪随机数列：设第一个数为 \(x\)，$f(x) = (x ^ 2 + c) \bmod N $， \(c\) 是一个随机参数。

其中，\(x, f(x), f(f(x))\cdots\) 就是一个伪随机数列。

这样子我们通过相邻两项做差，就是我们之前求得的东西，证明如下：

$|a - b| \equiv 0 (\bmod p) \leftrightarrow |f(a) - f(b)| = |a ^ 2 - b ^ 2| = |a - b| \times |a + b| \equiv 0(\bmod p) $

记两个变量分别为 \(t\)，\(r\)。

我们对于这个上这个式子不断的去求，最后由于 \(f(f(r))\) 的增长速度，显然是要快于 \(f(t)\)最终会与 \(t\) 在环上相遇。这样子就形成了一个 \(\rho\) 的形状。

这也就是 Pollard 为什么将这算法叫做 Pollard-Rho 算法。

这个时候如果数据足够随机，它的期望值是 \(O(n ^ {\frac{1}{4}} \log n)\)。

但是实际上，如果 \(gcd(d, n) > 1\)，那么就一定有 \(gcd(kd \bmod n, n) > 1\) 所以我们可以减少求公约数的次数，将这些数字乘起来，然后再统一与 \(n\) 求公因数。

最后反正就是一顿操作猛如虎，然后就给硬生生将一个 \(O(+\infty)\) 的做法砍成了 \(O(n ^ {\frac{1}{4}})\)

这里判读是否为素数我用的是 Miller_Rbin，不会的可以自行递归学习。（

code time

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define rl register ll
#define fom(i, a) for(rl i=a; i; -- i)
#define foa(i, a, b) for(rl i=a; i < b; ++ i)
#define fos(i, a, b) for(rl i=a; i <= b; ++ i)
#define fop(i, a, b) for(rl i=a; i >= b; -- i)

namespace IO
{
    int pp1=-1; const int pp2=(1<<21)-1; char buf[1<<21],*p1=buf,*p2=buf,buffer[1<<21];
    inline void flush() {fwrite(buffer,1,pp1+1,stdout); pp1=-1;}
    inline void pc(const char ch) {if(pp1==pp2) flush(); buffer[++pp1]=ch;}
    inline char gc(){ return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
    template <class T> inline void read(T &res){char ch=gc();bool f=0;res=0; for(;!isdigit(ch);ch=gc()) f|=ch=='-'; for(;isdigit(ch);ch=gc()) res=(res<<1)+(res<<3)+(ch^48); res=f?~res+1:res;}
    template <class T> inline void ww(T x) { if(!x) pc('0'); else { static int stk[21]; int top = 0; if(x<0) pc('-'),x=~x+1; while(x) stk[top++]=x%10, x/=10; while(top--) pc(stk[top]^48);}}
}

#define ws IO::pc(' ')
#define wl IO::pc('\n')
#define ww(x) IO::ww(x)
#define read(x) IO::read(x)
#define flush() IO::flush()

ll max_factor, n;

// mt19937 rand(time(0));

inline ll gcd(ll a, ll b) { return b ? gcd(b, a % b) : a; } // 求最大公因数

inline ll qmi(ll a, ll k, ll p) { ll res = 1; while(k) { if(k & 1) res = ((__int128)res * a) % p; a = ((__int128)a * a) % p, k >>= 1; } return res; } // 被压缩的快速幂

inline bool Miller_Rabin(ll p) // 判断是否为素数
{
    if(p < 2) return false;
    if(p == 2 || p == 3) return true;
    ll d = p - 1, r = 0;

    while(!(d & 1)) ++ r, d >>= 1;
    
    foa(k, 0, 10) 
    {
        ll a = rand() % (p - 2) + 2;
        ll x = qmi(a, d, p);
        if(x == 1 || x == p - 1) continue;
        foa(i, 0, r - 1)
        {
            x = (__int128)x * x % p;
            if(x == p - 1) break;
        }
        if(x != p - 1) return false;
    }
    return true;
}

inline ll Pollard_Rho(ll n) // 本题的核心算法
{
    ll s = 0, t = 0;
    ll c = (ll)rand() % (n - 1) + 1;
    ll step = 0, goal = 1, val = 1;

    for(goal = 1;; goal *= 2, s = t, val = 1)
    {
        for(step = 1; step <= goal; ++ step)
        {
            t = ((__int128)t * t + c) % n;
            val = ((__int128)val * abs(t - s)) % n;
            if((step % 127) == 0) { ll d = gcd(val, n); if(d > 1) return d; }
        }
        ll d = gcd(val, n);
        if(d > 1) return d;
    }
}

inline void fac(ll x) // 分解
{
    if(x <= max_factor || x < 2) return ;
    if(Miller_Rabin(x)) { max_factor = max(max_factor, x); return ; }

    ll p = x;
    while(p >= x) p = Pollard_Rho(x);
    while((x % p) == 0) x /= p;
    fac(x), fac(p); 
}

inline void solve() // 多测
{
    srand(time(0));
    max_factor = 0; read(n);
    fac(n);
    if(max_factor == n) flush(), puts("Prime");
    else ww(max_factor), wl;
}

int main()
{
    // freopen("1.in", "r", stdin);
    ll T; read(T); while(T --) solve(); 
    flush(); return 0;
}