[学习笔记]Pollard-rho算法

一.什么是Pollard-rho

​ 这是一个用来寻找一个合数的因子的算法。很显然的,我们可以使用试除法，1~\(\sqrt{n}\)之间一个一个试。很显然他很慢。

二.朴素的代码

​ 我们来看一个沙雕代码。

	n=read();
	int a=rand()%(n-2)+2;
	if(n%a==0)cout<<"I found "<<a;
	else cout<<"I'm such a stupid.";

​ 很显然的，这段代码很沙雕，完全就是在碰运气抽签。在这 𝑁 − 1 个数中，我们恰好只有两个因子 𝑝 和 𝑞。因此概率就是 \(\frac{n-1}{2}\), 如果𝑁<=\(10^{10}\)，那么成功的机会约为 \(\frac{1}{𝟓𝟎𝟎𝟎𝟎𝟎𝟎𝟎𝟎𝟎}\),这比购买一张彩票中奖的概率还要小，换句话说，我们不得不使用不同的随机数重复将近 𝑁 次来找到一个因子，显然，这比试除法还要糟。但是我们可以优化他。

三.Birthday-trick

​ Birthday-trick悖论被称为生日悖论。先不讨论他跟生日的关系是什么。但是他可以对上面的沙雕进行优化。

​ 很显然的，我们在\(10^{3}\)中，选一个28的概率就很小。但是如果选两个数的差是28，概率大大滴增加，如果是三个数任意两个差为28呢？k个数呢？

​ 概率大概如下（盗的图）

​ 这张图很直观的告诉我们在\(10^{3}\)下，k=30时，概率就超过一半，k=100时，美滋滋。

四.优化沙雕

​ 对于沙雕程序，我们我们不考虑选取一个单独的a,而是考虑选择两个数的差看看是否满足条件。又Birthday-trick我们看出当我们选的足够多，这个概率就很大了。但是直接除，再加上要选很多的数，这种操作是不行的。再加上答案只有p,q。我们很难找出答案。

​ 但是幸运的是，我们不止这一种方法。对于差为a，我们转向考虑gcd(a,n)>1.对于这个，出答案的可能性会大大增加。至于有多少种答案我们暂不分析。

五.Pollard-rho

​ 很遗憾，我们去要选取大约\(N^{\frac{1}{4}}\)个数字，当N足够大的时候，我们的内存扛不住。于是我们采用Pollard-rho算法。

​ 我们不一口气生成k个数，而是依次生成并进行比较。我们使用一个函数来生成（看起来是的）随机数。

\[𝑓(𝑥) = (𝑥^2 + 𝑎)\mod\;𝑁 \]

​ ( 我们可以自己指定𝑎，也可以用 rand()随即函数来生成，这不是重点 )
​ 我们从𝑥1 = 2 或者其他数开始，让𝑥2 = 𝑓(𝑥1), 𝑥3 = 𝑓(𝑥2), … … 生成规则为\(x_{i+1}=f(x_i)\) .

​ 假设 𝑁 = 55 , 𝑓(𝑥) = ( 𝑥2 + 2 ) 𝑚𝑜𝑑 55

​ 你可以发现对于大部分的数据这个算法能够正常运行，但是对于某些数据，它将会进入无限循环。为什么呢？这是因为存在 𝑓 环的原因。当它发生的时候，我们会在一个有限数集中进行无限循环
​ 例如，我们可以构造一个伪随机函数并生成如下伪随机数：
​ 2, 10, 16, 23, 29, 13, 16, 23, 29, 13 … …
​ 在这个例子中，我们最终将会在16, 23, 29, 13这个圈中无限循环，永远找不到因子

​ 然后我们大致可以将这个过程画为一个\(\rho\),这就是算法名称的来历。

六.关于环

​ 在上面我们发现了算法是可能进入死循环的，为了判定是否出现了循环，我们使用Floyd的方法（当然可以全部记录下来一个个找，不过内存不够）

​ 对于两个由f产生的数a,b来说，我们假设运行一次\(f()\)是走一步的话，我们让a一次走一步，b一次走两步，当b第一次赶上a的时候，就必定出现了环。

七.后记

​ 首先对于这个算法来讲，要特判一下因子为2的数，以及使用Miller-rabin可以提前判断一下是不是素数再进行算法。

八.一道例题

​ luogu的模板题。（超时）

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<fstream>
#include<ctime>
using namespace std;
#define m_for(i,a,b) for(int i=(a);i<=(b);++i)
#define ll long long
const int times=10;
ll multi(ll a,ll b,ll m)
{
    ll ans=0%m;
    while(b)
    {
        if(b&1)ans=(ans+a)%m,b--;
        b >>= 1;
        a=(a+a)%m;
    }
    return ans;
}
ll KSM(ll a,ll b,ll m)
{
    ll ans=1%m;
    while(b)
    {
        if(b&1)ans=multi(ans, a, m),b--;
        b>>=1;
        a=multi(a,a,m);
    }
    return ans;
}
inline bool Miller_Rabin(ll p){
	if(p<=2)return p==2;
	if(!(p&1))return 0;
	ll u=p-1;
	ll power=0;
	while(!(u&1))u>>=1,power++;
	m_for(i,1,times){
		ll a=rand()%(p-2)+2,x=KSM(a,u,p),y;
		for(int i=1;i<=power;++i,x=y){
			if((y=multi(x,x,p))==1&&x!=1&&x!=p-1)return 0;
		}
		if(x!=1)return 0;
	}
	return 1;
}
inline ll m_max(ll a,ll b){return a>b?a:b;}
ll ans=0;
ll mod,c;
ll f(ll x){return (multi(x,x,mod)+c)%mod;}
ll gcd(ll a, ll b){return b ? gcd(b, a % b) : a;}
inline ll m_abs(ll x){return x>0?x:-x;}
ll find(ll p){
	if(!(p&1))return 2;
	if(Miller_Rabin(p))return p;
	mod=p;
	ll a,b;
	a=b=rand()%(p-2)+2,c=rand()%(p-2)+2;
	do{
		a=f(a);b=f(f(b));
		ll x=gcd(m_abs(a-b),p);
		if(x>1&&x<p)return x;
	}while(a!=b);
	return find(p);
}
void Pollard_rho(ll x){
	if(x<=1)return;
	if(Miller_Rabin(x)){
		ans=m_max(x,ans);return ;
	}
	ll u=find(x);
	Pollard_rho(u);
	Pollard_rho(x/u);
}
int main(){
srand(time(NULL));
int cas;
ll p;
scanf("%d",&cas);
while(cas--){
	scanf("%lld",&p);
	if(Miller_Rabin(p))printf("Prime\n");
	else{
		ans=0;
		Pollard_rho(p);
		printf("%lld\n",ans);
	}
}
	return 0;
}