Pollar-Rho算法

Posted by

Pollar-Rho 算法是一种用于快速分解质因数的算法。

问题引入

给定一个正整数N \in \mathbb{N}_{+},试快速找到它的一个因数。

考虑朴素算法,因数是成对分布的,N的所有因数可以被分成两块,即[1,\sqrt N][\sqrt N+1,N]。只需要把[1,\sqrt N]里的数遍历一遍,再根据除法就可以找出至少两个因数了。这个方法的时间复杂度为O(\sqrt N)

N\ge10^{18}时,这个算法的运行时间我们是无法接受的,希望有更优秀的算法。一种想法是通过随机的方法,猜测一个数是不是N的因数,如果运气好可以在O(1)的时间复杂度下求解答案,但是对于N\ge10^{18}的数据,成功猜测的概率是\frac{1}{10^{18}},期望猜测的次数是10^{18}。如果是在[1,\sqrt N]里进行猜测,成功率会大一些。我们希望有方法来优化猜测。

生日悖论

不考虑出生年份,问:一个房间中至少多少人,才能使其中两个人生日相同的概率达到50\%?

解: 假设一年有 n 天,屋子中有 k 人,用整数 1, 2,\dots, k对这些人进行编号。假定每个人的生日均匀分布于 n 天之中,且两个人的生日相互独立。

设 k 个人生日互不相同为事件A, 则事件A的概率为

P(A)=\frac{n}{n} * \frac{n-1}{n} * \dots *\frac{n-k+1}{n}

至少有两个人生日相同的概率为P(\overline A)=1-P(A)。根据题意可知P(\overline A)\ge\frac{1}{2}
,那么就有1 * \frac{n-1}{n} * \dots *\frac{n-k+1}{n} \le \frac{1}{2}

由不等式1+x\le e^x可得

P(A) \le e^{-\frac{1}{n}}* e^{-\frac{2}{n}}*\dots * e^{-\frac{k-1}{n}}=e^{-\frac{k(k-1)}{2n}}\le\frac{1}{2}\
e^{-\frac{k(k-1)}{2n}}\le\frac{1}{2}

n=365代入,解得k=23。所以一个房间中至少23人,使其中两个人生日相同的概率达到50\%,但这个数学事实十分反直觉,故称之为一个悖论。

然而我们可以得到一个不等式方程,e^{-\frac{k(k-1)}{2n}}\le 1-p,其中p是一个概率。那么当k=57n=365时,可以求得p\approx 0.99

考虑一个问题,设置一个数据n,在[1,1000]里随机选取i个数(i=1时就是它自己),使它们之间有两个数的差值为k。当i=1时成功的概率是\frac{1}{1000},当i=2时成功的概率是\frac{1}{500}(考虑绝对值,num_2可以取num_1-knum_1+k),随着i的增大,这个概率也会增大最后趋向于1。

优化随机算法

最大公约数一定是某个数的约数,即\forall k \in\mathbb{N}_{+},gcd(k,n)|n,只要选适当的k使得1<\gcd(k,n)< n,就可以求得一个约数\gcd(k,n)。满足这样条件的k不少,k有若干个质因子,每个质因子及其倍数都是可行的。

根据生日悖论,取多组数据作差能够优化取数的精确性。那么不妨取一组数 x_1,x_2,\dots,x_n若有1<\gcd(|x_i-x_j|,n),则称gcd(|x_i-x_j|,N)N的一个因子(更严谨一些,应该是非平凡因子,即满足 1< d 的那些数)。

构造伪随机函数

构造一个伪随机数序列,取相邻的两项来求\gcd(|x_i-x_j|,N)。通过f(x)=(x^2+c)\%n来生成随机数,其中c=rand(),是一个随机的常数。

随机取一个x_1,令x_2=f(x_1),x_3=f(x_2),\dots,x_i=f(x_{i-1}),在一定范围内可以认为这个数列是“随机”的。

举个例子,设N=50,c=2,x_1=1 f(x)生成的数据为

1,3,11,23,31,11,23,31,\dots

可以发现数据在3以后都在11,23,31之间循环,这也是f(x)被称为伪随机函数的原因。

Floyd判环

考虑两个人在赛跑,A的速度快,B的速度慢,经过一定时间后,A总是会和B相遇,且相遇时A跑过的总距离减去B跑过的总距离一定是圈长的n倍。

初始a=f(1),b=f(f(1)),每一次更新a=f(a),b=f(f(b)),只要检查在更新过程中a、b是否相等,如果相等了,那么就出现了环。

