数论模板

exgcd 求解逆元

// ax%b==1
exgcd(a,b,x,y);
ans = (x+b)%b;

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);
}

CRT

void exgcd(ll a,ll b,ll &x,ll &y){
    if(b==0){
        x = 1;
        y = 0;
    }
    else{
        exgcd(b,a%b,y,x);
        y -= a/b*x;
    }
}

ll crt(int k,ll* a,ll* r){
    ll n=1,ans=0;
    for(int i=1;i<=k;i++){
        n = n*a[i];
    }
    for(int i=1;i<=k;i++){
        ll m=n/a[i],b,y;
        exgcd(m,a[i],b,y);
        ans = (ans+r[i]*m*b%n)%n;
    }
    return (ans%n+n)%n;
}

Lagrange Interpolation

#include <vector>

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

ll inv(ll a,ll P){
    return qpow(a,P-2,P);
}

vector<int> lagrangeInterpolation(const vector<int> &x,
                                  const vector<int> &y,
                                  ull P){
    int n=x.size();
    vector<int> M(n+1),px(n,1),f(n);
    M[0] = 1;

    for(int i=0;i<n;i++){
        for(int j=i;j>=0;j--){
            M[j+1] = (M[j]+M[j+1])%P;
            M[j] = (ll)M[j]*(P-x[i])%P;
        }
    }

    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            if(i!=j){
                px[i] = (ll)px[i]*(x[i]-x[j]+P)%P;
            }
        }
    }

    for(int i=0;i<n;i++){
        ll t=(ll)y[i]*inv(px[i],P)%P,k=M[n];
        for(int j=n-1;j>=0;j--){
            f[j] = (f[j]+k*t)%P;
            k = (M[j]+k*x[i])%P;
        }
    }

    return f;
}

Lucas

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

ll inv(ll a,ll P){
    return qpow(a,P-2,P);
}

ll fac[100005],ifac[100005];

void init(int P){
    fac[0] = 1;
    for(int i=1;i<P;i++){
        fac[i] = fac[i-1]*i%P;
    }
    ifac[P-1] = inv(fac[P-1],P);
    for(int i=P-1;i>=1;i--){
        ifac[i-1] = ifac[i]*i%P;
    }
}

ll C(int n,int m,int P){
    if(n<m){
        return 0;
    }
    if(n<P&&m<P){
        return fac[n]*ifac[n-m]%P*ifac[m]%P;
    }
    return C(n/P,m/P,P)*C(n%P,m%P,P)%P;
}

// -------------------------------------------

class BinomModPrime{
    int p;
    vector<int> fac,ifac;

    int calc(int n,int k){
        if(n<k){
            return 0;
        }
        ll res=fac[n];
        res = (res*ifac[k])%p;
        res = (res*ifac[n-k])%p;
        return res;
    }

public:
    BinomModPrime(int p):p(p),fac(p),ifac(p){
        fac[0] = 1;
        for(int i=1;i<p;i++){
            fac[i] = (ll)fac[i-1]*i%p;
        }
        ifac[p-1] = p-1;
        for(int i=p-1;i>=1;i--){
            ifac[i-1] = (ll)ifac[i]*i%p;
        }
    }

    int binomial(ll n,ll k){
        ll res=1;
        while(n>0||k>0){
            res = (res*calc(n%p,k%p))%p;
            n /= p;
            k /= p;
        }
        return res;
    }

    int operator()(const ll &n,const ll &k){
        return binomial(n,k);
    }
};