luogu P4718 【模板】Pollard-Rho演算法(貼程式碼)
阿新 • • 發佈:2019-01-12
嘟嘟嘟
從標題中能看出來,我只是想貼一個程式碼。
先扯一會兒。
前幾天模擬考到了這東西,今天有空就學了一下。
到網上找資料,發現前置技能是miller-rabin篩法,於是我還得先學這麼個東西。
學miller-rabin的話不得不推薦這兩篇文章:
大數質因解:淺談Miller-Rabin和Pollard-Rho演算法
素數與素性測試(Miller-Rabin測試)
但是我還是看了好幾遍才懂……
至於pollard-rho這東西,還是上面的第一篇部落格,以及這一篇A Quick Tutorial on Pollard's Rho Algorithm(不是英文)。
但是某谷的板子特別毒瘤,因為他卡常。
首先得把龜速
然後這種floyd判圈法也過不了,於是我又找到了一個
剩下的就是一些
#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 n, ll a, ll d) { ll t = quickpow(a, d, n); while(d != n - 1 && t != n - 1 && t != 1) t = mul(t, t, n), d <<= 1; return t == n - 1 || (d & 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 x = rand() % (n - 2) + 2;//x屬於[2, n] if(!test(n, x, d)) 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 pollar_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 = pollar_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; }