专栏文章

PollardRho-快速分解质因数

zzyzzyzzyzzyz
算法·理论参与者 29已保存评论 34

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
34 条
当前快照
1 份
快照标识符
@mhz5rxr9
此快照首次捕获于
2025/11/15 01:55
4 个月前
此快照最后确认于
2025/11/29 05:24
3 个月前
查看原文
鸣谢rvalue对此题解提出的宝贵意见

Pollard_Rho

PollardRhoPollard Rho (在此简称PR)可以用来在 O(N14)O(N^{\frac{1}{4}}) 的时间内分解质因数.
(这个算法是PollardPollard提出来的;算法中会涉及到一个环,它的形状为′′ρ′′''\rho'',所以叫PollardRhoPollard Rho )

题面(可以点)

求一个数的最大质因数.
这题不需要卡常,不需要卡常,不需要卡常!!!

前置知识:

Miller_Rabin,快速幂,快速乘,gcd

1.Miller_Rabin

Miller_Rabin(在此简称MR)是一个高效(O(log2N)O(log_2 N))判素数的随机算法,是PR的一部分,十分重要
而MR也有两个前置知识:费马小定理,二次探测定理

(1)费马小定理

这个应该比较简单吧.
若 pp 为质数,有
ap≡a(modp)a^p \equiv a\pmod{p}
我们现在要验证的数为pp ,那么可以选取一些数作为aa,带入上式,分别检验是否满足费马小定理.
只要有一个aa 不满足 ap≡a(modp)a^p \equiv a\pmod{p} ,那么p为合数.
但是这只是必要不充分条件.存在这样一类强伪素数pp,满足
∀a<p,ap≡a(modp)\forall a<p, a^p\equiv a\pmod{p}
这也就意味着无法使用费马小定理将它判断为合数.

(2)二次探测

若x2≡1(modp)x^2\equiv1\pmod{p},且pp为质数,则x≡1(modp)x\equiv1\pmod{p}或x≡p−1(modp)x\equiv p-1\pmod{p}
证明:
因为x2≡1(modp)x^2\equiv1\pmod{p}
所以p∣(x−1)(x+1)p\mid(x-1)(x+1)
所以p∣x−1orp∣x+1p\mid x-1 \quad or\quad p\mid x+1
所以x≡1(modp)x\equiv1\pmod{p}或x≡p−1(modp)x\equiv p-1\pmod{p}
证毕.
那么要如何利用它呢?
我们要检测的数仍然为pp,然后再选定一个数xx.
此时pp已经满足了费马小定理(否则会被直接判合数),即:xp−1≡1(modp)x^{p-1} \equiv 1\pmod{p}
(注意,这个式子中的模数在递归求解的过程中是始终不变的,但指数会改变)
这个式子的形式是不是和二次探测定理中的形式很相似?
实际上,我们可以利用二次探测来继续判断质数.
若2∣p−12\mid p-1 ,则(xp−12)2≡1(modp)(x^\frac{p-1}{2})^2\equiv1\pmod{p},讨论xp−12x^\frac{p-1}{2} 模pp意义下的值.
a. 若既不是1,也不是p-1,那么说明p是合数,返回false.
b. 若是1,则继续递归 xp−12x^\frac{p-1}{2}
c.若为p-1,那么不能利用二次探测继续递归,说明目前无法验证p为合数,返回true.
多取几个x来测试正确率就很高了.
Code:
CPP
IL int qpow(int x,int a,int mod)//快速幂
{
    int ans = 1;
    while(a)
    {
        if(a&1) ans=ans*x%mod;
        a>>=1,x=x*x%mod;
    }
    return ans;
}
IL bool test(int x,int p)//费马小定理
{
    return qpow(x,p-1,p)==1;
}
IL bool Miller_Rabin(int x,int p)
{
    if(!test(x,p)) return false;
    int k=p-1;
    while(!(k&1))
    {
        k>>=1;
        RG int t=qpow(x,k,p);
        if(t!=1&&t!=p-1) return false;
        if(t==p-1) return true;
    }
    return true;
}
IL bool Test(int p)
{
    if(p==1) return false;
    if(p==2||p==3||p==5||p==7||p==11) return true;
    return Miller_Rabin(2,p)&&Miller_Rabin(3,p)&&Miller_Rabin(5,p)&&Miller_Rabin(7,p);
    //取不同的x来测试,提高正确率
}

