数论模板

Miller-Rabin 素性测试 & Pollard-Rho

#include <random>

const int TEST_TIME=16;
random_device rd;
mt19937_64 gen(rd());

ll bmul(ll a,ll b,ll m){
  ull c=(ull)a*(ull)b-(ull)((ld)a/m*b+0.5L)*(ull)m;
  if(c<(ull)m){
      return c;
  }
  return c+m;
}

ll qpow(ll a,ll b,ll P){
    ll res=1;
    while(b>0){
        if(b&1){
            res = bmul(res,a,P);
        }
        a = bmul(a,a,P);
        b >>= 1;
    }
    return res;
}

bool MillerRabin(ll n){
    if(n<2){
        return false;
    }
    if(n==2||n==3){
        return true;
    }

    ll u=n-1,t=0;
    while(u%2==0){
        u /= 2;
        t++;
    }

    uniform_int_distribution<ll> dist(0,n-4);

    for(int i=1;i<=TEST_TIME;i++){
        ll a=dist(gen)+2,v=qpow(a,u,n),s=0;
        if(v==1){
            continue;
        }
        for(s=0;s<t;s++){
            if(v==n-1){
                break;
            }
            v = bmul(v,v,n);
        }
        if(s==t){
            return false;
        }
    }

    return true;
}

ll gcd(ll a,ll b){
    return b==0?a:gcd(b,a%b);
}

ll PollardRho(ll x){
    ll s=0,t=0,c=uniform_int_distribution<>(0,x-2)(gen)+1,val=1;

    for(int gl=1;;gl*=2,s=t,val=1){
        for(int stp=1;stp<=gl;stp++){
            t = (bmul(t,t,x)+c)%x;
            val = bmul(val,abs(t-s),x);
            if((stp%127)==0){
                ll d=gcd(val,x);
                if(d>1){
                    return d;
                }
            }
        }
        ll d=gcd(val,x);
        if(d>1){
            return d;
        }
    }
}

配套

ll mxfac;
void maxFac(ll x,bool fst=true){
    if(fst){
        mxfac = 0;
    }

    if(x<=mxfac||x<2){
        return;
    }

    if(MillerRabin(x)){
        mxfac = max(mxfac,x);
        return;
    }

    ll p=x;
    while(p>=x){
        p = PollardRho(x);
    }
    while((x%p)==0){
        x /= p;
    }

    maxFac(x,false);
    maxFac(p,false);
}

#include <vector>

vector<ll> facs;
void getFacs(ll x,bool fst=true){
    if(fst){
        facs.clear();
    }

    if(x==1){
        return;
    }

    if(MillerRabin(x)){
        facs.push_back(x);
        return;
    }

    ll p=x;
    while(p>=x){
        p = PollardRho(x);
    }
    while((x%p)==0){
        x /= p;
    }

    getFacs(x,false);
    getFacs(p,false);
}