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} floor * 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;
}