2.快速幂,快速乘

前者就不说了,主要是快速乘.(不过一般快速幂里面的乘法也要用到快速乘)
它的用处是计算两个longlong级别的数的模意义下的乘法
原理:
a%b=a−[a/b]×ba\%b=a-[a/b]\times b
网上搜O(1)快速乘,这里不解释了.
Code:
CPP
IL int qmul(int x,int y,int mod)
{
    return (x*y-(long long)((long double)x/mod*y)*mod+mod)%mod;
}

3.gcd

我用的是辗转相除,但听说用二进制更快?但是好像也只是常数级别的优化(后文会提到一个很重要的东东)
CPP
IL int gcd(int x,int y)
{
    return y?gcd(y,x%y):x;	
}
了解了上述知识后,你就可以开始做这题了QAQ

做法

方法一:试除法

这个很简单,就不说了.

方法二:rand

我要分解的数为NN,我在区间[1,N][1,N]中rand一个数,看它是不是N的因数.
(很搞笑吧)

方法三:Birthday Trick 优化的rand(正解)

嗯,我们从[1,N][1,N]中rand两个数,那么它们的差为N的因数的可能性会大大提升.
但是因为要两两作差,所以没有优化.
但是我们可以退而求其次,使这两个数的差不一定要满足与N的因数完全相等,而是它们的最大公因数不等于一.
那么这个时候我们成功的几率就很高了.
原因大概如下 （此处为感性理解，具体解释看后文）:
若一个数N=p∗qN=p*q (p,q均为质数)
那么满足(N,x)≠1(x<N)(N,x) \ne 1(x<N)的个数是p+q−2p+q-2.
所以找一次找出p,qp,q 成功的概率为p+q−1N\frac{p+q-1}{N}.
推广到求N=∏i=1nprimeiaiN=\prod_{i=1}^{n}prime_i^{a_i}的非平凡因子(不是1与本身的因数),找ii次能找到的概率.
pi=1−∏j=0iϕ(N)−jNp_i=1-\prod_{j=0}^i \frac{\phi(N)-j}{N}
-----口胡结束-----
如果这样做,要选N14N^{\frac{1}{4}}个数,无法存储.
因此,我们必须每次运算只记录两个数.
我们构造一个伪随机数,推荐以下代码
CPP
 t1=(qmul(t1,t1,x)+b);
就是xi=xi−12+Cx_i=x_{i-1}^2+C(C为常数).
比较xx数列的相邻两项.
但是,会出现环.因为我们的xx数列是模意义下生成的,所以可能∃j≠i,xi=xj\exists j\ne i,x_i=x_j.
那么如何解决呢?
用hash吗?不不不,那太麻烦了.
以下内容为引用
那么,如何探测环的出现呢? 一种方法是记录当前产生过的所有的数x1,x2,......xnx_1 , x_2 , ... ... x_n ,并检测是否存在xl=xn(l<n)x_l = x_n (l < n)。 在实际过程中,当 n 增长到一定大小时,可能会造成的内存不够用的情况。 另一种方法是由Floyd发明的一个算法,我们举例来说明这个聪明而又有趣的算法。 假设我们在一个很长很长的圆形轨道上行走,我们如何知道我们已经走完了一圈呢? 当然,我们可以像第一种方法那样做,但是更聪明的方法是让 A 和 B 按照 B 的速度是 A 的速度的两倍从同一起点开始往前走,当 B 第一次敢上 A 时(也就是我们常说的套圈) , 我们便知道,B 已经走了至少一圈了。
while ( b != a )
a = f(a); // a runs once
b = f(f(b)); // b runs twice as fast.
p = GCD( | b - a | , N);
就这样,我们可以较快的找出NN的两个非平凡因子p,qp,q.
递归求解,直到本身为素数返回即可.
代码如下:
CPP

