组合数
递推
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; }