6.1.5 质因数分解
质因数分解
根据唯一分解定理
任何一个大于 1 的整数都可以被分解成有限个质数的乘积的形式
\[n = \prod_{i = 1}^m p_i^{c_i}=p_1^{c_1} \times p_2^{c_2} \times \cdots \times p_m^{c_m} \]
我们有两种算法来分解一个数
试除法
直接枚举 \(2 \sim \sqrt n\) 的数，如果能整除，就一直除，直到把这个因子除完，时间复杂度 \(O(\sqrt n)\)
int p[N], c[N], cnt;
void divide (int n) {
cnt = 0;
for (int i = 2; i * i <= n; i++) {
if (n % i == 0) {
p[++cnt] = i, c[cnt] = 0;
while (n % i == 0) n /= i, c[cnt]++;
}
}
if (n > 1) p[++cnt] = n, c[cnt]++;
}
\(Pollard-Rho\) 算法
这个神人算法让我们能够在 \(O(\sqrt p)\) 也就是 \(O(\sqrt[4] n)\) 的时间复杂度内求解出 \(n\) 的所有质因子，它和 \(Miller-Rabin\) 一样，都是随机算法，好消息是准确率同样高得离谱
我们先理解一下大概的思路，我们现在的目的是将 \(n\) 分解质因数
现在的问题是，我们没办法快速的找到 \(n\) 的某个质因数 \(p\)
所以 \(Pollard\) 发明了一种算法，利用数列（在 \(\bmod p\) 意义下）循环的性质，找到 \(n\) 的某个因数，然后递归求解，即可找到 \(N\) 的所有质因数
生日悖论
在讲具体算法之前，我们先说说这个反直觉的生日悖论。其实生日悖论不是一个逻辑悖论，只是很反直觉罢了
一个房间里有 23 个人，他们中两人生日相同的概率超过一半（不考虑闰年，他们的生日完全随机）
证明：正难则反，生日相同的概率不好统计，但生日不同的概率较好统计，设 \(A\) 为这 23 个人生日都互不相同
如何理解？我们考虑第前两个人生日不在同一天的概率，自然就是除去第一个人的生日，第二个人在剩下的日子里随便选，后面的也一样
生日悖论能告诉我们，如果我们不断在某个范围内生成随机整数，那么很快就会出现重复的数，期望大概是根号级别的。更准确的说，对于一个 \([1,N]\) 范围内整数的理想随机数生成器，生成序列中第一个重复数字前，期望有 \(\sqrt {\frac {\pi N} 2}\) 个数
应用到题目上，即对于最坏情况 \(n = p^2\) 如果我们不断在 \(1 \sim n - 1\) 之间生成随机数，那么期望在生成 \(\sqrt p = \sqrt[4] n\) 个数后出现在 \(\bmod p\) 意义下相同的数，那么这两个数的差的绝对值 \(d\) 就满足 \(d \equiv 0 \ (\bmod p)\) ，于是也满足 \(\gcd(d, n) > 1\)
但这个东西并不能直接用，因为它是在 两两比较 的情况下才成立的，如果我们枚举两个数进行比较，那么时间复杂度退化到 \(O(\sqrt n \log n)\) ，甚至不如试除法
伪随机数列
\(Pollard\) 代佬想出了一种很神奇的伪随机数生成器来生成 \([0, N - 1]\) 间的伪随机数序列，设序列第一个数为 \(x\) ，\(f(x) = (x^2 + c) \bmod N\) ，则 \(x,f(x),f(f(x)),\cdots\) 为一个伪随机数序列
这个序列有个简单性质，就是后一项的值是根据前一项的值推导过来的，那么一旦出现重复的数，就会进入循环，如下图（\(c = 1\)）
这里我们假定这个伪随机数列足够随机，那么期望时间复杂度就是 \(O(\sqrt p)\) 的
那么这个序列怎么就能够降低时间复杂度了呢？
它还有个非常重要的性质
如果 \(|i - j| \equiv 0 \ (\bmod p)\) ，那么一定有 \(|f(i) - f(j)| = |i^2 - j^2| = |i - j| \cdot |i + j| \equiv 0 \ (\bmod p)\)
由此可得，只要距离为 \(d\) 的两个数满足条件，那么所有距离为 \(d\) 的数均满足条件，那么我们每算一个 \(d\) 其实相当于算了这个序列上所有这个距离的数的差值，所以我们每个距离只需要算一次即可，下一次就去算下一个距离
如果伪随机数序列足够随机，那么期望的时间复杂度为\(O(\sqrt[4]n \log n)\)
inline ll Pollard_rho (ll n) {
static mt19937_64 eg ((unsigned int)time (0));
uniform_int_distribution <ll> u0 (1, n - 1);
ll c = u0 (eg);
auto f = [&](ll t) { return ((i8)t * t + c) % n; };
ll x = f(0), y = f(x);
while (x != y) {
ll d = gcd (abs (x - y), n);
if (d != 1) return d;
x = f(x), y = f(f(y));
}
return n;
}
但实际上还可以进一步优化，如果 \(\gcd(d, n) > 1\) 那么 \(\gcd (kd \bmod n, n) > 1\) ，所以我们可以减少 \(\gcd\) 的次数，即先将待选数乘起来，再统一求 \(\gcd\)
我们可以倍增去求，每隔 \(1,2,4,8,\cdots\) 求一次 \(\gcd\) ，真阳只需要求期望 \(\log(\sqrt[4] n)\) 次公因数，总期望时间复杂度为 \(O(\sqrt[4] n + \log(\sqrt[4] n)\log n)\) ，可以看做 \(O(\sqrt[4] n)\)
但这样有个缺点，后面间隔很大，可能已经满足条件却迟迟无法退出，所以我们固定一个间隔，每隔这么多次就求一次 \(\gcd\)
code
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <vector>
#include <ctime>
#include <random>
using namespace std;
namespace michaele {
typedef long long ll;
typedef __int128 i8;
int pri[] = {2, 3, 5, 7, 13, 29, 37};
vector <ll> fac;
inline ll pw (ll a, ll b, ll mod) {
ll res = 1, t = a;
while (b) {
if (b & 1) {
res = (i8)res * t % mod;
}
t = (i8)t * t % mod;
b >>= 1;
}
return res;
}
inline ll gcd (ll a, ll b) {
if (a == 0 || b == 0) return a | b;
if (a == b) return a;
int az = __builtin_ctzll (a);
int bz = __builtin_ctzll (b);
int z = min (az, bz);
ll diff;
b >>= bz;
while (a) {
a >>= az;
diff = a - b;
if (!diff) break;
az = __builtin_ctzll (diff);
b = a < b ? a : b;
a = diff < 0 ? -diff : diff;
}
return b << z;
}
inline bool MR (ll n) {
if (n < 3) return n == 2;
if (!(n & 1)) return false;
ll d = n - 1, k = 0;
while (!(d & 1)) d >>= 1, k++;
for (int i = 0; i < 7; i++) {
ll v = pw (pri[i], d, n);
if (v <= 1 || v == n - 1) continue;
for (int j = 1; j <= k; j++) {
v = (i8)v * v % n;
if (v == n - 1 && j != k) {
v = 1;
break;
}
if (v == 1) return false;
}
if (v != 1) return false;
}
return true;
}
inline ll Pollard_rho (ll n) {
static mt19937_64 eg ((unsigned int)time (0));
uniform_int_distribution <ll> u0 (1, n - 1);
ll c = u0 (eg);
auto f = [&](ll t) { return ((i8)t * t + c) % n; };
ll x = 0, y = 0, s = 1, k = 1;
while (1) {
for (int i = 1; i <= k; i++) {
x = f(x);
s = (i8)s * abs (x - y) % n;
if ((i & 127) == 0) {
ll d = gcd (s, n);
if (d > 1) return d;
}
}
ll d = gcd (s, n);
if (d > 1) return d;
y = x;
s = 1;
k <<= 1;
}
while (x != y) {
ll d = gcd (abs (x - y), n);
if (d != 1) return d;
x = f(x), y = f(f(y));
}
return n;
}
ll mafac;
void get_fac (ll n) {
if (n == 1) return;
if (n == 4) {
fac.push_back (2);
fac.push_back (2);
return;
}
if (MR (n)) {
mafac = max (mafac, n);
return;
}
ll tmp = n;
while (tmp == n) tmp = Pollard_rho (n);
get_fac (tmp), get_fac (n / tmp);
}
inline void rd (ll &res) {
res = 0; char c = getchar ();
while (c < '0' || c > '9') c = getchar ();
while (c >= '0' && c <= '9') {
res = (res << 1) + (res << 3) + (c ^ 48);
c = getchar ();
}
}
void main () {
ios :: sync_with_stdio(0);
cin.tie(0);
cout.tie (0);
ll T, x;
rd(T);
while (T--) {
fac.clear ();
rd(x);
if (MR (x)) {
cout << "Prime" << "\n";
continue;
}
mafac = 0;
get_fac (x);
cout << mafac << "\n";
}
}
}
int main () {
michaele :: main ();
return 0;
}
浙公网安备 33010602011771号