我们每次判断d=gcd(|x_i-x_j|,n),判断d是否满足1< d< n,若满足则可直接返回d。由于x_i是一个伪随机数列,必定会形成环,在形成环时就不能再继续操作了,直接返回n本身,并且在后续操作里调整随机常数c,重新分解。

    ll PR(ll N)
    {
        ll c = rand() % (N - 1) + 1;
        ll t = f(0, c, N);
        ll r = f(f(0, c, N), c, N);
        while(t != r)
        {
            ll d = gcd(abs(t - r), N);
            if(d > 1) return d;
            t = f(t, c, N);
            r = f(f(r, c, N), c, N);
        }
        return N;
    }

倍增优化

使用gcd求解的时间复杂度为O(log N),频繁地调用会使算法运行地很慢,可以通过乘法累积来减少求gcd的次数。如果1< gcd(a,b),则有1< gcd(ac,b)c\in \mathbb{N}_{+},并且有1< gcd(ac \%b,b)=gcd(a,b)

考虑每过一段时间将这些差值进行gcd运算,设s=\prod|x_0-x_j|\%n,如果某一时刻得到s=0那么表示分解失败,退出并返回n本身。每隔2^k个数,计算是否满足1< gcd(s, n)

    ll PR(ll x)
    {
        ll s = 0,t = 0;
        ll c= rand() % (x - 1) + 1;
        int step = 0, goal = 1;
        ll val = 1;
        for( goal=1; ;goal <<= 1, s = t, val = 1)
        {
            for(step = 1; step <= goal; ++step)
            {
                t = f(t, c, x);
                val= val * abs(t - s) % x;
                if( (step % 127) == 0 )
                {
                    ll d = gcd(val,x);
                    if (d > 1) return d;
                }
            }
            ll d = gcd(val, x);
            if(d > 1) return d;
        }
    }

例题:P4718 【模板】Pollard-Rho算法

对于一个数 n ,用Miller Rabin算法判断是否为素数,如果是就可以直接返回了,否则用Pollard-Rho算法找一个因子 p ,将 n 除去因子 p 。再递归分解 n 和 p ,用Miller Rabin判断是否出现质因子,并用max_factor更新就可以求出最大质因子了。由于这个题目的数据过于庞大,用Floyd判环的方法是不够的,这里采用倍增优化的方法。

    #include<bits/stdc++.h>

    using namespace std;

    typedef long long ll;
    #define lll __int128

    int t;
    ll max_factor, n;

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

    inline ll qp(ll x,ll p,ll mod)
    {
        ll ans = 1;
        while(p)
        {
            if(p & 1)
                ans = (lll)ans * x % mod;
            x = (lll) x * x % mod;
            p >>= 1;
        }
        return ans;
    }

    inline bool Miller_rabin(ll p){
        if (p < 2) return 0;
        if (p == 2) return 1;
        if (p == 3) return 1;
        ll d = p - 1, r = 0;
        while(!(d&1)) ++r, d >>=1;
        for(ll k =0; k < 10; ++k) {
            ll a = rand() % (p - 2) + 2;
            ll x = qp(a, d, p);
            if(x == 1 || x == p-1) continue;
            for(int i = 0; i < r - 1; ++ i) {
                x = (lll)x * x % p;
                if(x == p - 1) break;
            }
            if(x != p - 1) return 0;
        }
        return 1;
    }

    inline ll f(ll x,ll c,ll n)
    {
        return ((lll)x * x + c) % n;
    }

    inline ll PR(ll x)
    {
        ll s = 0,t = 0;
        ll c = (ll)rand() % (x - 1) + 1;
        int step = 0, goal = 1;
        ll val = 1;
        for( goal=1; ;goal <<= 1, s = t, val = 1)
        {
            for(step = 1; step <= goal; ++step)
            {
                t = f(t, c, x);
                val= (lll)val * abs(t - s) % x;
                if( (step % 127) == 0 )
                {
                    ll d = gcd(val,x);
                    if (d > 1) return d;
                }
            }
            ll d = gcd(val, x);
            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 = PR(x);
        while((x % p) == 0) x /= p;
        fac(x),fac(p);
    }

    int main()
    {
        scanf("%d", &t);
        while(t --)
        {
            srand((unsigned)time(NULL));
            max_factor=0;
            scanf("%lld", &n);
            fac(n);
            if(max_factor == n) printf("Prime\n");
            else printf("%lld\n",max_factor);
        }
        return 0;
    }

Leave a Reply

电子邮件地址不会被公开。 必填项已用*标注