syqwq code template我常用的代码模板，现在主要是数学方面的算法内容，带了一点点的知识点讲解，建议配合o

syqwq code template

[TOC]

ff. TODO

杜教筛
整除分块
广义莫比乌斯反演
异或方程组板子 bitset

0. general

#include <bits/stdc++.h>
#define rep(i, a, b) for (ll i = a; i <= b; i++)
#define per(i, a, b) for (ll i = a; i >= b; i--)
#define pb push_back
#define vc vector
#define fi first
#define se second
#define all(x, n) (x) + 1, (x) + 1 + n
#define NL "\n"
#define NN " "
using namespace std;
typedef long long ll;
typedef vc<ll> vi;
typedef pair<ll, ll> PII;
typedef vc<PII> vpii;
template<class T,class S>
bool chmax(T& x,S y){ return x<y?x=y,1:0; }
template<class T,class S>
bool chmin(T& x,S y){ return x>y?x=y,1:0; }
// #define CF
// #define int long long
//  ================================================

void solve() {
}
//  ================================================
signed main() {
    ios::sync_with_stdio(0);
    cin.tie(0), cout.tie(0);
    cout.setf(ios::fixed), cout.precision(5);
    int T = 1;
#ifdef CF
    cin >> T;
#endif // CF
    while (T--) solve();
    return 0;
}

1. Graph theory

1.0 basic build graph

vpii e[N];
void add(int a,int b,int c){e[a].pb({b,c});}
// ==============================================
int idx=0,h[N],ne[M],to[M],w[M];
void add(int a,int b,int c){
    w[++idx]=c,to[idx]=b,ne[idx]=h[a],ha[a]=idx;
}

1.1 tarjan SCC

有向图中的scc要求当前点能够到达所有点，并且其他点能到当前点，于是利用这个思想，在dfs遍历的时候按照时间戳压栈，那么在上面的点就是他能到的点（如果维护一定性质的话）直到遍历完，如果发现当前点就是这个连通块的最高的点（也就是dfn最小的点），那就把这个scc中的点全部出栈。 scc缩点之后，得到的是拓扑图。

vi e[N];
int dn=0,dfn[N],low[N],stc[N],top=0;
int col[N],cn=0,sz[N];

void scc(int id){
    low[id]=dfn[id]=++dn;
    stc[++top]=id,ins[id]=1;
    for(int it:e[id]){
      if(!dfn[it]){
          scc(it),chmin(low[id],low[it]);
      }else if(ins[it]) chmin(low[id],dfn[it]);
    }
    if(low[id]==dfn[id]){
        cn++;int x;
        do{
          col[x=stc[top--]]=cn,ins[x]=0,sz[cn]++;
        }while(x!=id);
    }
}

1.2 tarjan eDCC

解决方法是记录边的编号。对于 vector，在 push_back 时将当前边的编号一并压入。对于链式前向星，使用成对变换技巧：初始化 cnt = 1，每条边及其反边在链式前向星中存储的编号分别为 2k 和 2k+1 ，将当前边编号异或 1 即得反边编号。时间复杂度 $O(n+m)$ 边双缩点之后，得到的是连通分量作为节点的树。

int dfn[N],low[N],dn=0;
int stc[N],top=0;
int cn=0,col[M];

void form(int id){
    cn++;
    for(int x=0;x!=id;) col[x=stc[top--]]=cn;
}
void dcc(int id,int eid){
    stc[++top]=id,dfn[id]=low[id]=++dn;
    for(auto _:e[id]){
        if(_.se==eid) continue;
        int it=_.fi;
        if(!dfn[it]){
            dcc(it,_.se);
            chmin(low[id],low[it]);
            if(low[it]>low[id]) form(it);
        }else chmin(low[id],dfn[it]);
    }
    if(!eid) form(id);
}

1.3 tarjan vDCC

还不会呜呜呜

1.4 LCA

倍增法，先预处理深度和倍增祖先，然后先跳到一起，再一起跳。时间复杂度 $O(n)$ 预处理， $O(\log n)$ 查询

