luogu P4718 【模板】Pollard-Rho算法(贴代码)

时间:2023-03-08 23:07:36
luogu P4718 【模板】Pollard-Rho算法(贴代码)

嘟嘟嘟



从标题中能看出来,我只是想贴一个代码。

先扯一会儿。



前几天模拟考到了这东西,今天有空就学了一下。

到网上找资料,发现前置技能是miller-rabin筛法,于是我还得先学这么个东西。

学miller-rabin的话不得不推荐这两篇文章:

大数质因解:浅谈Miller-Rabin和Pollard-Rho算法

素数与素性测试(Miller-Rabin测试)

但是我还是看了好几遍才懂……



至于pollard-rho这东西,还是上面的第一篇博客,以及这一篇A Quick Tutorial on Pollard's Rho Algorithm(不是英文)。



但是某谷的板子特别毒瘤,因为他卡常。

首先得把龟速快速乘改成\(O(1)\)的。

顺便说一下,\(O(1)\)的快速乘的原理就是用\(a * b \ \ mod \ \ n = a * b - \lfloor \frac{a * b}{n} \rfloor * n\)这个式子。然后因为什么自然溢出后溢出的位的差只可能是0或1,所以判断一下这个结果是否小于0(这里我是真不懂)。

还有一点就是这种floyd判圈法也过不了,于是我又找到了一个相对好理解的倍增的pollard-rho,快的吓人。具体看代码吧。



剩下的就是一些玄学的细节了,代码里也有注释。

#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<vector>
#include<stack>
#include<queue>
#include<ctime>
using namespace std;
#define enter puts("")
#define space putchar(' ')
#define Mem(a, x) memset(a, x, sizeof(a))
#define In inline
typedef long long ll;
typedef double db;
const int INF = 0x3f3f3f3f;
const db eps = 1e-8;
//const int maxn = ;
inline ll read()
{
ll ans = 0;
char ch = getchar(), last = ' ';
while(!isdigit(ch)) {last = ch; ch = getchar();}
while(isdigit(ch)) {ans = (ans << 1) + (ans << 3) + ch - '0'; ch = getchar();}
if(last == '-') ans = -ans;
return ans;
}
inline void write(ll x)
{
if(x < 0) x = -x, putchar('-');
if(x >= 10) write(x / 10);
putchar(x % 10 + '0');
} ll n, Max = 0; In ll mul(ll a, ll b, ll mod)
{
ll d = (long double)a / mod * b + 1e-8; //还必须是long double……double精度不够
ll r = a * b - d * mod;
return r < 0 ? r + mod : r;
}
In ll quickpow(ll a, ll b, ll mod)
{
ll ret = 1;
for(; b; b >>= 1, a = mul(a, a, mod))
if(b & 1) ret = mul(ret, a, mod);
return ret;
} In bool test(ll a, ll d, ll n)
{
ll t = quickpow(a, d, n);
if(t == 1) return 1;
while(d != n - 1 && t != n - 1 && t != 1) t = mul(t, t, n), d <<= 1;
return t == n - 1; //这里就不用判1了,因为只可能在while前出现1的情况
}
int a[] = {2, 3, 5, 7, 11};
In bool miller_rabin(ll n)
{
if(n == 1) return 0;
for(int i = 0; i < 5; ++i) //先粗筛一遍
{
if(n == a[i]) return 1;
if(!(n % a[i])) return 0;
}
ll d = n - 1;
while(!(d & 1)) d >>= 1;
for(int i = 1; i <= 5; ++i) //搞五遍
{
ll a = rand() % (n - 2) + 2;//x属于[2, n - 1]
if(!test(a, d, n)) return 0;
}
return 1;
} In ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
In ll f(ll x, ll a, ll mod) {return (mul(x, x, mod) + a) % mod;}
const int M = (1 << 7) - 1; //我也不知道为啥M是这个数……
ll pollard_rho(ll n) //倍增版!减少gcd调用次数。(好像不用判环?)
{
for(int i = 0; i < 5; ++i) if(n % a[i] == 0) return a[i];
ll x = rand(), y = x, t = 1, a = rand() % (n - 2) + 2;
for(int k = 2;; k <<= 1, y = x)
{
ll q = 1;
for(int i = 1; i <= k; ++i)
{
x = f(x, a, n);
q = mul(q, abs(x - y), n);
// if(i >= M) //不等价!
if(!(i & M)) //超过了2 ^ 7,再用gcd
{
t = gcd(q, n);
if(t > 1) break; //找到了
}
}
if(t > 1 || (t = gcd(q, n)) > 1) break; //之所以再求一遍,是处理刚开始k < M的情况
}
return t;
}
In void find(ll x)
{
if(x == 1 || x <= Max) return;
if(miller_rabin(x)) {Max = max(Max, x); return;}
ll p = x;
while(p == x) p = pollard_rho(x);
while(x % p == 0) x /= p;
find(p); find(x);
} int main()
{
freopen("1.in", "r", stdin);
freopen("1.out", "w", stdout);
srand(time(0));
int T = read();
while(T--)
{
n = read();
Max = 0;
find(n);
if(Max == n) puts("Prime");
else write(Max), enter;
}
return 0;
}