void Pollard_Rho(int x) 
{
    if(Test(x))//素数测试
    {
        Ans=max(x,Ans);
        return; 
    }
    int t1=rand()%(x-1)+1;
    int t2=t1,b=rand()%(x-1)+1;
    t2=(qmul(t2,t2,x)+b)%x;//要用快速乘
    int i=0;
    while(t1!=t2)
    {
        ++i;
        int t=gcd(abs(t1-t2),x);
        if(t!=1&&t!=x) 
        {
            Pollard_Rho(t),Pollard_Rho(x/t);	
        }
        t1=(qmul(t1,t1,x)+b)%x;
        t2=(qmul(t2,t2,x)+b)%x;
        t2=(qmul(t2,t2,x)+b)%x;
    }
}

正解优化

我在这里提一个比较重要的优化,加上后基本不需要卡常就可以跑进2500ms2500ms.
我们要频繁地求gcd,可不可以更快地求呢?
可以!
我们可以将若干个差值累积乘到一起,再求gcd.这是不影响正确性的.
因为若gcd(a,N)>1gcd(a,N)>1,那么gcd(a×b,N)>1gcd(a\times b,N)>1
这样就可以少求很多gcd.
至于累乘多少个差值,这就比较玄学了.这题取的是127(可能是面向数据编程?QAQ)
改进代码如下:
CPP
void Pollard_Rho(int x) 
{
    if(Test(x))
    {
        Ans=max(x,Ans);
        return; 
    }
    int t1=rand()%(x-1)+1;
    int t2=t1,b=rand()%(x-1)+1;
    t2=(qmul(t2,t2,x)+b)%x;
    int p=1;
    int i=0;
    while(t1!=t2)
    {
        ++i;
        p=qmul(p,abs(t1-t2),x);
        if(p==0) //①
        {
            int t=gcd(abs(t1-t2),x);
            if(t!=1&&t!=x) 
            {
                Pollard_Rho(t),Pollard_Rho(x/t);	
            }
            return;
        }
        if(i%127==0)
        {
            p=gcd(p,x);
            if(p!=1&&p!=x)
        	{
            	Pollard_Rho(p),Pollard_Rho(x/p);
            	return;
        	}
            p=1;
    	}
        t1=(qmul(t1,t1,x)+b)%x;
        t2=(qmul(t2,t2,x)+b)%x;
        t2=(qmul(t2,t2,x)+b)%x;
    }
    p=gcd(p,x);
    if(p!=1&&p!=x)
    {
        Pollard_Rho(p),Pollard_Rho(x/p);
        return;//②
    }
}
还是有两个要注意的地方:
①:p为0,说明乘上当前差值后变为了x的倍数,那么当前差值与N的gcd一定为x的因子.
②:ρ\rho环的长度可能小于127,所以需要在跳出循环时判断.

Code:

请自行忽略#define int long long
CPP
#include<bits/stdc++.h>
#define gc getchar
#define R register int
#define RG register
#define int long long
#define IL inline 
using namespace std;
int rd()
{
    int ans = 0,flag = 1;
    char ch = gc();
    while((ch>'9'||ch<'0')&&ch!='-') ch=gc();
    if(ch == '-') flag=-1;
    while(ch>='0'&&ch<='9') ans+=(ans<<2)+(ans<<1)+ch-48,ch=gc();
    return flag*ans;
}
int Ans,n;

IL int qmul(int x,int y,int mod)
{
    return (x*y-(long long)((long double)x/mod*y)*mod+mod)%mod;
}
IL int qpow(int x,int a,int mod)
{
    RG int ans = 1;
    while(a)
    {
        if(a&1) ans=qmul(ans,x,mod)%mod;
        a>>=1,x=qmul(x,x,mod)%mod;
    }
    return ans;
}
IL bool test(int x,int p)
{
    return qpow(x,p-1,p)==1;
}
IL bool Miller_Rabin(int x,int p)
{
    if(!test(x,p)) return false;
    int k=p-1;
    while(!(k&1))
    {
        k>>=1;
        RG int t=qpow(x,k,p);
        if(t!=1&&t!=p-1) return false;
        if(t==p-1) return true;
    }
    return true;
}
IL bool Test(int p)
{
    if(p==1||p==2152302898747) return false;
    if(p==2||p==3||p==5||p==7||p==11) return true;
    return Miller_Rabin(2,p)&&Miller_Rabin(3,p)&&Miller_Rabin(5,p)&&Miller_Rabin(7,p)&&Miller_Rabin(11,p);
}