int dep[N],anc[N][20];
namespace ZX{   // zu xian
  	void bfs(int root){
    	  queue<int> q;
    	  dep[root]=1,q.push(root);
    	  while(q.size()){
	      	  int id=q.front();q.pop();
	      	  for(auto it:e[id]){
	        	    if(dep[it]) continue;
	        	    dep[it]=dep[id]+1;
	        	    anc[it][0]=id;
	        	    rep(i,1,19) anc[it][i]=anc[anc[it][i-1]][i-1];
	        	    q.push(it);
	      	  }
    	  }
    }
    int lca(int x,int y){
        if(dep[x]<dep[y]) swap(x,y);
		    per(i,19,0) if(dep[anc[x][i]]>=dep[y]) x=anc[x][i];
		    if(x==y) return x;
		    per(i,19,0)
		        if(anc[x][i]!=anc[y][i]) x=anc[x][i],y=anc[y][i];
		    return anc[x][0];
	  }
}

2. Math

2.1 Number theory

2.1.1 Eular sieve

线性筛，时间复杂度 $O(n)$

vi primes;
bool st[N];
void ss(int maxn){
    for(int i=2;i<=maxn;i++){
        if(!st[i]) primes.pb(i);
        for(int j:primes){
          if(j>maxn/i) break;
          st[i*j]=1;
          if(i%j==0) break;
      }
    }
}

2.1.2 Eratothenes sieve

埃式筛，时间复杂度 $O(n\log \log n)$ ，修改内层循环可以做区间筛

vi primes;
bitset<N> st;
void ass(int maxn){
    for(int i=2;i<=maxn;i++){
        if(st[i]) continue;
        primes.pb(i);
        for(int j=i*2;j<=maxn;j+=i) st[j]=1;
    }
}

2.1.3 gcd

$(a,b)=(a,b-a)=(b,a-b)$ ，时间复杂度 $O(\log \min\{ a,b \}))$

ll gcd(ll x,ll y){ return y ? x : gcd(y, x%y); }

也可以，cpp内置 __gcd(x,y)

2.2.4 exgcd

方程 $ax+by=\operatorname{gcd}(a,b)$ 的一组特解，通解是 $(x,y)=(x_0+k\cdot \frac{b}{(a,b)},y_0-k\cdot \frac{a}{(a,b)}),k\in \mathbb{Z}$

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

2.2.5 lowbit

ll lowbit(ll x){ return x&(-x); }

2.2.6 popcount

ll count(ll x){ return __builtin_popcountll(x); }
int count(int x){ return __builtin_popcount(x); }

2.2.7 $\phi(n)$

欧拉函数 $\phi(n)=\sum_{i=1}^{n}[\operatorname{gcd}(i,n)=1]$ .

由容斥原理可推出： $n=p_1^{\alpha_1}\cdot p_2^{\alpha_2}\cdots p_j^{\alpha_j}\implies \phi(n)=n\cdot \prod_{i=1}^{j}(1-p_{i}^{-1}) =(p_1-1)(p_2-1)\cdots(p_j-1)\cdot p_1^{\alpha_1-1}\cdot p_2^{\alpha_2-1}\cdots p_j^{\alpha_j-1}$ . 且若 $p$ 为素数，有 $\phi(p)=p-1$ .

可以在 $O(\log n)$ 求单点欧拉函数，可以在 $O(n)$ 递推求 $1$ 到 $n$ 欧拉函数.

// $O(\log n)$ 求单点欧拉函数
ll get_phi(ll x){
	ll phi=1;
	for(ll i=2;i<=x/i;i++){
		if(x%i==0){
			phi*=(i-1);
			x/=i;
			while(x%i==0) phi*=i,x/=i;
		}
	}
	if(x>1) phi*=(x-1);
	return phi;
}
// $O(n)$  递推求 $1$ 到 $n$ 欧拉函数
bool st[N];
vi primes;
ll phi[N];
void get_phi(ll maxn){
	phi[1]=1;
	for(ll i=2;i<=maxn;i++){
		if(!st[i]) primes.pb(i),phi[i]=i-1;
		for(ll j:primes){
			if(j>maxn/i) break;
			st[i*j]=1;
			if(i%j==0) { phi[i*j]=phi[i]*j; break; }
			phi[i*j]=phi[i]*(j-1);
		}
	}
}

2.2.8 $a\cdot a^{-1}\equiv 1(\operatorname{mod} p)$

单点求乘法逆元：可以直接解不定方程， $ab\equiv 1(\operatorname{mod} p) \iff ab=mp+1 \iff ab-mp=1$ ，有解的充要条件是 $(a,p)=1$ . 也可以由欧拉定理 $(n,a)=1 \implies a^{\phi(n)}\equiv 1(\operatorname{mod} n)$ ，特别的，若 $p$ 为素数，则 $a^{p-1}\equiv 1(\operatorname{mod} p)$ .

