问题引入

给定一个正整数 aZ+a \in \mathbb{Z}^+,请找到其一个因子。

朴素算法

要找到大整数的一个因子,最朴素的想法是试除法。也就是这样:

int Find_Factor(int x) {
    for(int i = 2; i * i <= N; i++)
        if(n % i == 0)
        	return i;
    return n;
}

这样的时间复杂度是 Θ(n)\Theta(\sqrt n) 的,不能够处理 101810^{18} 及以上的大数。

基于随机的算法

我们可以先用 Miller-Rabbin 算法检测质数,如果这个数不是质数,每次从 [2n1][2,n - 1] 中随机一个数字进行判断,也就是这样:

int RandInt(int l , int r) {
	static mt19937 Rand(time(0));
	uniform_int_distribution<int> dis(l, r);
	return dis(Rand);
}

int Find_Factor(int n) {
	if(Miller_Rabbin(n))
		return n;
	int x;
	do {
		x = RandInt(2 , n - 1);
	} while(n % x);
	return x;
}

上述算法最差时间复杂度为 O(n)O(n)

这个糟糕的随机算法可以进行改进:

int RandInt(int l , int r) {
	static mt19937 Rand(time(0));
	uniform_int_distribution<int> dis(l, r);
	return dis(Rand);
}

int Find_Factor(int n) {
	if(Miller_Rabbin(n))
		return n;
	int x , y;
	do {
		x = RandInt(2 , n - 1);
		y = __gcd(n , x);
	} while(y == 1);
	return y;
}

分析最差时间复杂度:当 nnp2p^2pp 为质数),显然当随机数 xxp,2p,3p,,(p1)pp , 2p , 3p , \cdots , (p-1)p 时都可以找到 gcd(x,n)=p>1\gcd(x , n) = p > 1,所以时间复杂度为 O(nlogn)O(\sqrt n \log n) 的,连暴力都打不过,但这是 Pollard-Rho 算法的基础。

生日悖论

这个悖论是 Pollard-Rho 算法的前置知识。

在一个房间中至少有多少人可以使得同一天出生的人出生的概率为 50%50\%

假设一年有 nn 天,房间里共有 kk 个人,编号为 1k1 \sim k,假设生日随机分布。

kk 个人生日不相同为事件 AA,则有:

Pr{A}=i=nk+1nin\operatorname{Pr}\{A\} = \sum_{i=n-k+ 1}^n \frac{i}{n}

同一天出生的人出生的概率为 Pr{A}=1Pr{A}\operatorname{Pr}\{\overline A\} = 1-\operatorname{Pr}\{A\},当 Pr{A}12\operatorname{Pr}\{ \overline A\} \ge \frac{1}{2} 时显然有 Pr{A}12\operatorname{Pr}\{A\} \le \frac{1}{2}

使用 Mathemaica 计算出来 n=365n = 365k=23k = 23

实际上可以使用随机变量指示器进行近似分析:

设第 ii 个人和第 jj 个人同时在第 ww 天出生为事件 BB,有

Pr{N}=1n2\operatorname{Pr}\{N\}=\frac{1}{n^2}

由于一年有 nn 天,所以有

Pr{i,j 在同一天出生}=1n\operatorname{Pr}\{i,j \ \text{在同一天出生}\} = \frac{1}{n}

设随机变量指示器 Xi,jX_{i,j},定义:

