数论模板

组合数

　　递推

void get_C()
{
    for(int i=0;i<=2005;i++)
    {
        C[0][i]=C[i][i]=1;
        for(int j=1;j<i;j++)
        {
            C[j][i]=(C[j-1][i-1]+C[j][i-1])%MOD;
        }
    }
}

　　虎爷版

int f[10000+10];LL pf[10000+10];
LL qmod(LL a, LL b)
{
    if (b == 0) return 1;
    LL r = a % MOD;
    LL k = 1;
    while (b > 1){
        if ((b & 1)!=0)
            k = (k * r) % MOD;
        r = (r * r) % MOD;
        b >>= 1;
    }
    return (r * k) % MOD;
}
void init(){
    f[0]=1;
    for(int i=1;i<=maxn;i++){
        f[i]=(long long)f[i-1]*i%MOD;
        pf[i]=qmod(f[i],MOD-2);
    }
}

int C(int n,int m){
    if(n==m||m==0)return 1;
    return ((long long)f[n]*qmod(f[m],MOD-2)%MOD)*qmod(f[n-m],MOD-2)%MOD;
}

 　　杜教版

int fac[maxn], inv[maxn];
long long qmod(int a, int b)
{
    long long ret = 1;
    for(; b; b >>= 1, a = (long long)a * a % MOD)
        if(b & 1)
            ret = ret * a % MOD;
    return ret;
}
void init()
{
    fac[0] = 1;
    for(int i = 1; i < maxn; i++)
        fac[i] = (long long)fac[i - 1] * i % MOD;
    inv[maxn - 1] = qmod(fac[maxn - 1], MOD - 2);
    for(int i = maxn - 1; i > 0; i--)
        inv[i - 1] = (long long) inv[i] * i % MOD;
}
long long C(int n, int m)
{
    if(n < 0 || m < 0 || m > n)
        return 0;
    return (long long)fac[n] * inv[m] % MOD * inv[n - m] % MOD;
}

埃氏筛法

const int maxn = 1e+6 + 7;
bool prime[maxn];
int rec[maxn], cnt;
void init_prime()
{
    cnt = 0;
    memset(prime, true, sizeof(prime));
    prime[0] = prime[1] = false;///表明true为质数
    for (int i = 2; i <= maxn; ++i){
        if (prime[i]) rec[cnt++] = i;
        //此处边界判断为rec[j] <= maxn / i,如果写成i * rec[j] <= maxn,需要确保i * rec[j]不会溢出int
        for (int j = 0; j < cnt && rec[j] <= maxn / i; ++j){
            prime[i * rec[j]] = false;
            if (i % rec[j] == 0)
                break;
        }
    }
}

快速幂

　　取模

LL qmod(LL a, LL b)
{
    if (b == 0) return 1;
    LL r = a % MOD;
    LL k = 1;
    while (b > 1){
        if ((b & 1)!=0)
            k = (k * r) % MOD;
        r = (r * r) % MOD;
        b >>= 1;
    }
    return (r * k) % MOD;
}

　　不取模

LL qmod(LL a, LL b)
{
    if (b == 0) return 1;
    LL r = a;
    LL k = 1;
    while (b > 1){
        if ((b & 1)!=0)
            k = (k * r);
        r = (r*r);
        b >>= 1;
    }
    return r*k;
}

　　传模数

LL qmod(LL a, LL b,LL p)
{
    if (b == 0) return 1;
    LL r = a % p;
    LL k = 1;
    while (b > 1){
        if ((b & 1)!=0)
            k = (k * r) % p;
        r = (r * r) %p;
        b >>= 1;
    }
    return (r * k) % p;
}

快速乘

　　lgn

LL mult(LL A,LL B)
{
    LL z = 0;
    if (B == 0) return z;
    z = mult(A,B >> 1);
    z = (z << 1) % mod;
    if (B & 1) z = (z + A) %mod;
    return z;
}

　　O(1)

//O(1)快速乘
inline LL quick_mul(LL x,LL y,LL MOD){
    x=x%MOD,y=y%MOD;
    return ((x*y-(LL)(((long double)x*y+0.5)/MOD)*MOD)%MOD+MOD)%MOD;
}