递推求 $1$ 到 $n$ 的乘法逆元：求 $i$ 的逆元，考虑 $p=\left\lfloor \frac{p}{i} \right\rfloor\cdot i+p \% i\equiv 0(\operatorname{mod} p)$ ，其中注意到 $p\%i<i$ ，于是可以递推。 $\left\lfloor \frac{p}{i} \right\rfloor\cdot i\equiv -(p\%i) \iff i^{-1}\equiv -(p\%i)^{-1}\cdot \left\lfloor \frac{p}{i} \right\rfloor(\operatorname{mod} p)$ .

2.2.9 extended eular theorem

欧拉定理： $(a,m)=1\implies a^{\phi(m)}\equiv 1(\operatorname{mod} m)$ 扩展欧拉定理： $a^{c}\equiv a^{c\%\phi(m)+\phi(m)}(\operatorname{mod} m) \text{if} c\ge \phi(m)$ ，在这里不要求 $(a,m)=1$

2.2.10 linear mod equation system

中国剩余定理和两两相消。考虑方程组

\begin{cases} x \equiv a_1(\operatorname{mod} m_1) \\ x \equiv a_2(\operatorname{mod} m_2) \\ \cdots \\ x \equiv a_k(\operatorname{mod} m_k) \end{cases}

中国剩余定理：如果 $m_1,m_2, \ldots ,m_k$ 两两互素，则有 $x\equiv \sum_{i=1}^{k}M_i' M_i a_i(\operatorname{mod} M)$ ，其中， $M=\prod_{i=1}^{k} m_i$ ， $M_i=M / m_i$ ， $M_i' M_i\equiv 1(\operatorname{mod} m_i)$

2.2.11 Inclusion-Exclusion principle

最简单的形式，便于理解：

\left\vert S-A \cup B \right\vert=\left\vert S \right\vert - \left\vert A \right\vert -\left\vert B \right\vert +\left\vert A \cap B \right\vert

推广：我们将满足某种性质记作 $a_i$ ，不满足某种性质记作 $1-a_i$ ，则有

N((1-a_1)(1-a_2)\cdots (1-a_n))=\sum_{k=0}^{n}(-1)^{k}\sum_{1\le j_1\le \cdots\le j_k\le n}N(a_{j_1}\cdots a_{j_k})

如果选 $k$ 种性质进行容斥是相互等价的，还可以写成：

N((1-a_1)(1-a_2)\cdots (1-a_n))=\sum_{k=0}^{n}(-1)^{k}\binom{n}{k}N(a_1a_2\cdots a_k)

可以考虑 dp ，如果要枚举集合可以使用 dfs，注意记录 $-1$ 的符号

2.2.12 multiplicative function

积性函数：若对于 $f(x)$ ，有 $p \bot q \implies f(pq)=f(p)f(q)$ ，则称 $f$ 为积性函数

如果不要求 $p \bot q$ ，则称为完全积性函数

对于积性函数，可以利用欧拉筛， $O(n)$ 递推他们的值

// 完全积性函数
bool st[N];
vi primes;
ll f[N];
void ss(int maxn){
    for(int i=2;i<=maxn;i++){
        if(!st[i]) primes.pb(i),f[i]=initial_value(i);
        for(int j:primes){
          if(j>maxn/i) break;
          st[i*j]=1,f[i*j]=f[i]*f[j];
          if(i%j==0) break;
      }
    }
}
// 求 p \bot q 的积性函数
// 单点求
ll get_f(ll n){
	ll ans=1;
	for(int i=2;i<=n/i;i++){
		int cnt=0;
		while(n%i==0) cnt++,n/=i;
		ans=ans*f(i,cnt)%mod; // f(p,k)=f(p^k)
	}
	if(n>1) ans=ans*f(n,1)%mod;
	return ans;
}
// 也可以利用最小的质数来递推 1 ~ n
// cnt[] 记录最小质数出现的次数
bool st[N];
vi primes;
ll f[N],cnt[N];
void ss(int maxn){
	f[1]=1;
    for(int i=2;i<=maxn;i++){
        if(!st[i]) primes.pb(i),cnt[i]=1,f[i]=calc_f(i,1);
        for(int j:primes){
            if(j>maxn/i) break;
            st[i*j]=1;
            if(i%j==0){
                cnt[i*j]=cnt[i]+1;
                f[i*j]=f[i]/calc_f(j,cnt[i])*calc_f(j,cnt[i]+1);
                break;
            }
            cnt[i*j]=1;
            f[i*j]=f[i]*calc_f(j,1);
      }
    }
}

