Pollard-Rho算法

prelogue

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

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

Pollard-Rho 算法

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

main body

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

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

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

我们如果单纯的用这个随机数来进行匹配，效率极其低下，甚至可以说是 $O(+\mathrm{\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(modp)$ 那么就同时满足$gcd(d,n)>1$。但是此时我们需要对这些数字一一验证，从整体上来说，时间复杂度并没有更优。

---

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

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

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

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

记两个变量分别为 $t$，$r$。

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

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

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

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

最后反正就是一顿操作猛如虎，然后就给硬生生将一个 $O(+\mathrm{\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;
}

__EOF__

本文作者：carp
本文链接：https://www.cnblogs.com/carp-oier/p/Pollard_Rho.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！