Pollard-Rho

算法简介

Pollard-Rho 主要用于质因数分解，是一个随机化算法。

前置算法

Miller-Rabin

伪随机函数

此算法用 (f_x=f_{x-1}^2+c mod n) 生成随机数用于寻找质因数。
（其中，(f_1) 和 (c) 可以随机生成， (n) 为正在分解的数）
但是，这个函数会产生环，所以叫伪随机函数。

判环

考虑用两个指针，一个每次走一步，另一个每次走两步。
如果两个指针的值在某一时刻相等，则找到了环。
因为两个指针位置的差每次增加1，所以不用担心找不到环。

算法主体

如果当前的数是质数，记录一下，直接返回。
判质数可用 Miller-Rabin 实现。
用上面的伪随机函数生成 (f) 数组。
根据生日悖论，随机两个数 (x) ， (y) ，使得 (gcd(|x-y|,n)>1) 的概率比只随机一个数的概率要大。
所以，用判环使用的两个指针的值作为 (x) ， (y) ，
设 (d=gcd(|x-y|,n)) 。
若 (d>1) 且 (d<n)，则递归处理 (d) 和 (frac{n}{d}) 两个部分。
同时，如果找到了环，且还没有找到合法的因数，就重新随机 (f_1) 和 (c) 。

现在的代码

#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<ctime>
#define int long long
#define LD long double
#define ull unsigned long long
using namespace std;
int T,n,ans;
int MUL(int a,int b,int p) //a*b%p
{
	int x=(LD)a/p*b;
	return ((ull)a*b-(ull)x*p+p)%p;
}
int POW(int a,int b,int p) //a^b%p
{
	if(!b) return 1;
	if(b==1) return a;
	int sum=POW(a,b/2,p);
	if(b%2) return MUL(MUL(sum,sum,p),a,p);
	return MUL(sum,sum,p);
}
int MAX(int x,int y)
{
	return x>y?x:y;
}
int ABS(int x)
{
	return x>0?x:-x;
}
int gcd(int x,int y)
{
	if(!y) return x;
	return gcd(y,x%y);
}
int f(int x,int c,int p)
{
	return (MUL(x,x,p)+c)%p;
}
bool MR(int x)
{
	if(x==0||x==1) return false;
	if(x==2) return true;
	if(x%2==0) return false;
	int p=x-1,q=0;
	while(p%2==0) q++,p/=2;
	for(int i=1;i<=10;i++)
	{
		int a=rand()%(x-1)+1;
		if(POW(a,x-1,x)!=1) return false;
		int lst=1;
		for(int j=0;j<q;j++)
		{
			int t=POW(a,(1ll<<j)*p,x);
			if(t==1&&lst!=1&&lst!=x-1) return false;
			lst=t;
		}
		if(lst!=1&&lst!=x-1) return false;
	}
	return true;
}
int find(int x)
{
	if(x%2==0) return 2;
	if(MR(x)) return x;
	int t=rand()%(x+1);
	int a=t,b=t;
	int c=rand()%(x+1);
	while(1)
	{
		a=f(a,c,x),b=f(f(b,c,x),c,x);
		int d=gcd(ABS(a-b),x);
		if(d>1&&d<x) return d;
		if(a==b) return find(x);
	}
}
void PR(int x)
{
	if(x<=1) return;
	if(MR(x))
	{
		ans=MAX(ans,x);
		return;
	}
	int y=find(x);
	PR(y),PR(x/y);
}
signed main()
{
	srand(time(0));
	scanf("%lld",&T);
	while(T--)
	{
		scanf("%lld",&n);
		if(MR(n)) puts("Prime");
		else
		{
			ans=0,PR(n);
			printf("%lld
",ans);
		}
	}
	return 0;
}

黑科技 gcd

上面的代码因为效率过慢过不了洛谷上的 模板题。
观察到代码中会使用很多次 gcd 函数，所以考虑优化计算 gcd 的方法。

新算法
设现在要计算 (a) 和 (b) 两数的 gcd。
若 (a) 和 (b) 都是偶数，则递归计算 (2 imes gcd(frac{a}{2},frac{b}{2})) 。
若 (a) 是偶数，且 (b) 是奇数，则递归计算 (gcd(frac{a}{2},b)) 。
若 (a) 是奇数，(b) 是偶数同理。
若 (a) 和 (b) 都是奇数，则递归计算 (gcd(|a-b|,a)) 。
证明：
设原来 (a) ， (b) 的 gcd 为 (k) 。
则 (a=kA) ， (b=kB) 。
(|a-b|=k|A-B|) 。
因为 (A) 与 (B) 必然互质，
所以，现在要证明的就是 (A) 与 (|A-B|) 互质。
为了方便，先假设 (A>B) ，去掉绝对值。相反的情况的证明也是类似的。
考虑反证法。
设 (A) 与 (A-B) 不互质，
则 (A=qx) ， (A-B=qy) 。
则 (B=A-(A-B)=qx-qy=q(x-y)) ，所以 (A) 与 (B) 不互质，不成立。
所以前面的递归是成立的。

把原算法中的 gcd 部分换成这种写法，并加上一些比较常见的卡常技巧即可轻松通过本题。
贴一个用迭代实现的板子。

//__builtin_ctzll(x) 函数用于求x在二进制中末尾0的个数
int gcd(int x,int y)
{
	if(!x) return y;
	if(!y) return x;
	int t=__builtin_ctzll(x|y);
	x>>=__builtin_ctzll(x);
	do
	{
		y>>=__builtin_ctzll(y);
		if(x>y) swap(x,y);
		y-=x;
	}while(y);
	return x<<t;
}