Xi,j={1当 i 和 j 在同一天出生0OtherwiseX_{i,j} =\left\{ \begin{array}{lc} 1 &\text{当} \ i \ \text{和} \ j\ \text{在同一天出生}\\ 0 &Otherwise \end{array} \right.

那么有:

E[Xi,j]=Pr{i,j 在同一天出生}=1n\mathbb{E}[X_{i,j}] = \operatorname{Pr}\{i,j \ \text{在同一天出生}\} = \frac{1}{n}

XX 为两个人在同一天出生的人的无序对数,即 X=i=1kj=i+1kXi,jX = \sum\limits_{i=1}^k \sum\limits_{j = i + 1}^kX_{i,j} 则有:

E[X]=E[i=1kj=i+1kXi,j]\mathbb{E}[X] = \mathbb{E}[\sum\limits_{i=1}^k \sum\limits_{j = i + 1}^kX_{i,j}]

由数学期望的线性性质得:

E[i=1kj=i+1kXi,j]=i=1kj=i+1kE[Xi,j]=i=1kj=i+1k1n=k(k1)2×1n \mathbb{E}[\sum\limits_{i=1}^k \sum\limits_{j = i + 1}^kX_{i,j}] = \sum\limits_{i=1}^k \sum\limits_{j = i + 1}^k \mathbb{E}[X_{i,j}] = \sum\limits_{i=1}^k \sum\limits_{j = i + 1}^k \frac{1}{n} = \frac{k(k-1)}{2} \times \frac{1}{n}

解得当 n=365n = 365 时,房间中有 2828 人就可以期望有两人生日相同。

于是可以这样优化:

原来从 [1,n][1,n] 中随机出我们想要的数字 kk 的概率为 1n\frac{1}{n},现在我们随机两个数 i,ji,j 使得 i,j=k|i,j|=k 的概率大约为 2n\frac{2}{n},概率扩大了两倍左右。

伪随机数序列

Pollard-Rho 使用一种特别的伪随机数生成器来生成 [1,n1][1,n-1] 间的伪随机数序列:设序列第一个数为 xxf(x):=(x2+c)modnf(x) := (x^2+c) \mod n,其中 cc 为可以自行指定的常数,则 x,f(x),f(f(x)),f(f(f(x)))x,f(x),f(f(x)),f(f(f(x))) 为一个伪随机数序列。

显然的,每个数都由前一个数字决定,以生成的数又是有限的,那么迟早会进入循环,而且大概率为混循环,所以生成的序列类似一个 ρ\rho 形,故名 ρ\rho 算法。

遇到环后我们接下来做的随机都没有意义了,所以需要更换 cc 重新开始。

Floyd 判环算法

设置两个变量 t,rt,r,每次判断是否有 gcd(tr,n)>1\gcd(|t-r|,n) > 1,如果没有,就令 t=f(t)t=f(t)r=f(f(r))r=f(f(r))。因为 rr 跑的更快,那么它们最后会相遇,这时候就换一个 cc 重新生成随机数。

这个伪随机数生成器有一个性质:如果 ij0(mod p)|i-j| \equiv 0 (\operatorname{mod} \ p),那么有 f(i)f(j)=i2j2=ij×i+j0(mod p)|f(i)-f(j)| = |i^2 - j^2| = |i - j| \times|i + j| \equiv 0 (\operatorname{mod} \ p)

由此可得,只要环上距离为 xx 的两个数满足条件,那么所有距离为 xx 的数都满足条件。在Floyd判环的过程中,每次移动都相当于在检查一个新的距离 xx ,这样就不需要进行两两比较了。

这个算法的复杂度依赖于这个伪随机数生成器的随机程度,还没有被严格证明。如果它是足够随机的,那么期望复杂度显然是 Θ(n14logn)\Theta(n^{\frac{1}{4}}\log n)

ll RandInt(ll l , ll r) {
	static mt19937 Rand(time(0));
	uniform_int_distribution<ll> dis(l, r);
	return dis(Rand);
}

ll Pollard_Rho(ll n) {
	if(n == 4) {
		return 2; 
	} 
	if(Miller_Rabin(n)) {
		return n;
	}
	while(true) {
		ll c = RandInt(1 , n - 1);
		auto f = [=](ll x) {
			return ((__int128)x * x + c) % n;
		};
		ll t = f(0) , r = f(f(0));
		while(t != r) {
			ll d = __gcd(abs(t - r) , n);
			if(d > 1)
				return d;
			t = f(t);
			r = f(f(r));
		}
	}
}

倍增优化

由于 gcd(a,b)\gcd(a,b) 的运算会花费 O(logn)O(\log n) 的时间,从而使我们的时间复杂度挂上了 log\log,这使人非常不爽,现在我们想办法去掉这个 log\log

我们使用 乘法累计 来减少求 gcd\gcd 的次数。如果 gcd(a,b)>1\gcd(a,b) > 1,那么有 gcd(ac,b)>1,cN\gcd(ac,b)>1,c\in \mathbb{N}^*,并且有 gcd(ac mod b,b)=gcd(a,b)\gcd(ac\ \operatorname{mod} \ b,b) = \gcd(a,b)

我们每过一段时间对这些差值进行 gcd\gcd 运算,设 s=x0x1 mod ns = \prod|x_0-x_1| \ \operatorname{mod} \ n,如果过程中存在 s=0s = 0 测表示分解失败, nn 是质数。每隔 2k12^k-1 个数,计算是否满足 1<gcd(s,n)<n1 < \gcd(s,n) < n。此处 k=7k=7,可以更具需要自行调节。

typedef long long ll;
typedef __int128 lll;

ll RandInt(ll l , ll r) {
	static mt19937 Rand(time(0));
	uniform_int_distribution<ll> dis(l, r);
	return dis(Rand);
}

ll Pollard_Rho(ll n) {
	if(Miller_Rabin(n)) {
		return n;
	}
	ll s = 0 , t = 0;
	ll c = RandInt(1 , n - 1);
	int step = 0 , goal = 1;
	ll value = 1;
	auto f = [=](ll x) {
		return ((lll)x * x + c) % n;
	};
	for(goal = 1;; goal <<= 1, s = t , value = 1) {
		for(step = 1; step <= goal; step++) {
			t = f(t);
			value = ((lll)value * abs(t - s)) % n;
			if(step % 127 == 0) {
				ll d = __gcd(value , n);
				if(d > 1) return d;
			}
		}
		ll d = __gcd(value , n);
		if(d > 1) return d;
	}
}

【模板】Pollard-Rho算法

题目传送门

思路

对于一个数 nn,如果用 Miller_Rabin 判断出来他是质数,就可以直接返回。

否则找到它的一个真因子 pp,把 nn 除去所有的 pp 因子,再递归分解 nnpp,其中使用 Miller_Rabin 判断是否存在质因子,并打擂台取最大值。

代码

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef __int128 lll;

ll RandInt(ll l , ll r) {
	static mt19937 Rand(time(0));
	uniform_int_distribution<ll> dis(l, r);
	return dis(Rand);
}

ll Quick_Power(ll a , ll b , ll p) {
	ll res = 1;
	while(b) {
		if(b & 1) res = (lll)res * a % p;
		a = (lll)a * a % p;
		b >>= 1;
	}
	return res;
}

bool Miller_Rabin(ll n) {
	if(n < 3 || n % 2 == 0) return n == 2;
	ll a = n - 1 , b = 0;
	while(a % 2 == 0) {
		a /= 2;
		b++;
	}
	for(int i = 1 , j; i <= 10; i++) {
		ll x = RandInt(2 , n - 1);
		ll v = Quick_Power(x , a , n);
		if(v == 1) continue;
		for(j = 0; j < b; j++) {
			if(v == n - 1) break;
			v = (lll) v * v % n;
		}
		if(j == b) return 0;
	}
	return 1; 
}

ll Pollard_Rho(ll n) {
	ll s = 0 , t = 0;
	ll c = RandInt(1 , n - 1);
	int step = 0 , goal = 1;
	ll value = 1;
	auto f = [=](ll x) {
		return ((lll)x * x + c) % n;
	};
	for(goal = 1;; goal <<= 1, s = t , value = 1) {
		for(step = 1; step <= goal; step++) {
			t = f(t);
			value = ((lll)value * abs(t - s)) % n;
			if(step % 127 == 0) {
				ll d = __gcd(value , n);
				if(d > 1) return d;
			}
		}
		ll d = __gcd(value , n);
		if(d > 1) return d;
	}
}

ll Ans;

void Fac(ll n) {
	if(n <= Ans || n < 2) return;
	if(Miller_Rabin(n)) {
		Ans = max(Ans , n);
		return;
	}
	ll p = n;
	while(p == n) p = Pollard_Rho(n);
	while((n % p) == 0) n /= p;
	Fac(n);
	Fac(p); 
}

ll T , N;

int main() {
	ios::sync_with_stdio(false);
	cin.tie(0);
	cout.tie(0);
	cin >> T;
	while(T--) {
		Ans = 0;
		cin >> N;
		Fac(N);
		if(Ans == N) {
			cout << "Prime" << endl;
		}
		else {
			cout << Ans << endl;
		}
	}
	return 0;
}