2.2.13 Möbius inversion

对于数论函数，常见的两种莫比乌斯反演的两种形式：

对因子反演： $f(n)=\sum_{d \mid n}g(d) \iff g(n)=\sum_{d \mid n} \mu(\frac{n}{d})f(d)=\sum_{d \mid n} \mu(d)f(\frac{n}{d})$ （后面的等号是由于因数的成对出现）
对倍数反演： $f(d)=\sum_{d \mid n , n\le N}g(n) \iff g(d)=\sum_{d \mid n, n\le N}\mu(\frac{n}{d})f(n)$

本质是在整除的意义上划分出的集合上进行容斥。其中 $\mu$ 为莫比乌斯函数：

$\mu(1)=1$
$x=p_1 p_2 \cdots p_k \implies \mu(x)=(-1)^{k}$
$x=p_1^{\alpha_1}p_2^{\alpha_2}\cdots p_k^{\alpha_k} \implies \mu(x)=0$

$\mu$ 为积性函数，利用欧拉筛可以 $O(n)$ 递推求

bool st[N];
vi primes;
int mu[N];
void get_mu(int maxn){
	mu[1]=1;
    for(int i=2;i<=maxn;i++){
        if(!st[i]) primes.pb(i),mu[i]=-1;
        for(int j:primes){
            if(j>maxn/i) break;
            st[i*j]=1;
            if(i%j==0){ mu[i*j]=0; break; }
            mu[i*j]=-mu[i];
      }
    }
}

一般单项不好求，但是如果对于倍数的式子相加，或者因数的式子相加的和（就是考虑一个集合的时候）好求，可以考虑整体求，然后对求和的函数进行反演。

2.2.14 Dirichlet convolution

设 $f\colon \mathbb{N} \to \mathbb{R}$ , $g\colon \mathbb{N} \to \mathbb{R}$ ，则定义他们的狄利克雷卷积为 $(f * g)(n)=\sum_{d \mid n}f(d)g(\frac{n}{d})$ .

若 $f$ 和 $g$ 为积性函数，则他们的卷积也为积性函数，且满足交换律，结合律。

为了方便，我们定义如下函数：

$1(n)=1$ ，在狄利克雷卷积的乘法中与 $\mu$ 互为逆元
$\epsilon(n)=[n=1]$ ，狄利克雷卷积的乘法单位元
$\text{id}(n)=n$

于是，我们有：

$f=g*1 \iff g=f*\mu$ (aka. Möbius inversion)
$\epsilon = \mu * 1 \iff \mu=\epsilon*\mu$ ，证明可以考虑右侧反演单位元或者左侧的话其实就是 $(1+(-1))^{k}$ 的二项式展开
$\text{id}=\phi * 1 \iff \phi = \text{id} * \mu$

一些技巧：

$\sum_{i=1}^{n}i\cdot [i \bot n]=\frac{1}{2}(n\phi(n)+\epsilon(n))$ Proof: 考虑到 $\operatorname{gcd}(i,n)=\operatorname{gcd}(n-i,n)$ ，所以他们是成对出现的，一共的对数就是 $\phi(n)$ ，每一对贡献的和是 $n$ ，特判 $n=1$ 的情况。

2.2.15 matrix multiplication

可以用来加速线性 dp 的递推。数据结构中，将维护值扩展成成维护矩阵。

$A \in M_{n \times t}(\mathbb{R})$ , $B\in M_{t\times m}(\mathbb{R})$ ,那么 $AB[i,j]=\sum_{k=1}^{t}A[i,k]\cdot B[k,m]$ .

