1. 程式人生 > 實用技巧 >P4718-[模板]Pollard-Rho演算法

P4718-[模板]Pollard-Rho演算法

正題

題目連結:https://www.luogu.com.cn/problem/P4718


題目大意

給出一個數\(n\),如果它是質數則輸出\(Prime\),否則輸出它的最大質因子。


解題思路

\(\text{Pollard-Rho}\)演算法的前置知識是\(\text{Miller-Rabin}\)。在使用\(\text{Miller-Rabin}\)判掉質數之後,\(\text{Pollard-Rho}\)使用基於隨機的思想能夠較快的求出一個大數的因子之一。

樸素的隨機演算法就是隨機一個數判斷它是不是因子,我們先使用一個較為優秀的隨機方式,\(f(x)=f(x-1)^2+c\)(其中\(c\)

為一個常數)。

然後我們利用在這個函式上“跑”的距離來判斷,也就是每次拿某兩個\(i,j\),判斷\(|f(i)-f(j)|\)是否為它的因數。

但是如果列舉的話\(f\)函式上會出現一些“環”,我們需要快速的判掉“環”的方法。每次拿\(s,t\),令\(t=2s\),若環長為\(c\),那麼有\(f(x)=f(x+c)\),當某一時刻\(f(t)=f(s)\)那麼環長一定是\(s\)的整數倍。

然後判到環就退出,如果沒有找到就換一個常數重新做,這樣的我們的演算法雛形就形成了。

但是這樣還是跑的很慢,發現我們在過程中大量呼叫了\(gcd(d,p)\)導致時間變慢。考慮優化,我們可以每次先做一堆,然後在把這一堆拿過去一起搞定。首先我們有\(gcd(ac,b)=gcd(a,b)\)

,然後根據\(gcd\)的原理,我們有\(gcd(a,b)=gcd(a\% b,b)\)那麼也就是我們有\(gcd(a,b)=gcd(ac\%b,b)\)

那麼假設我們有若干個間隔\(a_1,a_2,a_3,...\)那麼我們把這數乘起來模\(p\),然後把得到的結果\(k\)\(p\)\(gcd\)就等價於拿\(a\)中逐個取與\(p\)\(gcd\)

所以我們的優化方法就是第\(i\)次拿\(2^i\)個間隔去一起與\(p\)判斷,但是因為\(i\)後面會很大導致副作用,所以將\(i\)設一個上界\(8\)即可。
時間複雜度期望是\(O(n^{2.5})\),但跑的飛快

回到這題來,我們先對\(n\)

\(\text{MR}\)判斷一次質數,然後跑\(\text{Pr}\)弄出一個因子\(d\),之後將\(n\)的因子\(d\)都去光後分別把\(n\)\(d\)丟下去遞迴繼續跑。可以記錄一個目前最大質因子來剪去一些不優狀態。


code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#define ll long long
using namespace std;
const ll pri[10]={2,3,5,7,11,13,17,19,23,27};
ll T,n,ans;
ll ksc(ll a,ll b,ll p){
    ll c=(long double)a*b/p;
    long double ans=a*b-c*p;
    if(ans<0)ans+=p;
    else if(ans>=p)ans-=p;
    return ans;
}
ll power(ll x,ll b,ll p){
    ll ans=1;
    while(b){
        if(b&1)ans=ksc(ans,x,p);
        x=ksc(x,x,p);b>>=1;
    }
    return ans;
}
bool Mr(ll p){
    if(p==2)return 1;
    if(p<2||!(p&1))return 0;
    ll t=p-1,s=0;
    while(!(t&1))t>>=1,s++;
    for(ll i=0;i<10&&pri[i]<p;i++){
        ll x=power(pri[i],t,p),k;
        for(ll j=0;j<s;j++){
            k=ksc(x,x,p);
            if(k==1&&x!=1&&x!=p-1)
                return 0;
            x=k;
        }
        if(x!=1)return 0;
    }
    return 1;
}
ll gcd(ll a,ll b)
{return (!b)?a:gcd(b,a%b);}
ll Pr(ll p){
    ll s=0,t=0,c=1ll*rand()%(p-1)+1;
    for(ll g=1,val=1,d;;g<<=1,s=t,val=1){
        for(ll j=0;j<g;j++){
            t=(ksc(t,t,p)+c)%p;
            val=ksc(val,abs(t-s),p);
            if(j%127==0&&(d=gcd(p,val))>1)
                return d;
        }
        d=gcd(p,val);
        if(d>1)return d;
    }
    return p;
}
void solve(ll n){
    if(n<ans||n<2)return;
    if(Mr(n)){ans=n;return;}
    ll d=0;
    while((d=Pr(n))>=n);
    while(n%d==0)n/=d;
    solve(n);solve(d);
    return;
}
signed main()
{
    srand(998244353);
    scanf("%lld",&T);
    while(T--){
        scanf("%lld",&n);
        if(Mr(n)){
            printf("Prime\n");
            continue;
        }
        ans=0;solve(n);
        printf("%lld\n",ans);
    }
    return 0;
}