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=57,n=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-k或num_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)。通过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;
}
}
对于一个数 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;
}