1. 程式人生 > >luogu P4718 【模板】Pollard-Rho演算法(貼程式碼)

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)\)的。
然後這種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 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;
}