[科技]\(Miller\_Rabin\) 和 \(Pollard\_Rho\) 及各种玄学优化
[科技] \(Miller\_Rabin\) 和 \(Pollard\_Rho\)
先讲\(Miller\_Rabin\)吧,\(Miller\_Rabin\)是用来检验素数的高效算法。
我们先要知道两个定理:
- 费马小定理:当\(p\)为质数时,\(x^{p - 1} \equiv 1 \ \ (mod \ \ p)\)。但这只是一个充分条件,但不是必要条件。即就算\(x\)和\(p\)互质,那么\(x^{p - 1}\)不一定在模\(p\)意义下同余于\(1\)。
- \(Fermat\)定理:若\(p\)为奇素数,且\(0 < x < p\),那么\(x ^ 2 \equiv 1 \ \ (mod \ \ p)\)的解为\(x = 1\)或\(x = p - 1\)。这个还是比较好证明的,移项可得\(x ^ 2 - 1 \equiv 0\),即\((x - 1)(x + 1) \equiv 0\),可得\(p | (x - 1)(x + 1)\)。而\(p\)是质数,那么可知\(x = 1\)或\(x = p - 1\)。
首先我们根据费马小定理就可以排除大量的合数了,即如果存在\(x\),不满足\(x ^ {p - 1} \equiv 1\),则\(p\)就不是质数了。但是之前也说过,这个并不是一个必要条件,所以也存在满足费马小定理的合数,这种数叫做\(Carmichael\)数。最小的\(Carmichael\)数是\(561\),对于任意的数,都有\(x ^ {560} \equiv 1\),而\(561 = 3 * 11 * 17\)。\(Carmichael\)数使得我们无法使用费马小定理来判断一个数是否为素数。
然后\(Miller\_Rabin\)算法就是优化了这一点,尝试使用多次判定来提高正确性。当我们判断了一个数满足了费马小定理之后,我们继续判断\(x ^ { (p - 1) / 2} \equiv 1\)是否正确。于是我们就得出了一个算法:我们将\(p - 1\)表示成\(2 ^ s * t\)的形式,然后选取几个比较小的质数\(pri\),从\(pri ^ t\)开始,依次判断\(pri ^ t, pri ^ {2t},pri ^ {4t} \dots\)等等是否满足这个同余式。如果对于所有取的素数都满足时,就可以确定这个数大概率是个素数了。经过测试发现,对于\(int\)范围以内的数,取出小于\(30\)的素数来验证就可以保证完全正确的判断出是否为素数,在\(long\ \ long\)范围内的错误的概率也是可以忽略不计的。
再来讲讲\(Phollard\_Rho\),他是基于\(Miller\_Rabin\)的一种分解质因数的方法,同时也是基于随机化的方法,期望复杂度为\(O(n ^ {\frac{1}{4}})\),但实际情况下的复杂度也是个玄学,似乎表现的非常优秀。
假设我们需要对于\(x\)分解质因数,我们先通过\(Miller\_Rabin\)素性测试判断给定的\(x\)是否为素数。然后开始\(Pollard\_Rho\)进行分解。
我们先要了解的是\(Pollard\_Rho\)是将某合数\(x\)分解为两个非平凡因子\(a, b\)(\(1\) 和 \(x\)为平凡因子)的算法。如果需要对\(x\)分解质因数,那么递归下去做就行了。
先思考一下平时的\(O(\sqrt n)\)分解质因数的方法,即枚举\(2\)到\(\sqrt n\)的质因数,判断是否能被整除。实际上一个数的质因数是\(log\)级别的,所以我们其中的枚举很多都是无用枚举。考虑如何来优化这个枚举方法。
有一个悖论叫做生日悖论,意思是在一个人数为\(23\)人的班级中,存在两个人生日相同的概率接近\(0.5\),而在人数大约为\(60\)人的班级中,存在两个人生日相同的概率已经非常接近\(1\)了。这虽然与实际非常不相符,但是简单计算一下发现的确是这样。以\(23\)人为例,我们计算没有两个人生日相同的概率\(P = \frac{365}{365} * \frac{364}{365} * \frac{363}{365} * \dots * \frac{365 - 22}{365}\),此时的\(P\)大约为\(0.49\)。
举一个简单的例子就是假设我们需要在\([1, 100]\)之间选一个数,那么选到\(39\)的概率为\(1 \%\),但是如果我们选择两个数\(a, b\),使得\(| a - b | = 39\),那么概率就变成了\(2 \%\)。
这给了我们启发,我们在选择测试因子时,也采用这种组合随机采样方法。我们随机两个数\(a, b\),判断\(gcd( | a - b |, n)\)是否大于1,如果大于\(1\),则我们找到了一个因子,这样会大大提高我们找到因子的概率。
那么我们该如何生成这些随机数呢?\(Pollard\_Rho\)通过伪随机数的方法来提供待测因子,即选定一个起始数字\(x_0\),以及一个常数\(c\),通过\(x_i = x_{i - 1} ^ 2 + c\ \ (mod \ \ n)\)的方法来生成一系列的随机因子。至于为什么选用这个函数呢?似乎是这个函数在带入复数之后迭代得出的点集叫曼德勃罗集,这个集合又与混沌系统有关,即这个集合中的一系列数字非周期又不收敛,这使得这个函数生成的随机数非常优秀。并且这些生成的随机数,每一个都完全依赖于前一个数字,那么在迭代了一定次数之后,一定会出现循环。所以当我们出现循环的时候,就重新设置\(x_0\)与\(c\),继续检测,这是不用系统的\(rand\)函数的一个原因,因为系统的随机并不一定会出现循环。而这个循环也构成了一个像\(\rho\)一样的形状,这也是其名称后缀\(Rho\)的来源。
然后我们就可以设计一个基于倍增的随机算法了,每次在\(\rho\)的路径上取\([2 ^ {k - 1}, 2 ^ k]\)这段区间,然后对于所有的\(i \in [2 ^ {k - 1}, 2 ^ k]\),判断\(gcd(x_i - x_{2^{k - 1}}, n)\)是否大于1,如果大于1,则找到了一个因子,返回即可。这样既可以不在\(\rho\)的环上停留太久,也可以减少\(gcd\)的次数。
于是我们得出了\(Pollard\_Rho\)的大致算法:
- 对于需要检验的数\(n\),用\(Miller\_Rabin\)进行素性测试,如果是素数,则直接返回即可,否则进入第二步。
- 随机生成数\(x_0, c\),然后开始倍增的判断\(x_i - x_l\)与\(n\)是否有公因数,如果找到了公因数,返回即可。若已经出现了循环,则返回\(n\),表示该次探查失败,并重复2步骤。
这样我们就可以非常快速的写出\(Pollard\_Rho\)的算法了,这里放出例题
Code:
#pragma GCC optimize (2, "inline", "Ofast")
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int pri[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
int T;
ll n, ans = 0;
ll QMul(ll x, ll y, ll Md) {
ll ans = 0;
for(; y; y >>= 1, x = (x + x) % Md) if(y & 1) ans = (ans + x) % Md;
return ans;
}
ll Qpow(ll x, ll y, ll Md) {
ll ans = 1;
for(; y; y >>= 1, x = QMul(x, x, Md)) if(y & 1) ans = QMul(ans, x, Md);
return ans;
}
bool Miller_Rabin(ll x) {
if(x == 2) return 1;
if(!(x & 1) || x == 1) return 0;
ll s = 0, t = x - 1;
while(!(t & 1)) {
t >>= 1;
s ++;
}
for(int i = 0; i < 10 && pri[i] < x; i++) {
ll a = pri[i], b = Qpow(a, t, x);
ll k;
for(int j = 1; j <= s; j++) {
k = QMul(b, b, x);
if(k == 1 && b != 1 && b != x - 1) return 0;
b = k;
}
if(b != 1) return 0;
}
return 1;
}
ll Gcd(ll a, ll b) {
return !b ? a : Gcd(b, a % b);
}
ll Pollard_Rho(ll n, ll c) {
ll x = 1ll * rand() * rand() % (n - 2) + 1, y = x;
int i = 1, k = 2;
while(1) {
i++;
x = QMul(x, x, n);
x = (x + c) % n;
ll G = Gcd(y - x, n);
if(G > 1 && G < n) return G;
if(y == x) return n;
if(i == k) {
y = x;
k <<= 1;
}
}
}
void Find(ll x, ll c) {
if(x == 1 || x < ans) return ;
if(Miller_Rabin(x)) return (void) (ans = max(ans, x));
ll p = x;
while(p == x) p = Pollard_Rho(p, c--);
while(x % p == 0) x /= p;
Find(p, c); Find(x, c);
}
int main() {
srand(2333);
scanf("%d", &T);
while(T--) {
scanf("%lld", &n);
ans = 0;
bool f = Miller_Rabin(n);
if(f) {
puts("Prime");
continue;
}
Find(n, 1000000);
printf("%lld\n", ans);
}
return 0;
}
好的……本科技讲解到此结束
\(nmdwsm!!!\)
为什么交到\(luogu\)的模板题上\(T\)飞了!!本机\(0.7s\) \(luogu\)上居然\(T\)了,难以置信……
于是我们开始漫漫的卡常优化之旅……
首先注意到我们在乘法的时候为了防止\(long\ \ long\)相乘导致溢出的情况,我们使用了龟速乘……实际上它确实是龟速,我们在乘法时强转成\(\_\_int128\),就可以避免溢出的情况了。
Code:
#pragma GCC optimize (2, "inline", "Ofast")
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int pri[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
int T;
ll n, ans = 0;
ll Qpow(ll x, ll y, ll Md) {
ll ans = 1;
for(; y; y >>= 1, x = x * x % Md) if(y & 1) ans = ans * x % Md;
return ans;
}
bool Miller_Rabin(ll x) {
if(x == 2) return 1;
if(!(x & 1) || x == 1) return 0;
ll s = 0, t = x - 1;
while(!(t & 1)) {
t >>= 1;
s ++;
}
for(int i = 0; i < 10 && pri[i] < x; i++) {
ll a = pri[i], b = Qpow(a, t, x);
ll k;
for(int j = 1; j <= s; j++) {
k = (__int128)b * b % x;
if(k == 1 && b != 1 && b != x - 1) return 0;
b = k;
}
if(b != 1) return 0;
}
return 1;
}
ll Gcd(ll a, ll b) {
return !b ? a : Gcd(b, a % b);
}
ll Pollard_Rho(ll n, ll c) {
ll x = 1ll * rand() * rand() % (n - 2) + 1, y = x;
int i = 1, k = 2;
while(1) {
i++;
x = (__int128)x * x % n;
x = (x + c) % n;
ll G = Gcd(y - x, n);
if(G > 1 && G < n) return G;
if(y == x) return n;
if(i == k) {
y = x;
k <<= 1;
}
}
}
void Find(ll x, ll c) {
if(x == 1 || x < ans) return ;
if(Miller_Rabin(x)) return (void) (ans = max(ans, x));
ll p = x;
while(p == x) p = Pollard_Rho(p, c--);
while(x % p == 0) x /= p;
Find(p, c); Find(x, c);
}
int main() {
srand(2333);
scanf("%d", &T);
while(T--) {
scanf("%lld", &n);
ans = 0;
bool f = Miller_Rabin(n);
if(f) {
puts("Prime");
continue;
}
Find(n, 1000000);
printf("%lld\n", ans);
}
return 0;
}
于是我们修改之后再次交一发:
\(nmdwsm!!!\)为什么这样就快了那么多啊……真是
然后发现最后一个点我们死活跑不过去md真毒瘤。然后我们继续考虑优化算法。发现整个算法的复杂度的瓶颈在于\(Pollard\_Rho\)中\(gcd\)上,虽然我们用了倍增算法减少了\(gcd\)的次数,但是这个复杂度仍然是无法承受的。但是需要注意到一件事情,我们只需要知道\(x\)的一个因子即可,并不需要做出规定一定是某个。然后一个显然的结论就是,根据欧几里得算法,\(gcd(a, b) > 1 \to gcd(a * c \% b, b) > 1\)。所以我们可以把一次倍增中所有的随机因子都乘起来之后,一起做一次\(gcd\),这样原本是互质的,现在仍然互质,原本存在公因数的,现在仍然存在。这样就可以大大减少\(gcd\)的次数,唯一需要处理的边界条件就是所有随机因子的乘积与\(n\)做\(gcd\)之后,\(gcd\)为\(n\),这个时候虽然存在公因数,但是却不能直接返回\(n\)(因为\(n\)是我们用来判断本次探测不成功的标记),这个时候我们再按照原方法枚举一遍,找到某个非平凡因子返回即可。
同时\(gcd\)的地方也是可以优化很多的,大致就是仿照高精度\(gcd\)那样用二进制处理,常数会小很多。
Code:
#pragma GCC optimize(3, "inline", "Ofast")
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int pri[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
int T;
ll n, ans = 0;
ll Qpow(ll x, ll y, ll Md) {
ll ans = 1;
for(; y; y >>= 1, x = x * x % Md) if(y & 1) ans = ans * x % Md;
return ans;
}
bool Miller_Rabin(ll x) {
if(x == 2) return 1;
if(!(x & 1) || x == 1) return 0;
ll s = 0, t = x - 1;
while(!(t & 1)) {
t >>= 1;
s ++;
}
for(int i = 0; i < 10 && pri[i] < x; i++) {
ll a = pri[i], b = Qpow(a, t, x);
ll k;
for(int j = 1; j <= s; j++) {
k = (__int128)b * b % x;
if(k == 1 && b != 1 && b != x - 1) return 0;
b = k;
}
if(b != 1) return 0;
}
return 1;
}
ll Gcd(ll a, ll b) {
if(!a || !b) return a | b;
int t = __builtin_ctzll(a | b); //__builtin_ctzll是得到某个数二进制下末尾0的个数
a >>= __builtin_ctzll(a);
do {
b >>= __builtin_ctzll(b);
if(b < a) swap(a, b);
b -= a;
}while(b);
return a << t;
}
ll Pollard_Rho(ll n, ll c) {
ll x = 1ll * rand() * rand() % (n - 2) + 1, y = x;
int i = 1, k = 2;
ll Mu = 1, z = y;
while(1) {
i++;
x = (__int128) x * x % n;
x = (x + c) % n;
Mu = (__int128) Mu * (y > x ? y - x : x - y) % n;
if(y == x) return n;
if(i == k) {
ll G = Gcd(Mu, n);
if(G == 1) { y = x; k <<= 1; continue; }
if(G == n) {
x = y;
for(int t = 1; t <= k >> 1; t++) {
x = (__int128) x * x % n;
x = (x + c) % n;
ll g = Gcd(y > x ? y - x : x - y, n);
if(g > 1 && g < n) return g;
}
}
return G;
}
}
}
void Find(ll x, ll c) {
if(x == 1 || x < ans) return ;
if(Miller_Rabin(x)) return (void) (ans = max(ans, x));
ll p = x;
while(p == x) p = Pollard_Rho(p, c--);
while(x % p == 0) x /= p;
Find(p, c); Find(x, c);
}
int main() {
srand(2333);
scanf("%d", &T);
while(T--) {
scanf("%lld", &n);
ans = 0;
bool f = Miller_Rabin(n);
if(f) {
puts("Prime");
continue;
}
Find(n, 100000);
printf("%lld\n", ans);
}
return 0;
}
尝试再次交一发,愉快的发现它过了。爆OJ真舒服
然后还有一个关于\(127\)这个数的优化,现在还没有看懂,咕咕咕~
原文地址:https://www.cnblogs.com/Apocrypha/p/10672354.html