Pollard_Rho大数分解学习小记

前言

• 该算法用于对大数进行质因子分解。
• 它的核心思想是：对于一个数N，若N不为质数，设N=c*d（c、d均大于1），我们找到一个c，递归分解c和d。

生日悖论

• 一个班有k个同学，至少有两名学生的生日是同一天的概率P是多少？
• 正难则反，考虑容斥。P=1-任意两名学生的生日都不同的概率。
• $k=1,P=0$
• $k=2,P=1-1\ast \frac{364}{365}$
• $k=3,P=1-1\ast \frac{364}{365}\ast \frac{363}{365}$
• $\dots \dots$
• $k=23,P=0.507$
• 当k>55时，P就很接近于1了。
• 当k=100时，P=0.999999692751072。
• 其实，如果说一年有N天，那么只要$k\ge \sqrt{N}$，他们中存在相同生日的概率就会大于50%。（详见《算导》黑书）

这有什么关联吗

• 估计大家刚看到生日悖论会和我产生同样的疑问：这有什么关联吗？
• 实际上，关联很大。
• 我们从[2,N]范围内随机k个数，刚好刷出c或d的概率很小，但如果它们两两作差呢？
• 拿c来说，我们选出k个数，如果它们有两个数的差的绝对值=c，那就perfect了！

---

• 用上面的例子来说，一个班有k个同学，我们想要知道至少有两名同学的生日恰好相差c天的概率是多少。
• 而且，这次的形势比生日悖论要乐观许多。毕竟恰好相差c天，早c天、晚c天都可以；而且这里的c代指任意一个N的非平凡因子，若N不为质数，那少说也有两个c满足条件。
• 根据生日悖论的分析，当K取到$\sqrt{N}$时，概率就>50%了。

---

• 虽然很有道理，但是我们可以更稳一点。
• 在这k个数中，如果有两个数a、b，满足gcd(|a-b|,N)>1，那么我们就可以直接选这个gcd了。
• 这么一来，概率就翻了好多倍了。
• 因为N=c*d，c和d都很稀有，但是c和d的倍数便是一波大军。

流水线

• 如果我们要直接取k个数，那么要取到$N^{\frac{1}{4}}$个数，那样很耗空间。
• 可以不断生成伪随机数，然后像流水线一样判断。
• 设随机函数$f(x)=(x^2+c)modN$，其中c为一个定值，可以随机一个。

int find_factor(int N) {
    a = rand();
    for（ int i= 1; i <= 1000000; i++ ) {
        b = f(a);
        p = GCD( abs( b - a ) , N);
        if( p > 1 ) return p;//Found factor: p
        a = b;
    }
    return 0;//Failed
}

判断循环

Floyd判圈

• 发现上面的伪随机函数会循环。我们应如何判环呢？
• 可以傻逼逼地开个布尔数组，记录一下出现过的数。但是f(x)的取值范围可是0~N-1的；而当N超过$10^8$的时候，你就gg了。
• 刘汝佳先生说：想象一下，假设有两个小孩子在一个“可以无限向前跑”的跑道上赛跑，同时出发。但其中一个小孩的速度是另一个的两倍。如果跑道是直的，跑得快的小孩永远在前面；但如果跑道有环，则跑得快的小孩将“追上”跑得慢的小孩。
• 于是我们yy出了传说中的“Floyd判圈”算法。

int find_factor(int N) {
    a = rand();
    b = a;
    do {
        a = f(a);//a runs once
        b = f(f(b));//b runs twice as fast
        p = GCD( abs( b - a ) , N);
        if( p > 1 ) return p;//Found factor: p
    } while( b != a );
    return 0;//Failed
}

• 注：据维基百科说，即使n是合数，该算法也可能无法找到非平凡因子。在这种情况下，可以使用除2或不同的起始值再次尝试该方法。

另一种方法——Brent判圈

• 还有一种做法，是由Richard P. Brent提出的一种判圈法，又叫Brent判圈。
• 这种做法就是记录一下当前跑到第几个数了，如果跑到2的某个次幂，就将b=a，每次也是求GCD( abs( b - a ) , N)，若某次新随机出的a=b则return。

ll get(ll k)
{
	ll a=rand()%(k+1), b=a, c=rand()%k,r,i=1,n=2;
	while(1)
	{
		i++;
		a=f(a);
		if(a==b) return k;
		r=gcd(abs(a-b),k);
		if(1<r&&r<k) return r;
		if(i==n) b=a, n<<=1;
	}
}

• Brent声称，平均而言，他的循环搜索算法比Floyd的运行速度快36％左右，并且它将Pollard rho算法加速了大约24％。

失败了怎么破

• 如果算法失败了，我们只需重新弄一个种子，再来一发即可。
• 不过请放心，大多数时候你并不会失败。

与Miller_Rabin有机结合♂♂

• 我们可以发现pollard rho直到现在都还没有与Miller-Rabin有任何联系，但马上就不是了。
• 对于pollard rho，它可以在Θ(sqrt§)的时间复杂度内找到N的一个小因子p，这一点我们曾论证过。可见，如果N的因子很多、因子值很小的整数N来说，效率是很优异的。
• 但是，如果反过来呢？如果说N是大整数，恰好因子很少、因子值很大？

---

• 譬如，如果N是一个超大的质数，会发生什么结果？
• 我们会源源不断地生产随机数，但是每次我们都发现GCD( abs( b - a ) , N) =1，于是每次都失败；然后重设种子，继续奋战……直到最后，你都很难判断这个N是否真的不可约。
• 但是，一旦拥有Miller_Rabin，一切便都已解决。
• 我们在想要分解N的时候，预先对其进行素性判定，若其通过，则视其为素数。虽然这么做会降低我们的正确率，但却能显著提高我们的速率。

Code

ll ti(ll x,ll y,ll m)//返回x*y%m的黑科技，不会爆long long
{
	ld d=1;
	d=d*x/m*y;
	return ((x*y-((ll)d)*m)%m+m)%m;
}
ll gcd(ll x,ll y) 
{
	return y?gcd(y,x%y):x;
}
ll get(ll k)//找到k的一个非平凡因子
{
	ll a=rand()%(k+1), b=a, c=rand()%k,r,i=1,n=2;
	while(1)
	{
		i++;
		a=(ti(a,a,k)+c)%k;//a=f(a)
		if(a==b) return k;//如果a=b则遇到圈，算法失败
		r=gcd(abs(a-b),k);
		if(1<r&&r<k) return r;//找到非平凡因子
		if(i==n) b=a, n<<=1;//i达到2的某个次幂，则令b=a
	}
}
void pollard_rho(ll k)//对k进行质因数分解
{
	if(k==1) return;
	if(miller_rabin(k)) {pf[++pf[0]]=k; return;}//如果k为质数则直接return
	ll p=k;
	while(p==k) p=get(k);//不断寻找非平凡因子，直到算法成功
	pollard_rho(k/p),pollard_rho(p);//递归分解k/p和p
}

例题

【JZOJ4732】【NOIP2016提高A组模拟8.23】函数（线筛求欧拉函数+Pollard_Rho大数分解）