const int N=105;
const ll mod=1e9+7;
template<class T>
struct mat{
	int n,m;    // 0~n-1, 0~m-1
	T a[N][N];
	mat(int n,int m) { this->n=n,this->m=m; memset(a, 0, sizeof a); }
    T* operator[](int x) { return a[x]; }
    const T* operator[](int x) const { return a[x]; }
    mat operator-(const mat &rhs) const {
        mat<T> res(n,m);
        for(int i=0;i<n;i++) for(int j=0;j<m;j++)
            res[i][j]=(a[i][j]+mod-rhs[i][j])%mod;
        return res;
    }
    mat operator+(const mat &rhs) const {
        mat<T> res(n,m);
        for(int i=0;i<n;i++) for(int j=0;j<m;j++)
            res[i][j]=(a[i][j]+rhs[i][j])%mod;
        return res;
    }
    mat operator*(const mat &rhs) const {
    	// assert(m==rhs.n);
        mat<T> res(n,rhs.m);
		for(int i=0;i<n;i++) for(int j=0;j<rhs.m;j++) for(int k=0;k<m;k++)
			res[i][j]+=a[i][k]*rhs[k][j], res[i][j]%=mod;
        return res;
    }
    mat qmi(int k) {
    	// assert(n==m);
        mat<T> res(n, n), t=*this;
        for(int i=0;i<n;i++) res[i][i]=1;
        for(;k;t=t*t,k>>=1) if(k&1) res=res*t;
        return res;
    }
};

2.2.16 Gauss elimination

线性方程组 $Ax=b$ 的解。设 $A$ 的增广矩阵记作 $\tilde{A}=[A \mid b]$ ，将解集记作 $S=\{ x \mid Ax=b \}$ ，则有

$\left\vert S \right\vert =0 \iff \operatorname{rank}A<\operatorname{rank}\tilde{A}$
$\left\vert S \right\vert =1\iff \operatorname{rank}A=\operatorname{rank}\tilde{A}$
$\left\vert S \right\vert = \aleph_1 \iff \operatorname{rank}A>\operatorname{rank}\tilde{A}$

constexpr int N=105;
constexpr double eps=1e-6;
template<class T=double>
struct gmat{
	int n,m;	// 0~n-1, 0~m-1
	T a[N][N];
	gmat(int n,int m) { this->n=n,this->m=m; memset(a, 0, sizeof a); }
    T* operator[](int x) { return a[x]; }
    const T* operator[](int x) const { return a[x]; }
    // 0: no solution || 1: one solution || -1: inf solution
	int gauss(){
		int c,r;
		for(c=0,r=0;c<n;c++){
			int t=r;
			for(int i=r+1;i<n;i++) if(fabs(a[t][c])<fabs(a[i][c])) t=i;
			if(fabs(a[t][c])<eps) continue;
			for(int i=c;i<m;i++) swap(a[t][i],a[r][i]);
			for(int i=m-1;i>=c;i--) a[r][i]/=a[r][c];
			for(int i=r+1;i<n;i++) 
				if(fabs(a[i][c])>eps) 
					for(int j=m-1;j>=c;j--) a[i][j]-=a[i][c]*a[r][j];
			r++;
		}
		if(r<n){
			for(int i=r;i<n;i++) if(fabs(a[i][m-1])>eps) return 0;
			return -1;
		}
		for(int i=n-1;i>=0;i--) for(int j=0;j<i;j++) a[j][m-1]-=a[i][m-1]*a[j][i];
		return 1;
	}
};

求解异或方程组在时间复杂度要求不严格的情况下，可以使用 bit 结构体：

// 记得修改 gauss 中关于精度的判断
struct bit {
    bool val;
    bit(bool val=0):val(val) {}
    bit operator+(const bit& x) const { return bit(val^x.val); }
    bit operator-(const bit& x) const { return bit(val^x.val); }
    void operator-=(const bit& x) { val^=x.val; }
    bit operator*(const bit& x) const { return bit(val&x.val); }
    bit operator/(const bit& x) const { return bit(val); }
    void operator/=(const bit& x) { val=val/x.val; }
    bool operator<(const bit& x) const { return !val&&x.val; }
    bool operator==(const int x) const { return val==x; }
    bool operator!=(const int x) const { return !(*this==x); }
    operator bool() const { return val; }
    friend istream& operator>>(istream& in,bit& b){ return in>>b.val,in; }
    friend ostream& operator<<(ostream& out,const bit& b){ return out<<b.val,out; }
};

或者，使用 bitset 优化

struct bmat{
    int n,m;

}

2.2.27 linear basis

布尔域 $\mathbb{Z}_{2}$ 和异或（加法），逻辑与（数乘）构成线性空间，其中的极大线性无关组称为线性基。 $p[i]$ 表示最高位是 $i$ 位的基向量。

typedef unsigned long long ull;
constexpr int N=105;
constexpr int B=50;
int rk=0;
ull p[N];
void insert(ull x){
	for(int i=B;i>=0;i--){
		if(!(x>>i&1)) continue;
		if(!p[i]) return p[i]=x,void();
		x^=p[i];
	}
}

syqwq code template