【快速因数分解】Pollard's Rho 算法
Pollard-Rho 是一个很神奇的算法,用于在
简单来说:Pollard-Rho算法是 John Pollard发明的一种能 快速找到大整数的一个非1、非自身的因子 的算法。
前缀知识
PollardRhoPollardRho 算法:
问题模型:
给一个数
对于一个比较小的数(
对于一个非常大的合数
1.我们求 gcd
如果我们直接随机求
2.玄学随机法:
其实网上大多都叫这种方法为生日悖论(悖论是不符合经验但符合事实的结论),有兴趣的可以去看看。但博主也不懂觉得这个方法跟我们的 PollardRho 的关联。在1中我们用gcd的方案,现在我们考虑有没有办法能够快速的找到这些与
这个随机生成方法最终被
为一次操作,我们发现
判环
这种随机生成方法虽然优秀,但也有一些需要注意之处,比如有可能会生成一个环,并不断在这个环上生成以前生成过一次的数,所以我们必须写点东西来判环:
- 我们可以让
根据生成公式以比 快两倍的速度向后生成,这样当 再次与 相等时, 一定跑完了一个圈且没重复跑多少! - 我们可以用倍增的思想,让
记住 的位置,然后 再跑当前跑过次数的一倍的次数。这样不断让 记住 的位置, 再往下跑,因为倍增所以当 跑到 时,已经跑完一个圈,并且也不会多跑太多(这个必须好好想一想,也可以看看代码如何实现的)(因为这种方法会比第一种用的更多,因为它在 二次优化时还可以起到第二个作用)(这里就给第二种的代码吧)
inline ll rho(ll n) {
ll x, y, c, g;
rg i, j;
x = y = rand();
c = rand(); //初始化
i = 0, j = 1; //倍增初始化
while (++i) {
x = (ksc(x, x, n) + c) % n; // x只跑一次
if (x == y) break; //跑完了一个环
g = gcd(Abs(y - x), n); //求gcd(注意绝对值)
if (g > 1) return g;
if (i == j) y = x, j <<= 1; //倍增的实现
}
}
PollardRho 算法的二次优化:
我们发现其实
二次优化:我们发现我们在每一次生成操作中,乘法之后所模的数就是我们的
对原算法的一些影响:任何一个优化都是要考虑其对原算法的影响的。这个优化也会有一些影响:首先,因为我们会等大概127次后再去
那怎么办呢?(在上文就已经提到过了:倍增时我们的
inline ll rho(ll p) { //求出p的非平凡因子
ll x, y, z, c, g;
rg i, j; //先摆出来(z用来存(y-x)的乘积)
while (1) { //保证一定求出一个因子来
y = x = rand() % p; //随机初始化
z = 1;
c = rand() % p; //初始化
i = 0, j = 1; //倍增初始化
while (++i) { //开始玄学生成
x = (ksc(x, x, p) + c) % p; //可能要用快速乘
z = ksc(z, Abs(y - x), p); //我们将每一次的(y-x)都累乘起来
if (x == y || !z)
break; //如果跑完了环就再换一组试试(注意当z=0时,继续下去是没意义的)
if (!(i % 127) ||
i == j) { //我们不仅在等127次之后gcd我们还会倍增的来gcd
g = gcd(z, p);
if (g > 1) return g;
if (i == j)
y = x, j <<= 1; //维护倍增正确性,并判环(一箭双雕)
}
}
}
}
然后是两道道例题,都有点难度的。
code1:
#include <bits/stdc++.h>
using namespace std;
#define rg register int
typedef long long ll;
typedef unsigned long long ull;
typedef long double lb;
int n;
ll ans, m;
inline ll qr() { //快读
char ch;
while ((ch = getchar()) < '0' || ch > '9')
;
ll res = ch ^ 48;
while ((ch = getchar()) >= '0' && ch <= '9') res = res * 10 + (ch ^ 48);
return res;
}
inline ll Abs(ll x) { return x < 0 ? -x : x; } //取绝对值
inline ll gcd(ll x, ll y) { //非递归求gcd
ll z;
while (y) {
z = x;
x = y;
y = z % y;
}
return x;
}
inline ll ksc(ull x, ull y, ll p) { // O(1)快速乘(防爆long long)
return (x * y - (ull)((lb)x / p * y) * p + p) % p;
}
inline ll ksm(ll x, ll y, ll p) { //快速幂
ll res = 1;
while (y) {
if (y & 1) res = ksc(res, x, p);
x = ksc(x, x, p);
y >>= 1;
}
return res;
}
inline bool mr(ll x, ll p) { // mille rabin判质数
if (ksm(x, p - 1, p) != 1) return 0; //费马小定理
ll y = p - 1, z;
while (!(y & 1)) { //二次探测
y >>= 1;
z = ksm(x, y, p);
if (z != 1 && z != p - 1) return 0;
if (z == p - 1) return 1;
}
return 1;
}
inline bool prime(ll x) {
if (x < 2) return 0; // mille rabin判质数
if (x == 2 || x == 3 || x == 5 || x == 7 || x == 43) return 1;
return mr(2, x) && mr(3, x) && mr(5, x) && mr(7, x) && mr(43, x);
}
inline ll rho(ll p) { //求出p的非平凡因子
ll x, y, z, c, g;
rg i, j; //先摆出来(z用来存(y-x)的乘积)
while (1) { //保证一定求出一个因子来
y = x = rand() % p; //随机初始化
z = 1;
c = rand() % p; //初始化
i = 0, j = 1; //倍增初始化
while (++i) { //开始玄学生成
x = (ksc(x, x, p) + c) % p; //可能要用快速乘
z = ksc(z, Abs(y - x), p); //我们将每一次的(y-x)都累乘起来
if (x == y || !z)
break; //如果跑完了环就再换一组试试(注意当z=0时,继续下去是没意义的)
if (!(i % 127) ||
i == j) { //我们不仅在等127次之后gcd我们还会倍增的来gcd
g = gcd(z, p);
if (g > 1) return g;
if (i == j)
y = x, j <<= 1; //维护倍增正确性,并判环(一箭双雕)
}
}
}
}
inline void prho(ll p) { //不断的找他的质因子
if (p <= ans) return; //最优性剪纸
if (prime(p)) {
ans = p;
return;
}
ll pi = rho(p); //我们一次必定能求的出一个因子,所以不用while
while (p % pi == 0) p /= pi; //记得要除尽
return prho(pi), prho(p); //分开继续分解质因数
}
int main() {
// freopen("in.txt","r",stdin);
// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
n = qr();
srand(time(0)); //随机数生成必备!!!
for (rg i = 1; i <= n; ++i) {
ans = 1;
prho((m = qr()));
if (ans == m)
puts("Prime");
else
printf("%lld\n", ans);
}
return 0;
}
code2:
题目链接:杭电3864 D_num
根据这道题的答案性质,我们知道这个数在除了1以外还有三个约数,仔细想一下不难发现这类数最多只能有两个质因子,而这两个质因子分别就是它的两个约数,而第三个约数自然就是它自己了!(当然还有可能出现质数的三次方这类特殊的数,所以我们要特判一下)
所以我们的思路就是,先求出
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define db double
#define rg register int
#define cut {puts("is not a D_num");continue;} //这个可以看一看
ll n;
inline ll qr() {
char ch;
while ((ch = getchar()) < '0' || ch > '9')
;
ll res = ch ^ 48;
while ((ch = getchar()) >= '0' && ch <= '9') res = res * 10 + (ch ^ 48);
return res;
}
inline ll Abs(ll x) { return x < 0 ? -x : x; }
inline ll gcd(ll x, ll y) {
ll z;
while (y) {
z = x, x = y, y = z % y;
}
return x;
}
inline ll ksm(ll x, ll y, ll p) {
ll res = 1;
while (y) {
if (y & 1) res = (__int128)res * x % p;
x = (__int128)x * x % p;
y >>= 1;
}
return res;
}
inline bool mr(ll x, ll p) {
if (ksm(x, p - 1, p) != 1) return 0;
ll y = p - 1, z;
while (!(y & 1)) {
y >>= 1;
z = ksm(x, y, p);
if (z != 1 && z != p - 1) return 0;
if (z == p - 1) return 1;
}
return 1;
}
inline bool prime(ll p) {
if (p < 2) return 0;
if (p == 2 || p == 3 || p == 5 || p == 7 || p == 43) return 1;
return mr(2, p) && mr(3, p) && mr(5, p) && mr(7, p) && mr(43, p);
}
inline ll rho(ll p) {
ll x, y, z, c, g;
rg i, j;
while (1) {
y = x = rand() % p;
z = 1, c = rand() % p;
i = 0, j = 1;
while (++i) {
x = ((__int128)x * x + c) % p;
z = (__int128)z * Abs(y - x) % p;
if (x == y) break;
if (i % 72 || i == j) {
g = gcd(z, p);
if (g > 1 && g != p) return g;
if (i == j) y = x, j <<= 1;
}
}
}
}
int main() {
// freopen("in.txt","r",stdin);
// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
srand(time(0));
ll p, q;
while (scanf("%lld", &n) != EOF) {
if (prime(n) || n <= 2) cut;
p = rho(n);
q = n / p;
if ((__int128)p * p == q || (__int128)q * q == p) //这题很卡细节的
{
printf("%lld %lld %lld\n", p, q, n);
continue;
}
if (!prime(p) || !prime(q) || p == q) cut; //注意两个相等也不行
printf("%lld %lld %lld\n", min(p, q), max(p, q), n);
}
return 0;
}
参考
https://blog.csdn.net/qq_39972971/article/details/82346390
· 探究高空视频全景AR技术的实现原理
· 理解Rust引用及其生命周期标识(上)
· 浏览器原生「磁吸」效果!Anchor Positioning 锚点定位神器解析
· 没有源码,如何修改代码逻辑?
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· Obsidian + DeepSeek:免费 AI 助力你的知识管理,让你的笔记飞起来!
· 分享4款.NET开源、免费、实用的商城系统
· 解决跨域问题的这6种方案,真香!
· 一套基于 Material Design 规范实现的 Blazor 和 Razor 通用组件库
· 5. Nginx 负载均衡配置案例(附有详细截图说明++)