自适应积分

double const eps=1e-6;
double F(double x)
{
    return 1.0/sqrt(1.0+x*x);
}
double simpson(double a,double b)
{
    double c = a + (b-a)/2;
    return (F(a)+4*F(c) + F(b)) * (b-a)/6;
}
double asr(double a,double b,double eps,double A)
{
    double c = a + (b-a)/2;
    double L = simpson(a,c),R = simpson(c,b);
    if(fabs(L + R - A) <= 15*eps)return L + R + (L + R - A)/15.0;
    return asr(a,c,eps/2,L) + asr(c,b,eps/2,R);
}
double asr(double a,double b,double eps)
{
    return asr(a,b,eps,simpson(a,b));
}
///printf("%.6f
",asr(0,n,eps));

欧拉函数

　　筛法

#define maxn 1000001
int euler[maxn];
void Init(){
     euler[1]=1;
     for(int i=2;i<maxn;i++)
       euler[i]=i;
     for(int i=2;i<maxn;i++)
        if(euler[i]==i)
           for(int j=i;j<maxn;j+=i)
              euler[j]=euler[j]/i*(i-1);//先进行除法是为了防止中间数据的溢出
}

　　求单值

//直接求解欧拉函数  
int euler(int n){ //返回euler(n)   
     int res=n,a=n;  
     for(int i=2;i*i<=a;i++){  
         if(a%i==0){  
             res=res/i*(i-1);//先进行除法是为了防止中间数据的溢出   
             while(a%i==0) a/=i;  
         }  
     }  
     if(a>1) res=res/a*(a-1);  
     return res;  
}  

扩展欧几里得

　　常规扩欧

LL ex_gcd(LL a,LL b,LL &x,LL &y)
{
    if(b==0){
        x=1;y=0;
        return a;
    }
    LL r=ex_gcd(b,a%b,y,x);
    y-=x*(a/b);
    return r;
}

　　扩欧求逆元

int exgcd(int a,int b,int &x,int &y)
{
    if(b==0){x=1,y=0;return a;}
    int ans=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return ans;
}
int inv(int a,int p)
{
    int x,y;
    int g=ecgcd(a,p,x,y);
    if(1%g!=0) return -1;
    x*=1/g;
    p=abs(p);
    int ans=x%p;
    if(ans<=0)ans+=p;
    return ans;
}

牛顿迭代求开方

const double eps=1e-9;
double sqr(double x)
{
    double k=x;
    while(k*k-x>eps)
        k=0.5*(k+x/k);
    return k;
}

错排公式

///f[i] pf[i]代表阶乘表和阶乘逆元表
LL solve(int x)
{
    LL res=0;
    for(int i=2;i<=x;i++){
        LL as=(f[x]*pf[i])%MOD;
        if(i&1){
            res=(res-as+MOD)%MOD;
        }
        else res=(res+as)%MOD;
    }
    return res;
}

第二类斯特林数

int s2[maxn][maxn];
void init_s2()
{///s2[i][j]代表把i个数划分到j个集合中
    for(int i=1;i<maxn;i++){
        s2[i][1]=s2[i][i]=1;
        for(int j=2;j<i;j++){
            s2[i][j]=(s2[i-1][j-1]+(long long)j*s2[i-1][j]%MOD)%MOD;
        }
    }
}

 质因数分解

void factor(LL n,LL a[maxn],LL b[maxn],int &tot) {
    LL temp,now,i;
    temp=(int)((double)sqrt(n)+1);
    tot=0;
    now=n;
    for(i=2;i<=temp;i++) {
        if(now%i==0) {
            a[++tot]=i;
            b[tot]=0;
            while(now%i==0) {
                ++b[tot];
                now/=i;
            }
        }
    }
    if(now!=1) {
        a[++tot]=now;
        b[tot]=1;
    }
}

FFT

　　namespace封装

