[学习笔记]Pollard-Rho

之前学的都是假的

%%zzt

Miller_Rabin:Miller-Rabin与二次探测

大质数分解：

找到所有质因子，再logn搞出质因子的次数

方法：不断找到一个约数d，递归d,n/d进行分解，直到n是质数

快速幂快速乘：

ll qk(ll a,ll b,ll m){
    ll d=((long double)a/m*b);
    ll r=a*b-d*m;
    return ((ull)r+m)%m;
}
ll qm(ll x,ll y,ll mod){
    ll ret=1;
    while(y){
        if(y&1) ret=qk(ret,x,mod);
        x=qk(x,x,mod);
        y>>=1;
    }
    return ret;
}

注意快速乘：((ull)r+m)%m由于r可能<0或者>m，这一步是必须的

Miller-Rabin：

bool M_R(ll p){
    if(p==1) return false;
    if(p==2||p==3||p==5||p==7||p==11||p==61) return true;
    if(p%2==0||p%3==0||p%5==0||p%7==0) return false;
    int s=ctz(p-1);
    for(reg i=0;i<6;++i){
        int a=pri[i];
        ll k=(p-1)>>s;
        k=qm(a,k,p);
        if(k==1) continue;
        ll las=k;
        for(reg j=1;j<=s;++j){
            k=qk(k,k,p);
            if(k==1&&las!=p-1&&las!=1) return false;
            las=k;
        }
        if(k!=1) return false;
    }
    return true;
}

二进制gcd

ll gcd(ll a,ll b)
{
    if(!a||!b) return a|b;
    #define ctz __builtin_ctzll
    int shift=ctz(a|b);
    b>>=shift;
    while(a)
    {
        a>>=ctz(a);
        if(a<b)
            swap(a,b);
        a-=b;
    }
    return b<<shift;
    #undef ctz
}

就是高精gcd才用的更相减损术

Pollard-Rho

主体1

ll P_R(ll p,ll c){
    ll x=0,y=0,d;
    int k=2,has=1;
    ll tmp=1;
    while(1){
        ++has;
        x=ad(qk(x,x,p),c,p);
        tmp=qk(tmp,abs(y-x),p);
        if(x==y) return p;
        if(k==has){
            k<<=1;
            y=x;
            d=gcd(tmp,p);
            tmp=1;
            if(d>1) return d;
        }
    }
}

注意tmp，把所有的abs(y-x)乘在一起，gcd显然不变，并且减少了求gcd次数

主体2

void fin(ll p,int cnt){
    if(p==1||p<=ans) return;
    if(M_R(p)) {
        ans=max(ans,p);return;
    }
    ll d=P_R(p,cnt);
    while(d==p) --cnt,d=P_R(p,cnt);
    while(p%d==0) p/=d;
    fin(p,cnt);fin(d,cnt);
}

如果要求质因子，开个vector，最后去重即可

模板：

【模板】Pollard-Rho算法 

// luogu-judger-enable-o2
#include<bits/stdc++.h>
#define reg register int
#define il inline
#define fi first
#define se second
#define mk(a,b) make_pair(a,b)
#define numb (ch^'0')
#define pb push_back
#define solid const auto &
#define enter cout<<endl
#define pii pair<int,int>
using namespace std;
typedef long long ll;
template<class T>il void rd(T &x){
    char ch;x=0;bool fl=false;while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);
    for(x=numb;isdigit(ch=getchar());x=x*10+numb);(fl==true)&&(x=-x);}
template<class T>il void output(T x){if(x/10)output(x/10);putchar(x%10+'0');}
template<class T>il void ot(T x){if(x<0) putchar('-'),x=-x;output(x);putchar(' ');}
template<class T>il void prt(T a[],int st,int nd){for(reg i=st;i<=nd;++i) ot(a[i]);putchar('
');}
namespace Modulo{
const int mod=998244353;
int ad(int x,int y){return (x+y)>=mod?x+y-mod:x+y;}
void inc(int &x,int y){x=ad(x,y);}
int mul(int x,int y){return (ll)x*y%mod;}
void inc2(int &x,int y){x=mul(x,y);}
int qm(int x,int y=mod-2){int ret=1;while(y){if(y&1) ret=mul(x,ret);x=mul(x,x);y>>=1;}return ret;}
}
//using namespace Modulo;
namespace Miracle{
ll ans;
ll qk(ll x,ll y,ll mod){
    ll d=((long double)x*y/mod);
    ll r=x*y-mod*d;
    return r<0?r+mod:(r>=mod?r-mod:r);
}
ll qm(ll x,ll y,ll mod){
    ll ret=1;
    while(y){
        if(y&1) ret=qk(ret,x,mod);
        x=qk(x,x,mod);
        y>>=1;
    }
    return ret;
}
ll ad(ll x,ll y,ll mod){
    return (x+y)>=mod?x+y-mod:x+y;
}
ll gcd(ll a,ll b){
    if(!a||!b) return a+b;
    #define ctz __builtin_ctzll
    int tmp=ctz(a|b);
    b>>=ctz(b);
    while(a){
        a>>=ctz(a);
        if(a<b) swap(a,b);
        a-=b;
    }
    return b<<tmp;
}
int pri[6]={2,3,5,7,11,61};
bool M_R(ll p){
    // cout<<"M_R "<<p<<endl;
    if(p==1) return false;
    if(p==2||p==3||p==5||p==7||p==11||p==61) return true;
    if(p%2==0||p%3==0||p%5==0||p%7==0) return false;
    int s=ctz(p-1);
    for(reg i=0;i<6;++i){
        ll a=pri[i];
        ll tmp=(p-1)>>s;
        ll now=qm(a,tmp,p);
        if(now==1||now==p-1) continue;
        for(reg j=1;j<=s;++j){
            ll las=now;
            now=qk(now,now,p);
            if(now==1&&las!=p-1&&las!=1) return false;
        }
        if(now!=1) return false;
    }
    return true;
}
ll P_R(ll p,ll c){
    // cout<<" P_R "<<p<<" "<<c<<endl;
    ll x,y,d;
    ll tmp=1,has=1,k=2;
    x=y=0;
    while(1){
        ++has;
        x=ad(qk(x,x,p),c,p);
        tmp=qk(tmp,abs(y-x),p);
        if(x==y) return p;
        if(has==k){
            k<<=1;
            d=gcd(tmp,p);
            if(d>1) {
                // cout<<" ret "<<d<<endl;
                return d;
            }
            tmp=1;
            y=x;
        }
    }
}
void sol(ll n,int c){
    if(n==1||n<=ans) return;
    if(M_R(n)){
        ans=max(ans,n);return;
    }
    ll d=P_R(n,c);
    while(d==n) --c,d=P_R(n,c);
    while(n%d==0) n/=d;
    sol(n,c);sol(d,c);
}
int main(){
    int t;
    rd(t);
    srand(19260817);
    while(t--){
        ans=0;ll n;rd(n);
        sol(n,19260817);
        if(ans!=n) printf("%lld
",ans);
        else printf("Prime
");
    }
    return 0;
}

}
signed main(){
    Miracle::main();
    return 0;
}

/*
   Author: *Miracle*
*/