Pollard_Rho 大数求因子

目的

求合数 n 的某个非平凡因子（除 1 和它本身以外能够整除它的数）。

即 快速找到 n 的一个非1非本身的因子。

时间复杂度

期望时间复杂度 $$O(N^{\frac{1}{4}})$$

另一个时间复杂度 $$O(\sqrt{p})$$

p 为 n 的某个最小因子，满足 p 与 n/p 互质。

原理

式子

存在x, c, y, p：

p 为给出的数，要求的是这个数的一个非平凡因子

x = 随机数

c = 随机数

$$y=x,x=(x\times x+c)\%p,g=gcd(abs(y-x),p)$$

g > 1，说明 g 为 p 的一个非平凡因子

这是一个随机的方法，出现 g > 1的概率很高。

这恐怕没什么好说的，因为概率高就拿来用了。

对于取模后出现循环的处理

$$y=x,x=(x\times x+c)\%p$$

这个生成的数可能会重复，最后出现一个环。

我学习的是倍增的思路，y = x 的操作在 x 运算 2 次，4 次，8 次，16 次……的时候后再执行。

y 储存的是环上的一个数，x 就一直在环上跑，如果 y == x 的话说明已经跑过了一个圈，但是也不会跑过太多。

$$\sum _{i=1}^{n-1}{2}^i<2^n$$

对 gcd 的优化

gcd 时间复杂度是 $$O(logN)$$，N 为所求的数。

如果每次运行这个式子都来一个 gcd 那么时间复杂度就要乘上一个 log。

优化方法：

存在 n，m

$$n>m$$

$$gcd(n,m)=x$$

$$nmodm\equiv 0(modx)$$

所以可以将 x 乘起来，结果肯定保留着 p 的因子。

一般乘 127 次后取一次 gcd，但有可能 x 在这之前成为了 p 的倍数，取模后等于 0。

一次还可以在每次倍增后取一次 gcd，这样时间复杂度就优化了。

代码

inline ll quick_multi(ull x, ull y, ull p)		// x，y，p 范围 1e18
{
    return (x * y - (ull)((lb) x / p * y) * p + p) % p;
}

inline ll Abs(ll x)
{
    return x < 0 ? -x : x;
}

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

inline ll rho(ll p)             //求出p的非平凡因子
{  
    ll x, y, z, c, g;
    int i, j;                 //先摆出来（z用来存（y-x）的乘积）
    while (1)                   //保证一定求出一个因子来
    {              
        y = x = rand() % p;                //随机初始化
        z = 1;
        c = rand() % p;                  //初始化
        i = 0, j = 1;                    //倍增初始化
        while (++i)                      //开始玄学生成
        {                    
            x = (quick_multi(x, x, p) + c) % p;  //可能要用快速乘
            z = quick_multi(z, Abs(y - x), p);  //我们将每一次的（y-x）都累乘起来
            if (x == y || !z) break;            //如果跑完了环就再换一组试试（注意当z=0时，继续下去是没意义的）
            if (!(i % 127) || i == j)           //我们不仅在等127次之后gcd我们还会倍增的来gcd
            {  
                g = gcd(z, p);
                if (g > 1) return g;
                if (i == j)
                    y = x, j <<= 1;  //维护倍增正确性，并判环（一箭双雕）
            }
        }
    }
}

inline ll prho(ll x)        // 求最大质因子
{
    if (x == 1) return 1;
    else if (prime(x)) return x;    // mille rabin 素数判定
    else 
    {
        ll p = rho(x);
        while (x % p == 0) x/= p;
        return max(prho(x), prho(p));
    }
}

备注

我根据这篇博客学的【快速因数分解】Pollard’s Rho 算法 - Koshkaaa (cnblogs.com)