namespace FFT {
    struct comp {
        double r, i;
        explicit comp(double r = 0.0, double i = 0.0) : r(r), i(i) {}
        inline friend comp operator +(comp a, comp b) {
            return comp(a.r + b.r, a.i + b.i);
        }
        inline friend comp operator -(comp a, comp b) {
            return comp(a.r - b.r, a.i - b.i);
        }
        inline friend comp operator *(comp a, comp b) {
            return comp(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
        }
        inline friend comp operator /(comp a, double b) {
            return comp(a.r / b, a.i / b);
        }
    };
    const double pi = acos(-1);
    const int N = 5e5 + 10;
    int rev[N];
    inline void fft(comp *p, int n, int idft) {
        for (int i = 0; i < n; i++)
            if (i < rev[i])
                std::swap(p[i], p[rev[i]]);
        for (int j = 1; j < n; j <<= 1) {
            static comp wn1, w, t0, t1;
            wn1 = comp(cos(pi / j), idft * sin(pi / j));
            for (int i = 0; i < n; i += j << 1) {
                w = comp(1, 0);
                for (int k = 0; k < j; k++) {
                    t0 = p[i + k];
                    t1 = w * p[i + j + k];
                    p[i + k] = t0 + t1;
                    p[i + j + k] = t0 - t1;
                    w = w * wn1;
                }
            }
        }
        if (idft == -1) {
            for (int i = 0; i < n; i++)
                p[i] = p[i] / n;
        }
    }
    template <typename T>
    inline T* fft_main(T *a, T *b, int n, int m) {
        static T res[N];
        static int nn, len;
        static comp aa[N], bb[N];
        len = 0;
        for (nn = 1; nn < m + n; nn <<= 1)
            len++;
        for (int i = 0; i < nn; i++) {
            aa[i] = comp(a[i], 0);
            bb[i] = comp(b[i], 0);
        }
        rev[0] = 0;
        for (int i = 1; i < nn; i++)
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
        fft(aa, nn, 1);
        fft(bb, nn, 1);
        for (int i = 0; i < nn; i++)
            aa[i] = aa[i] * bb[i];
        fft(aa, nn, -1);
        for (int i = 0; i < nn; i++)
            res[i] = aa[i].r + 0.5;
        return res;
    }
}

　　手撕FFT（kuangbin版，不对逆元进行打表）

const double PI=acos(-1.0);
struct cp
{
    double r,i;
    cp(){}
    cp(double _r,double _i)
    {
        r=_r;i=_i;
    }
    cp operator +(const cp &a)const
    {
        return cp(a.r+r,a.i+i);
    }
    cp operator -(const cp &a)const
    {
        return cp(r-a.r,i-a.i);
    }
    cp operator *(const cp &a)const
    {
        return cp(r*a.r-i*a.i,r*a.i+i*a.r);
    }
    cp conj()
    {
        return cp(r,-i);
    }
};
void change(cp y[],int len)
{
    int i,j,k;
    for(i = 1, j = len/2;i < len-1;i++){
        if(i < j)swap(y[i],y[j]);
        k = len/2;
        while( j >= k){
            j -= k;
            k /= 2;
        }
        if(j < k)j += k;
    }
}
void fft(cp y[],int len,int on)
{///on==1代表正fft，on==-1代表逆fft
    change(y,len);
    for(int h = 2;h <= len;h <<= 1){
        cp wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j = 0;j < len;j += h){
            cp w(1,0);
            for(int k = j;k < j+h/2;k++){
                cp u = y[k];
                cp t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)
        for(int i = 0;i < len;i++)
            y[i].r /= len;
    /**
    for(int i=0;i<as;i++){
            num[i]=(long long)(x[i].r+0.5);
    } ///整数最后需要这样取整,fft中最后已经还原nn
    */
}

　　胡小兔（对逆元进行打表）

const double PI=acos(-1.0);
struct cp
{
    double r,i;
    cp(){}
    cp(double _r,double _i)
    {
        r=_r;i=_i;
    }
    cp operator +(const cp a)const
    {
        return cp(a.r+r,a.i+i);
    }
    cp operator -(const cp a)const
    {
        return cp(r-a.r,i-a.i);
    }
    cp operator *(const cp a)const
    {
        return cp(r*a.r-i*a.i,r*a.i+i*a.r);
    }
    cp conj()
    {
        return cp(r,-i);
    }
};
int n=1,m;
double q[N],res[N];
cp f[N],f1[N],g[N],omg[N],inv[N];
void FFT_init(){
    for(int i=0;i<n;i++){
        omg[i]=cp(cos(2*PI*i/n),sin(2*PI*i/n));
        inv[i]=omg[i].conj();
    }
}
/**
在n倍增为长度2m长度时while(n<2*m)n*=2;
才可以进行FFTinit；正向FFT时传omg,逆fft传inv
计算结束后需要手动除n还原
*/
void fft(cp *a,cp *omg){
    int lim=0;
    while((1 << lim) < n) lim++;
    for(int i=0;i<n;i++){
        int t=0;
        for(int j=0;j<lim;j++)
            if(i >> j&1) t |= 1 << (lim - j - 1);
        if(i<t) swap(a[i],a[t]);
    }
    for(int l=2;l<=n;l*=2){
        int m=l/2;
        for(cp *p=a;p!=a+n;p+=l){
            for(int i=0;i<m;i++){
                cp t=omg[n/l*i]*p[m+i];
                p[m+i]=p[i]-t;
                p[i]=p[i]+t;
            }
        }
    }
}

FWT

void fwt(ll f[], int len, int op) {
    int n = (1 << len);
    for (int i = 1; i <= len; ++i) {
        int m = (1 << i), len = m >> 1;
        for (int r = 0; r < n; r += m) {
            int t1 = r, t2 = r + len;
            for (int j = 0; j < len; ++j, ++t1, ++t2) {
                ll x1 = f[t1], x2 = f[t2];
                if (op == 1) {   //xor
                    f[t1] = x1 + x2;
                    f[t2] = x1 - x2;
                    //if(f[t1] >= MOD) f[t1] -= MOD;
                    //if(f[t2] < 0) f[t2] += MOD;
                }
                if (op == 2) {   //and
                    f[t1] = x1 + x2;
                    f[t2] = x2;
                    //if(f[t1] >= MOD) f[t1] -= MOD;
                }
                if (op == 3) {   //or
                    f[t1] = x1;
                    f[t2] = x2 + x1;
                    //if(f[t2] >= MOD) f[t2] -= MOD;
                }
            }
        }
    }
}
void ifwt(ll f[], int len, int op) {
    int n = (1 << len);
    for (int i = len; i >= 1; --i) {
        int m = (1 << i), len = m >> 1;
        for (int r = 0; r < n; r += m) {
            int t1 = r, t2 = r + len;
            for (int j = 0; j < len; ++j, ++t1, ++t2) {
                ll x1 = f[t1], x2 = f[t2];
                if (op == 1) {   //xor
                    f[t1] = (x1 + x2) / 2;
                    f[t2] = (x1 - x2) / 2;
                    //f[t1] = (x1 + x2) * inv2;
                    //f[t2] = (x1 - x2) * inv2;
                    //if(f[t1] >= MOD) f[t1] %= MOD;
                    //if(f[t2] >= MOD) f[t2] %= MOD;
                    //if(f[t2] < 0) f[t2] = f[t2] % MOD + MOD;
                }
                if (op == 2) {   //and
                    f[t1] = x1 - x2;
                    f[t2] = x2;
                    //if(f[t1] < 0) f[t1] += MOD;
                }
                if (op == 3) {   //or
                    f[t1] = x1;
                    f[t2] = x2 - x1;
                    //if(f[t2] < 0) f[t2] += MOD;
                }
            }
        }
    }
}

中国剩余定理

　　　互质

void exgcd(int a,int b,int &x,int &y)
{
    if (b==0){x=1;y=0;return ;}
    exgcd (b,a%b,x,y);
    int t=x;x=y;y=t-a/b*y;
}
int CRT (int a[],int b[],int nn)
{
///a为模数,b为余数，poj1006
    int res=0,x,y;
    for(int i=1;i<=3;i++){
        int as=nn/a[i];
        exgcd(as,a[i],x,y);
        res+=as*x*b[i];
    }
    return res;
}

　　　　不互质（未测试）

inline LL gcd(LL a,LL b){
    return b?gcd(b,a%b):a;
}
inline LL lcm(LL a,LL b){
    return a*b/gcd(a,b);
}
inline void exgcd(LL a,LL b,LL &x,LL &y){//a,b,x,y同ax+by=gcd(a,b)中的a,b,x,y
    if(!b){
        x=1,y=0;return;
    }
    LL t;
    exgcd(b,a%b,x,y);
    t=x,x=y,y=t-(a/b)*y;
}
inline bool merge(LL a1,LL m1,LL a2,LL m2,LL &a3,LL &m3){
    //将方程x=a1+k1m1和x=a2+k2m2合并为x=a3+k3m3;
    LL d=gcd(m2,m1),a=a2-a1,q,y;
    if(a%d){
        return false;
    }//无解
    m3=lcm(m1,m2);
    exgcd(m1/d,m2/d,q,y);
    a3=a/d*q*m1+a1;
    ((a3%=m3)+=m3)%=m3;
    return true;
}
inline LL CRT(LL a[],LL m[],int n){
    ///m为模数
    LL a1=a[1],m1=m[1];
    bool ok;
    for(int i=2;i<=n;i++){
        ok=merge(a1,m1,a[i],m[i],a1,m1);
        if(!ok)break;
    }
    if(!ok)return -1;
    return (a1%m1+m1)%m1;//返回最小正整数解
}

 线性求1~n中所有数%p下的逆元

void Inv(int p,int a[],int n){
//线性求<=n的数%p意义下的逆元
    a[1]=1;
    for(int i=2;i<=n;i++){
        a[i]=1LL*(p-p/i)*a[p%i]%p;
    }
}

 莫比乌斯反演

　　小红书版

int mu[1<<20];
LL n;
bool check[1<<20];
int prime[1<<20];

void Moblus()
{
    memset(check,false,sizeof(check));
    mu[1]=1;
    int tot=0;
    for(int i=2;i<=(1<<20);i++){
        if(!check[i]){
            prime[tot++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<tot;j++){
            if(i*prime[j]>(1<<20)) break;
            check[i*prime[j]]=true;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;break;
            }
            else{
                mu[i*prime[j]]=-mu[i];
            }
        }
    }
}

　　hzw大爷版

int tot;
int mupre[maxn],mu[maxn],pri[maxn];
bool mark[maxn];
void init_mu()
{
    mu[1]=1;
    for(int i=2;i<=maxn;i++){
        if(!mark[i]){mu[i]=-1;pri[++tot]=i;}
        for(int j=1;j<=tot&&i*pri[j]<=maxn;j++){
            mark[i*pri[j]]=1;
            if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
            else mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<=maxn;i++)
        mupre[i]=mupre[i-1]+mu[i];
}

 矩阵快速幂

///使用前要先对r,c赋值
struct mat{
    long long a[30][30];
    int r,c;
    mat operator *(const mat &b)const{
        mat ret;
        for (int i=0;i<r;i++){
            for (int j=0;j<b.c;j++){
                ret.a[i][j]=0;
                for (int k=0;k<c;k++)
                    ret.a[i][j]+=a[i][k]*b.a[k][j],ret.a[i][j]%=MOD;
            }
        }
        ret.r=r;
        ret.c=b.c;
        return ret;
    }
    void init_unit(int x)
    {
        r=c=x;
        for(int i=0;i<r;i++){
            for(int j=0;j<c;j++){
                if(i==j)a[i][j]=1;
                else a[i][j]=0;
            }
        }
    }
}unit;
mat qmod(mat p,LL n){
    unit.init_unit(p.c);
    mat ans=unit;
    while(n){
        if(n&1)ans=p*ans;
        p=p*p;
        n>>=1;
    }
    return ans;
}

求n！中有多少个m的乘积

LL solve(LL n,LL p)
{///返回个数，p为质数
    LL ans=0;
    while(n){
        n=n/p;
        ans+=n;
    }
    return ans;
}