inline int gcd(int x,int y)
{
    return y?gcd(y,x%y):x;  
}

void Pollard_Rho(int x) 
{
    if(Test(x))
    {
        Ans=max(x,Ans);
        return; 
    }
    int t1=rand()%(x-1)+1;
    int t2=t1,b=rand()%(x-1)+1;
    t2=(qmul(t2,t2,x)+b)%x;
    int p=1;
    int i=0;
    while(t1!=t2)
    {
        ++i;
        p=qmul(p,abs(t1-t2),x);
        if(p==0) 
        {
            int t=gcd(abs(t1-t2),x);
            if(t!=1&&t!=x) 
            {
                Pollard_Rho(t),Pollard_Rho(x/t);    
            }
            return;
        }
        if(i%127==0)
        {
            p=gcd(p,x);
            if(p!=1&&p!=x)
            {
                Pollard_Rho(p),Pollard_Rho(x/p);
                return;
            }
            p=1;
        }
        t1=(qmul(t1,t1,x)+b)%x;
        t2=(qmul(t2,t2,x)+b)%x;
        t2=(qmul(t2,t2,x)+b)%x;
    }
    p=gcd(p,x);
    if(p!=1&&p!=x)
    {
        Pollard_Rho(p),Pollard_Rho(x/p);
        return;
    }

}
signed main()
{
    //freopen("Divide.in","r",stdin);
    //freopen("Divide.out","w",stdout);
    ios::sync_with_stdio(false);
    cin>>n;
    for(R i=1;i<=n;i++)
    {
        RG int t;
        cin>>t;
        if(Test(t))
        {
            cout<<"Prime"<<endl;
            continue;
        }
        Ans=0;
        while(Ans==0)
            Pollard_Rho(t);
        cout<<Ans<<endl;
    }
    return 0;
}
基本没卡常.

复杂度

PollardRhoPollardRho 的复杂度是O(N14)O(N^\frac{1}{4}).
令N=a×bN=a\times b,显然a<Na<\sqrt N.
可以想到,我们在作差求gcd时,其实是在找模a意义下相等,在模N意义下不等的两个数.
那么我们的"天数"就是a,而要找到"相同生日的两人",只需a\sqrt a人概率就很大了.
所以,我们的算法复杂度为O(N4)O(\sqrt[4]{N}).

生日悖论的复杂度

首先是生日悖论. 一个 23 人的团体存在两人生日相同的概率要大于 50%. 一个推论是在 n 个值中随机选取若干个值, O(n)O(\sqrt{n})
次后就会有很大概率产生某个值与之前选的值重复的情况.
大概证明如下:
我们补集转化一下, 转而求选出的 k 个值都不同的概率. 显然它应该长这样: Pn(k)=∏i=0k−1(1−in)P_n(k)=\prod\limits_{i=0}^{k-1}(1-\frac{i}{n})
我们只要让这个值小于 12 就好了. 而由泰勒展开可得: exp(x)=1+x+x2/2!+x3/3!+⋯exp(x)=1+x+x^2/2!+x^3/3!+⋯
那么对于 x>0 有: 1+x<exp(x)1+x<exp(x)
于是就有: Pn(k)=∏i=0k−1(1−i/n)<∏i=0k−1exp(−i/n)=exp(−∑i=0k−1in)=exp(−k(k−1)2n)P_n(k)=\prod\limits_{i=0}^{k-1}(1-i/n)<\prod\limits_{i=0}^{k-1}exp(-i/n)=exp(-\sum\limits_{i=0}^{k-1}\frac{i}{n})=exp(\frac{-k(k-1)}{2n})
那么我们只要让不等式右边小于 12\frac{1}{2} 就好了. 那么我们有:
exp(−k(k−1)2n)<12exp(\frac{-k(k-1)}{2n}) < \frac{1}{2}
两边取对数解一下就有:
k2−k>2nln⁡2k^2-k>2n\ln2
又因为 ln2 是个常数, 于是k=O(k)k=O(\sqrt{k})
.
本文参考了*A Quick Tutorial on Pollard's Rho Algorithm(应该打不开)*,前文的翻译版,百度百科-生日悖论

评论

共 34 条评论，欢迎与作者交流。

正在加载评论...