发表于 | 分类于 | 0 | 阅读次数 945

问题引入

给定一个正整数 a∈Z+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 算法检测质数，如果这个数不是质数，每次从 [2，n−1][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;
}

分析最差时间复杂度：当 nn 为 p2p^2（pp 为质数），显然当随机数 xx 取 p,2p,3p,⋯ ,(p−1)pp , 2p , 3p , \cdots , (p-1)p 时都可以找到 gcd⁡(x,n)=p>1\gcd(x , n) = p > 1，所以时间复杂度为 O(nlog⁡n)O(\sqrt n \log n) 的，连暴力都打不过，但这是 Pollard-Rho 算法的基础。

生日悖论

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

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

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

设 kk 个人生日不相同为事件 AA，则有：

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

同一天出生的人出生的概率为 Pr⁡{A‾}=1−Pr⁡{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 = 365 时 k=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=1k∑j=i+1kXi,jX = \sum\limits_{i=1}^k \sum\limits_{j = i + 1}^kX_{i,j} 则有：

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

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

E[∑i=1k∑j=i+1kXi,j]=∑i=1k∑j=i+1kE[Xi,j]=∑i=1k∑j=i+1k1n=k(k−1)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,n−1][1,n-1] 间的伪随机数序列：设序列第一个数为 xx，f(x):=(x2+c)mod  nf(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⁡(∣t−r∣,n)>1\gcd(|t-r|,n) > 1，如果没有，就令 t=f(t)t=f(t)，r=f(f(r))r=f(f(r))。因为 rr 跑的更快，那么它们最后会相遇，这时候就换一个 cc 重新生成随机数。

这个伪随机数生成器有一个性质：如果 ∣i−j∣≡0(mod⁡ p)|i-j| \equiv 0 (\operatorname{mod} \ p)，那么有 ∣f(i)−f(j)∣=∣i2−j2∣=∣i−j∣×∣i+j∣≡0(mod⁡ p)|f(i)-f(j)| = |i^2 - j^2| = |i - j| \times|i + j| \equiv 0 (\operatorname{mod} \ p)

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

这个算法的复杂度依赖于这个伪随机数生成器的随机程度，还没有被严格证明。如果它是足够随机的，那么期望复杂度显然是 Θ(n14log⁡n)\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(log⁡n)O(\log n) 的时间，从而使我们的时间复杂度挂上了 log⁡\log，这使人非常不爽，现在我们想办法去掉这个 log⁡\log。

我们使用 乘法累计 来减少求 gcd⁡\gcd 的次数。如果 gcd⁡(a,b)>1\gcd(a,b) > 1，那么有 gcd⁡(ac,b)>1,c∈N∗\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=∏∣x0−x1∣ mod⁡ ns = \prod|x_0-x_1| \ \operatorname{mod} \ n，如果过程中存在 s=0s = 0 测表示分解失败， nn 是质数。每隔 2k−12^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 因子，再递归分解 nn 和 pp，其中使用 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;
}
  • 本文作者： CmsMartin
  • 本文链接： /archives/pollard-rho
  • 版权声明： 本博客所有文章除特别声明外，均采用CC BY-NC-SA 3.0 许可协议。转载请注明出处！