我的算法竞赛模板

本页是一个汇总页。

我的算法竞赛模板

本页是一个汇总页。

{% for post in site.tags[‘算法竞赛模板’] %}

{{ post.title }}

{{post.content}}

{% endfor %}

数据结构

以下数据结构均采用 ll 作为值类型,应用时根据需求调整。

typedef long long ll;
const ll INF = 1e9;  //表示(值)正无穷,且两个正无穷相加不会溢出
const int NPOS = -1; //表示(下标)不存在

离散化

在 vector 基础上的离散化,使用 push_back()向其中插值,init()排序并离散化,ask 查询离散化之后的值,at/[]运算符查离散前的值。

struct Ranker : vector<ll>
{
	void init() { sort(begin(), end()), resize(unique(begin(), end()) - begin()); }
	int ask(ll x) const { return lower_bound(begin(), end(), x) - begin(); }
};

并查集

struct UnionfindSet : vector<int>
{
	UnionfindSet(int n) : vector<int>(n)
	{
		for (int i = 0; i < n; ++i)
			at(i) = i;
	}
	void merge(int u, int w)
	{
		if (w = ask(w), u = ask(u), w != u)
			at(w) = u;
	}
	int ask(int u) { return at(u) != u ? at(u) = ask(at(u)) : u; }
};

单调队列和单调栈

使用示例

每次取队尾就是单调栈,取队头就是单调队列。

typedef pair<ll, int> pli;
struct Monotone : deque<pli>
{
	void push(const pli &p, int k)
	{
		while (!empty() && back().first >= p.first)
			pop_back();
		for (push_back(p); p.second - front().second >= k;)
			pop_front();
	}
};

ST 表

使用示例

$O(n\log n)$预处理,$O(1)$求静态区间最小值。

/*
//可选优化
#define log2(n) LOG2[n]
struct Log : vector<ll>
{
	Log(int N, ll E) : vector<ll>(N, -1)
	{
		for (int i = 1; i < N; ++i)
			at(i) = at(i / E) + 1;
	}
} LOG2(N, 2);
*/
struct SparseTable
{
	vector<vector<ll>> f;
	SparseTable(const vector<ll> &a) : f(log2(a.size()) + 1, a)
	{
		for (int k = 0; k + 1 < f.size(); ++k)
			for (int i = 0; i + (1 << k) < a.size(); ++i)
				f[k + 1][i] = min(f[k][i], f[k][i + (1 << k)]);
	}
	ll ask(int l, int r)
	{
		int k = log2(r - l + 1);
		return min(f[k][l], f[k][r + 1 - (1 << k)]);
	}
};

树状数组

模板中 Base 是对应的基础版本,支持单点修改区间查询。

一维

struct Fenwick
{
	struct BaseFenwick
	{
		vector<ll> v;
		BaseFenwick(int n) : v(n, 0) {}
		void add(int x, ll w)
		{
			for (; x < v.size(); x += x & -x)
				v[x] += w;
		}
		ll ask(int x)
		{
			ll ans = 0;
			for (; x; x -= x & -x)
				ans += v[x];
			return ans;
		}
	};
	pair<BaseFenwick, BaseFenwick> p;
	Fenwick(int n) : p(n, n) {}
	void add(int x, ll w) { p.first.add(x, w), p.second.add(x, x * w); }
	void add(int l, int r, ll w) { add(l, w), add(r + 1, -w); }
	ll ask(int x) { return (x + 1) * p.first.ask(x) - p.second.ask(x); }
	ll ask(int l, int r) { return ask(r) - ask(l - 1); }
};

二维

高维的数据结构只要每一维维护低一维的数据(树套树)即可。其余数据结构亦同理。

struct Fenwick2
{
	struct BaseFenwick2
	{
		vector<Fenwick> v;
		BaseFenwick2(int r, int c) : v(r, c) {}
		void add(int x, int b, int t, ll w)
		{
			for (; x < v.size(); x += x & -x)
				v[x].add(b, t, w);
		}
		ll ask(int x, int b, int t)
		{
			ll ans = 0;
			for (; x; x -= x & -x)
				ans += v[x].ask(b, t);
			return ans;
		}
	};
	pair<BaseFenwick2, BaseFenwick2> p;
	Fenwick2(int r, int c) : p(BaseFenwick2(r, c), BaseFenwick2(r, c)) {}
	void add(int x, int b, int t, ll w) { p.first.add(x, b, t, w), p.second.add(x, b, t, x * w); }
	void add(int l, int b, int r, int t, ll w) { add(l, b, t, w), add(r + 1, b, t, -w); } //(l,b)~(r,t)
	ll ask(int x, int b, int t) { return (x + 1) * p.first.ask(x, b, t) - p.second.ask(x, b, t); }
	ll ask(int l, int b, int r, int t) { return ask(r, b, t) - ask(l - 1, b, t); }
};

动态开点线段树

使用示例,支持区间线性变换、区间查询(最大值最小值区间和)。

这样写改可持久化也很方便,只要改down函数为每次都新建节点即可,示例

struct SegmentTree
{
	struct Seg
	{
		int l, r;
		ll min, max, sum;
		void upd(ll mul, ll add) { min = min * mul + add, max = max * mul + add, sum = sum * mul + add * (r - l + 1); }
		friend Seg operator+(const Seg &lc, const Seg &rc) { return {lc.l, rc.r, std::min(lc.min, rc.min), std::max(lc.max, rc.max), lc.sum + rc.sum}; }
	};
	struct Node : Seg
	{
		int lc, rc;
		ll mul, add;
	};
	vector<Node> v;
	SegmentTree(int l, int r) { build(l, r); }
	void build(int l, int r)
	{
		int rt = v.size();
		v.push_back({});
		v[rt].Seg::operator=({l, r, 0, 0, 0});
		v[rt].lc = v[rt].rc = NPOS;
		v[rt].mul = 1, v[rt].add = 0;
		//if (l < r) //动态开点的时候注释掉本行和下一行
		//down(rt), v[rt].Seg::operator=(v[v[rt].lc] + v[v[rt].rc]);
	}
	void down(int rt)
	{
		int m = v[rt].l + (v[rt].r - v[rt].l >> 1);
		if (v[rt].lc == NPOS)
			v[rt].lc = v.size(), build(v[rt].l, m);
		//else //非可持久化的时候注释掉本行和下一行
		//v.push_back(v[v[rt].lc]), v[rt].lc = v.size() - 1;
		if (v[rt].rc == NPOS)
			v[rt].rc = v.size(), build(m + 1, v[rt].r);
		//else //非可持久化的时候注释掉本行和下一行
		//v.push_back(v[v[rt].rc]), v[rt].rc = v.size() - 1;
		upd(v[v[rt].lc].l, v[v[rt].lc].r, v[rt].mul, v[rt].add, v[rt].lc);
		upd(v[v[rt].rc].l, v[v[rt].rc].r, v[rt].mul, v[rt].add, v[rt].rc);
		v[rt].mul = 1, v[rt].add = 0;
	}
	void upd(int l, int r, ll mul, ll add, int rt = 0)
	{
		if (l <= v[rt].l && v[rt].r <= r)
			return v[rt].mul *= mul, v[rt].add = v[rt].add * mul + add, v[rt].upd(mul, add);
		down(rt);
		if (r <= v[v[rt].lc].r)
			upd(l, r, mul, add, v[rt].lc);
		else if (l >= v[v[rt].rc].l)
			upd(l, r, mul, add, v[rt].rc);
		else
			upd(l, v[v[rt].lc].r, mul, add, v[rt].lc), upd(v[v[rt].rc].l, r, mul, add, v[rt].rc);
		v[rt].Seg::operator=(v[v[rt].lc] + v[v[rt].rc]);
	}
	Seg ask(int l, int r, int rt = 0)
	{
		if (l <= v[rt].l && v[rt].r <= r)
			return v[rt];
		down(rt);
		if (r <= v[v[rt].lc].r)
			return ask(l, r, v[rt].lc);
		if (l >= v[v[rt].rc].l)
			return ask(l, r, v[rt].rc);
		return ask(l, v[v[rt].lc].r, v[rt].lc) + ask(v[v[rt].rc].l, r, v[rt].rc);
	}
};

无旋 Treap

按子树大小分裂

使用示例

struct FhqTreap
{
	struct Node
	{
		int ch[2], siz, rev;
		ll key, val, min, add;
		void upd(ll v, int r)
		{
			val += v, min += v, add += v;
			if (r)
				rev ^= 1, swap(ch[0], ch[1]);
		}
	};
	vector<Node> v;
	int root;
	FhqTreap() : v(1), root(0) {}
	void down(int k)
	{
		if (!k)
			return;
		for (int i = 0, *ch = v[k].ch; i < 2; ++i)
			if (ch[i])
				v[ch[i]].upd(v[k].add, v[k].rev);
		v[k].add = v[k].rev = 0;
	}
	void up(int k)
	{
		if (!k)
			return;
		v[k].siz = 1, v[k].min = v[k].val;
		for (int i = 0, *ch = v[k].ch; i < 2; ++i)
			if (ch[i])
				v[k].siz += v[ch[i]].siz, v[k].min = min(v[k].min, v[ch[i]].min);
	}
	int merge(int a, int b)
	{
		if (!a || !b)
			return a + b;
		if (v[a].key < v[b].key)
			return down(a), v[a].ch[1] = merge(v[a].ch[1], b), up(a), a;
		return down(b), v[b].ch[0] = merge(a, v[b].ch[0]), up(b), b;
	}
	void split(int a, int s, int &l, int &r)
	{
		if (!s)
			l = 0, r = a;
		else if (v[v[a].ch[0]].siz < s)
			down(a), split(v[a].ch[1], s - v[v[a].ch[0]].siz - 1, v[a].ch[1], r), up(l = a);
		else
			down(a), split(v[a].ch[0], s, l, v[a].ch[0]), up(r = a);
	}
	void push_back(ll d) { v.push_back(Node{{0, 0}, 1, 0, rand(), d, d, d}), root = merge(root, v.size() - 1); }
	void insert(int x, ll d)
	{
		v.push_back(Node{{0, 0}, 1, 0, rand(), d, d, d});
		int a, b, c;
		split(root, x - 1, a, b), root = merge(merge(a, v.size() - 1), b);
	}
	void erase(int x)
	{
		int a, b, c;
		split(root, x, a, b), split(a, x - 1, a, c), root = merge(a, b);
	}
	Node ask(int l, int r)
	{
		int a, b, c;
		split(root, r, b, c), split(b, l - 1, a, b);
		Node ret = v[b];
		return root = merge(merge(a, b), c), ret;
	}
	void upd(int l, int r, ll add, int rev)
	{
		int a, b, c;
		split(root, r, b, c), split(b, l - 1, a, b), v[b].upd(add, rev), root = merge(merge(a, b), c);
	}
	void revolve(int l, int r, int d)
	{
		int a, b, c, e = r - l + 1;
		split(root, r, b, c), split(b, l - 1, a, b), split(b, (e - d % e) % e, b, e), root = merge(merge(a, merge(e, b)), c);
	}
};

按值大小分裂

使用示例,即排序树。

struct FhqTreap
{
	struct Node
	{
		int ch[2], siz;
		ll key, val;
	};
	vector<Node> v;
	int root;
	FhqTreap() : v(1), root(0) {}
	void up(int k) { v[k].siz = v[v[k].ch[0]].siz + v[v[k].ch[1]].siz + 1; }
	int merge(int a, int b)
	{
		if (!a || !b)
			return a + b;
		if (v[a].key < v[b].key)
			return v[a].ch[1] = merge(v[a].ch[1], b), up(a), a;
		return v[b].ch[0] = merge(a, v[b].ch[0]), up(b), b;
	}
	void splitVal(int a, ll w, int &l, int &r) //按值将树划分,使得左子树上的值恰小于w
	{
		if (!a)
			l = r = 0;
		else if (v[a].val > w)
			splitVal(v[a].ch[0], w, l, v[a].ch[0]), up(r = a);
		else
			splitVal(v[a].ch[1], w, v[a].ch[1], r), up(l = a);
	}
	void insert(ll x)
	{
		int a, b;
		v.push_back(Node{{0, 0}, 1, rand(), x}), splitVal(root, x, a, b), root = merge(merge(a, v.size() - 1), b);
	}
	void erase(ll x)
	{
		int a, b, c;
		splitVal(root, x, a, b), splitVal(a, x - 1, a, c), root = merge(merge(a, merge(v[c].ch[0], v[c].ch[1])), b);
	}
	ll kth(int k)
	{
		for (int u = root, ls;;)
		{
			if (ls = v[v[u].ch[0]].siz, ls + 1 == k)
				return v[u].val;
			if (ls < k)
				k -= ls + 1, u = v[u].ch[1];
			else
				u = v[u].ch[0];
		}
	}
	int lower_bound(ll x) { return upper_bound(x - 1); }
	int upper_bound(ll x)
	{
		int a, b, ret;
		return splitVal(root, x, a, b), ret = v[a].siz + 1, root = merge(a, b), ret;
	}
};

莫队

使用示例

struct Mo
{
	struct Query
	{
		int l, r, id;
		bool operator<(const Query &n) const
		{
			return l / BS != n.l / BS ? l < n.l : r < n.r;
		}
	};
	vector<Query> q;
	int L, R;
	void query(int l, int r) { q.push_back(Query{l, r, q.size()}); }
	void rev(int x) {}
	void cal(int id) {}
	void ask()
	{
		L = 0, R = -1;
		sort(q.begin(), q.end());
		for (int i = 0; i < q.size(); ++i)
		{
			while (L < q[i].l)
				rev(L++);
			while (L > q[i].l)
				rev(--L);
			while (R < q[i].r)
				rev(++R);
			while (R > q[i].r)
				rev(R--);
			cal(q[i].id);
		}
	}
};

带修莫队

使用示例

struct Mo
{
	struct Update
	{
		int pos, NEW, OLD;
	};
	struct Query
	{
		int t, l, r, id;
		bool operator<(const Query &n) const
		{
			return l / BS != n.l / BS ? l < n.l : r / BS != n.r / BS ? r < n.r : t < n.t;
		}
	};
	vector<Update> cq;
	vector<Query> q;
	int T, L, R;
	Mo() : cq(1) {}
	void query(int x, int y) { q.push_back(Query{cq.size() - 1, x, y, q.size()}); }
	void update(int x, int y) { cq.push_back(Update{x, y, t[x]}), t[x] = y; }
	void set(int x, int d)
	{
		if (vis[x])
			return rev(x), a[x] = d, rev(x);
		a[x] = d;
	}
	void rev(int x) {}
	void cal(int id) {}
	void ask()
	{
		T = L = 0, R = -1;
		sort(q.begin(), q.end());
		for (int i = 0; i < q.size(); ++i)
		{
			while (T < q[i].t)
				++T, set(cq[T].pos, cq[T].NEW);
			while (T > q[i].t)
				set(cq[T].pos, cq[T].OLD), --T;
			while (L < q[i].l)
				rev(L++);
			while (L > q[i].l)
				rev(--L);
			while (R < q[i].r)
				rev(++R);
			while (R > q[i].r)
				rev(R--);
			cal(q[i].id);
		}
	}
};

树上莫队

使用示例

按照欧拉序分块,使用 Tarjan 在生成欧拉序的同时预处理所有询问的 lca,预处理时间复杂度$O(n+q)$。 h 为查询图,即如果有一个询问(u,v),即在 h 上连$u\to v,v\to u$。多个询问边有序插入 h。

struct TreeMo : Graph
{
	struct Query
	{
		int l, r, lca, id;
		bool operator<(const Query &b) const
		{
			return l / BS != b.l / BS ? l < b.l : r < b.r;
		}
	};
	vector<Query> q;
	vector<int> dfp, dfi, dfo;
	UnionFindSet ufs;
	Graph h;
	int L, R;
	TreeMo(int n) : Graph(n), h(n), dfp(n * 2 + 1), dfi(n), dfo(n), ufs(n) {}
	void query(int x, int y)
	{
		h.add(Edge{x, y}), h.add(Edge{y, x});
		q.push_back(Query{0, 0, 0, q.size()});
	}
	void rev(int x) {}
	void cal(int id) {}
	void dfs(int u, int &cnt)
	{
		dfp[dfi[u] = ++cnt] = u;
		for (int i = 0, k, to; i < v[u].a.size(); ++i)
			if (k = v[u].a[i], to = e[k].second, !dfi[to])
				dfs(to, cnt), ufs.merge(u, to);
		dfp[dfo[u] = ++cnt] = u;
		for (int i = 0, k, to, id; i < h.v[u].a.size(); ++i)
			if (k = h.v[u].a[i], id = k / 2, to = h.e[k].second, dfo[to])
			{
				q[id].lca = ufs.fa(to);
				q[id].l = q[id].lca != u ? dfo[u] : dfi[u];
				q[id].r = dfi[to];
			}
	}
	void ask(int root = 1)
	{
		dfs(root, BS = 0), BS = sqrt(BS);
		sort(q.begin(), q.end());
		L = 0, R = -1;
		for (int i = 0; i < q.size(); ++i)
		{
			while (L < q[i].l)
				rev(dfp[L++]);
			while (L > q[i].l)
				rev(dfp[--L]);
			while (R < q[i].r)
				rev(dfp[++R]);
			while (R > q[i].r)
				rev(dfp[R--]);
			if (q[i].lca != dfp[L])
				rev(q[i].lca);
			cal(q[i].id);
			if (q[i].lca != dfp[L])
				rev(q[i].lca);
		}
	}
};

树上带修莫队

使用示例

struct CapitalTreeMo : Graph
{
	struct Update
	{
		int pos, NEW, OLD;
	};
	struct Query
	{
		int t, l, r, lca, id;
		bool operator<(const Query &b) const
		{
			return l / BS != b.l / BS ? l < b.l : r / BS != b.r / BS ? r < b.r : t < b.t; //在BZOJ4129上去掉r/BS还快100ms?
		}
	};
	vector<Update> cq;
	vector<Query> q;
	vector<int> dfp, dfi, dfo;
	UnionFindSet ufs;
	Graph h;
	int T, L, R;
	CapitalTreeMo(int n) : cq(1), Graph(n), h(n), dfp(n * 2 + 1), dfi(n), dfo(n), ufs(n) {}
	void query(int x, int y)
	{
		h.add(Edge{x, y}), h.add(Edge{y, x});
		q.push_back(Query{cq.size() - 1, 0, 0, 0, q.size()});
	}
	void update(int x, int y)
	{
		cq.push_back(Update{x, y, t[x]}), t[x] = y;
	}
	void dfs(int u, int &cnt)
	{
		dfp[dfi[u] = ++cnt] = u;
		for (int i = 0, k, to; i < v[u].a.size(); ++i)
			if (k = v[u].a[i], to = e[k].second, !dfi[to])
				dfs(to, cnt), ufs.merge(u, to);
		dfp[dfo[u] = ++cnt] = u;
		for (int i = 0, k, to, id; i < h.v[u].a.size(); ++i)
			if (k = h.v[u].a[i], id = k / 2, to = h.e[k].second, dfo[to])
			{
				q[id].lca = ufs.fa(to);
				q[id].l = q[id].lca != u ? dfo[u] : dfi[u];
				q[id].r = dfi[to];
			}
	}
	void set(int u, int d)
	{
		if (vis[u])
			return rev(u), a[u] = d, rev(u);
		a[u] = d;
	}
	void rev(int u) {}
	void cal(int id) {}
	void ask(int root = 1)
	{
		dfs(root, BS = 0), BS = sqrt(BS);
		sort(q.begin(), q.end());
		T = L = 0, R = -1;
		for (int i = 0; i < q.size(); ++i)
		{
			while (T < q[i].t)
				++T, set(cq[T].pos, cq[T].NEW);
			while (T > q[i].t)
				set(cq[T].pos, cq[T].OLD), --T;
			while (L < q[i].l)
				rev(dfp[L++]);
			while (L > q[i].l)
				rev(dfp[--L]);
			while (R < q[i].r)
				rev(dfp[++R]);
			while (R > q[i].r)
				rev(dfp[R--]);
			if (q[i].lca != dfp[L])
				rev(q[i].lca);
			cal(q[i].id);
			if (q[i].lca != dfp[L])
				rev(q[i].lca);
		}
	}
};

字符串/模式匹配

HashString

使用示例,如果要修改模数或者直接使用unsigned long long的自然溢出的话直接修改 Mod 即可。

使用unsigned long long的自然溢出快了 5 倍,但是容易被卡。

struct HashString : Mod
{
	vector<ll> f, p;
	HashString(const string &s, ll M = 1e9 + 7, ll P = 131) : Mod(M), f(s.size() + 1), p(s.size() + 1, 1)
	{
		for (int i = 0; i < s.size(); ++i)
		{
			f[i + 1] = add(mul(f[i], P), s[i]);
			p[i + 1] = mul(p[i], P);
		}
	}
	ll ask(int pos, int len) { return add(f[pos + len], -mul(f[pos], p[len])); } //从pos位置开始的长度为len的子串的hash值
};

KMP

struct KMP
{
	vector<int> next;
	string pattern;
	KMP(const string &pattern) : pattern(pattern), next(pattern.size() + 1, -1)
	{
		for (int i = 0, j = -1; i < pattern.size();)
		{
			if (j == -1 || pattern[i] == pattern[j])
				next[++i] = ++j;
			else
				j = next[j];
		}
	}
	int find_in(const string &text)
	{
		for (int i = 0, j = 0;;)
		{
			if (j == pattern.size())
				return i - j; //不return可得到所有匹配地址
			if (i == text.size())
				return -1;
			if (j == -1 || text[i] == pattern[j])
				++i, ++j;
			else
				j = next[j];
		}
	}
};

AC 自动机

struct AhoCorasick
{
	struct Node
	{
		int ch[26], val, f, last;
		int &to(char c)
		{
			return ch[c - 'a'];
		} //如果不确定c的范围,使用map
	};
	vector<Node> v;
	AhoCorasick() : v(1) {}
	void getFail()
	{
		for (deque<int> q(1, v[0].last = v[0].f = 0); !q.empty(); q.pop_front())
			for (char c = 'a'; c <= 'z'; ++c)
			{
				int r = q.front(), u = v[r].to(c), w = v[r].f;
				if (!r && u)
				{
					q.push_back(u);
					v[u].f = v[u].last = 0;
					continue;
				}
				if (!u)
				{
					v[r].to(c) = v[w].to(c);
					continue;
				}
				q.push_back(u);
				while (w && !v[w].to(c))
					w = v[w].f;
				v[u].f = v[w].to(c);
				v[u].last = v[v[u].f].val ? v[u].f : v[v[u].f].last;
			}
	}
	void add(const string &s, int val, int u = 0)
	{
		for (int i = 0; i < s.size(); u = v[u].to(s[i++]))
			if (!v[u].to(s[i]))
			{
				v[u].to(s[i]) = v.size();
				v.push_back(Node());
			}
		v[u].val = val;
	}
	bool find_in(const string &s, int u = 0) //调用需要调用`getFail()`生成失配函数。
	{
		for (int i = 0; i < s.size(); ++i)
			if (u = v[u].to(s[i]),
				v[u].val || v[u].last)
				return 1;
		return 0;
	}
};

暴力回文

使用示例

时间复杂度$O(n^2)$,常数低,但会被ababababa这样的数据卡。

int palindrome(const char *s)
{
	int ans = 0;
	for (int i = 0, b, e; s[i]; ++i)
	{
		for (b = i; s[i] == s[i + 1];)
			++i;
		for (e = i + 1; b && s[b - 1] == s[e];)
			--b, ++e;
		if (ans < e - b)
			ans = e - b; //此时[b,e)为最大回文区间
	}
	return ans;
}

线性回文

使用示例

对于一个位置 i,[i−f[i]+1,i+f[i]−1]是最长的以 i 为中心的奇回文串,g[i]−i 是最长的以 i 为开头的回文串长度。

struct Manacher
{
	vector<int> t, f, g;
	Manacher(const string &s) : t(s.size() + 1 << 1, 0), f(t), g(t) //t初始值为s中没有出现过的值,g开始为0
	{
		for (int i = 0; i < s.size(); ++i)
			t[i + 1 << 1] = s[i];
		for (int i = 1, p = 0, m = 0; i < t.size(); ++i)
		{
			for (f[i] = i < m ? min(f[2 * p - i], m - i) : 1;
				 0 < i - f[i] && i + f[i] < t.size() &&
				 t[i - f[i]] == t[i + f[i]];)
				++f[i];
			if (m < i + f[i])
				m = i + f[p = i];
		}
		for (int i = 2; i < t.size(); ++i)
			if (g[i - f[i] + 1] < i + 1)
				g[i - f[i] + 1] = i + 1;
		for (int i = 1; i < t.size(); ++i)
			if (g[i] < g[i - 1])
				g[i] = g[i - 1];
	}
	int ask(int l, int r) //多次询问可做一个ST表
	{
		int ans = 0;
		for (int i = l + 1 << 1, e = r + 1 << 1; i <= e; i += 2)
			if (ans < g[i] - i)
				ans = g[i] - i;
		return ans;
	}
};

后缀数组

使用示例

m:字符集大小。

s:字符串,其中最后一位为加入的 0。

sa[i]:字典序第 i 小的是哪个后缀。

rk[i]:后缀 i 的排名。

h[i]:lcp(sa[i],sa[i−1])。

struct SufArr
{
	vector<int> sa, rk, h;
	SufArr(const vector<int> &s, int m) : sa(s.size(), 0), rk(s), h(s.size(), 0)
	{
		vector<int> cnt(s.size() + m, 0);
		for (int i = 0; i < s.size(); ++i)
			++cnt[rk[i]];
		for (int i = 1; i < m; ++i)
			cnt[i] += cnt[i - 1];
		for (int i = 0; i < s.size(); ++i)
			sa[--cnt[rk[i]]] = i;
		for (int k = 1, j = 0; k <= s.size() && j < s.size() - 1; k <<= 1)
		{
			for (int i = 0; i < s.size(); ++i)
			{
				if (j = sa[i] - k, j < 0)
					j += s.size();
				h[cnt[rk[j]]++] = j;
			}
			cnt[0] = sa[h[0]] = j = 0;
			for (int i = 1; i < s.size(); ++i)
			{
				if (rk[h[i]] != rk[h[i - 1]] || rk[h[i] + k] != rk[h[i - 1] + k])
					cnt[++j] = i;
				sa[h[i]] = j;
			}
			swap(rk, sa), swap(sa, h);
		}
		for (int i = 0, k = 0, j = rk[0]; i < s.size() - 1; ++i, ++k)
			for (; ~k && s[i] != s[sa[j - 1] + k]; j = rk[sa[j] + 1], --k)
				h[j] = k;
	}
};

图论

这里用类似邻接表的方法存图。有的算法可能需要邻接矩阵,详见线性代数部分。

struct Graph
{
	struct Vertex
	{
		vector<int> o, i;		//相关出边和入边编号
		int siz, dep, top, dfn; //树链剖分中使用,依次代表子树节点数、深度、所在链的顶端节点、dfs序
	};
	struct Edge : pair<int, int>
	{
		ll len, cap; //边长、容量,图论算法使用
	};
	vector<Vertex> v; //点集
	vector<Edge> e;   //边集
	Graph(int n) : v(n) {}
	void add(const Edge &ed)
	{
		if (ed.first == ed.second)
			return; //如果有需要请拆点
		v[ed.first].o.push_back(e.size());
		v[ed.second].i.push_back(e.size());
		e.push_back(ed);
	}
	int ch(int u, int i = 0) { return e[v[u].o[i]].second; } //u的第i个孩子节点
	int fa(int u, int i = 0) { return e[v[u].i[i]].first; } //u的第i个父节点
};

最短路

Dijkstra 算法

使用示例,适用于边权为正的情况(无论有向图还是无向图),用于求单源最短路。

直接给出其优先队列优化的版本。另外,由于priority_queue并不提供修改优先级的操作,为避免重复扩展,这里选择将新元素直接插入队列并在运行时判断该点是否被处理,并不影响结果的正确性。

struct Dijkstra : Graph
{
	vector<ll> d;
	vector<int> p;
	Dijkstra(int n) : Graph(n) {}
	void ask(int s)
	{
		d.assign(v.size(), INF);
		p.assign(v.size(), e.size());
		priority_queue<pair<ll, int>> q;
		for (q.push(make_pair(d[s] = 0, s)); !q.empty();)
		{
			ll dis = -q.top().first;
			int u = q.top().second;
			if (q.pop(), d[u] < dis)
				continue;
			for (int i = 0, k, to; i != v[u].o.size(); ++i)
				if (k = v[u].o[i], to = e[k].second,
					d[to] > d[u] + e[k].len)
				{
					d[to] = d[u] + e[k].len, p[to] = k;
					q.push(make_pair(-d[to], to));
				}
		}
	}
};

BellmanFord 算法

使用示例,直接给出其队列优化、国内称之为 SPFA 算法的版本。较之 Dijkstra 算法,此算法不够快速稳定但是可以允许负边存在,当 s 到达负权回路时会直接返回 0。稀疏图上性能优秀。

struct BellmanFord : Graph
{
	vector<ll> d;
	vector<int> p;
	BellmanFord(int n) : Graph(n) {}
	bool ask(int s)
	{
		d.assign(v.size(), INF);
		p.assign(v.size(), e.size());
		vector<int> cnt(v.size(), 0), flag(v.size(), d[s] = 0);
		for (deque<int> q(cnt[s] = flag[s] = 1, s); !q.empty(); q.pop_front())
			for (int u = q.front(), i = flag[u] = 0, k, to; i < v[u].o.size(); ++i)
				if (k = v[u].o[i], to = e[k].second,
					d[to] > d[u] + e[k].len)
				{
					d[to] = d[u] + e[k].len, p[to] = k;
					if (!flag[to])
					{
						if (v.size() == ++cnt[to])
							return 0;
						flag[to] = 1, q.push_back(to);
					}
				}
		return 1;
	}
};

差分约束系统

按如下方式建图、跑 SPFA:

对每个不等式$x_i−x_j\leq d$,从$j$向$i$连一条边,边长为$d$。

若不等号的方向相反,即$x_i−x_j\geq d$,则在不等式两边同时乘以$-1$,变成$x_j−x_i\leq -d$,即从$i$到$j$连一条边,边长为$d$。

Floyed 求多源最短路

不连通的权置 INF。

void floyed(Mat &a, int n) // 其实一般手打即可,没必要套这个鬼模板…
{
	for (int k = 0; k < n; ++k)
		for (int i = 0; i < n; ++i)
			for (int j = 0; j < n; ++j)
				if (a[i][j] > a[i][k] + a[k][j])
					a[i][j] = a[i][k] + a[k][j];
}

Astar 求 k 短路

使用示例,在这个比较坑的例子中需要在调用前补一句if(s==t)++k;k 下标从 0 开始,即最短路==第 0 短路

朴素的想法是使用priority_queue从原点出发向外探索,当取出终点 t 第 k 次时就得到第 k 短路,类似 bfs 的思想,缺陷是越往后状态数越多。

我们在反向图上从$t\to s$跑 Astar 算法,通过优先展开到 s 近的状态,使搜索方向靠近答案,而不是一层一层全都展开,估价函数$f\approx g+h$,f 是估计的 s 到 t 的距离,g 是到达当前点已经点的花费,h 是预计剩下的花费。这里 h 取当前点的距离到 s 距离,可通过从 s 跑一遍 Dijkstra 可以预处理得出。

Astar 算法是只有到达终点的时候才能统计答案,这导致可能拓展很多个状态才能得到一个用来更新答案的有效状态。例如一个 n 元环,当我们到达终点之后,可能还要拓展 n 次才能得到下一个状态,于是时间复杂度就被卡到$O(nk)$。

struct Astar : Dijkstra
{
	vector<ll> ans;
	Astar(int n) : Dijkstra(n) {}
	void ask(int s, int t, int k)
	{
		Dijkstra::ask(s);
		ans.assign(k, INF);
		if (d[t] == INF)
			return;
		vector<int> cnt(v.size(), 0);
		priority_queue<pair<ll, int>> q;
		for (q.push(make_pair(-d[t], t)); cnt[s] < k && !q.empty();)
		{
			ll dis = -q.top().first;
			int u = q.top().second;
			if (u == s)
				ans[cnt[s]] = dis;
			if (q.pop(), ++cnt[u] > k)
				continue;
			for (int i = 0, k; i < v[u].i.size(); ++i)
				k = v[u].i[i], q.push(make_pair(d[u] - d[e[k].first] - e[k].len - dis, e[k].first));
		}
	}
};

网络流

ISAP 求最大流

使用示例

struct ISAP : Graph
{
	ll flow;
	vector<ll> f;
	vector<int> h, cur, gap;
	ISAP(int n) : Graph(n) {}
	void add(Edge ed)
	{
		Graph::add(ed);
		swap(ed.first, ed.second), ed.cap = 0;
		Graph::add(ed);
	}
	ll dfs(int s, int u, int t, ll r)
	{
		if (r == 0 || u == t)
			return r;
		ll _f, _r = 0;
		for (int &i = cur[u], k; i < v[u].o.size(); ++i)
			if (k = v[u].o[i], h[u] == h[e[k].second] + 1)
			{
				_f = dfs(s, e[k].second, t, min(r - _r, e[k].cap - f[k]));
				f[k] += _f, f[k ^ 1] -= _f, _r += _f;
				if (_r == r || h[s] >= v.size())
					return _r;
			}
		if (!--gap[h[u]])
			h[s] = v.size();
		return ++gap[++h[u]], cur[u] = 0, _r;
	}
	void ask(int s, int t)
	{
		h.assign(v.size(), 0);
		cur.assign(v.size(), 0);
		gap.assign(v.size() + 2, 0);
		/*
		for (deque<int> q(h[t] = gap[t] = 1, t); !q.empty(); q.pop_front()) //优化,加了能快一点
			for (int i = 0, u = q.front(), k, to; i < v[u].o.size(); ++i)
				if (to = e[v[u].o[i]].second, !h[to])
					++gap[h[to] = h[u] + 1], q.push_back(to);
		*/
		for (f.assign(e.size(), flow = 0); h[s] < v.size();)
			flow += dfs(s, s, t, INF);
	}
};

PrimalDual 求费用流

使用示例,定义一条边的费用为流量*边长,求总费用最小的最大流。性能优秀,只能跑非负权图。

struct PrimalDual : Graph
{
	ll flow, cost;
	vector<ll> f;
	PrimalDual(int n) : Graph(n) {}
	void add(Edge ed)
	{
		Graph::add(ed);
		swap(ed.first, ed.second), ed.cap = 0, ed.len *= -1;
		Graph::add(ed);
	}
	void ask(int s, int t) //询问s到t的最小费用最大流,答案存在flow、cost中
	{
		vector<int> p(v.size(), e.size());
		vector<ll> d(v.size(), INF), h(v.size(), 0);
		for (f.assign(e.size(), flow = cost = 0);;)
		{
			priority_queue<pair<ll, int>> q;
			for (q.push(make_pair(d[s] = 0, s)); !q.empty();)
			{
				ll dis = -q.top().first;
				int u = q.top().second;
				if (q.pop(), d[u] < dis)
					continue;
				for (int i = 0, k, to; i < v[u].o.size(); ++i)
					if (k = v[u].o[i], to = e[k].second,
						e[k].cap > f[k] && d[to] > d[u] + e[k].len + h[u] - h[to])
					{
						d[to] = d[u] + e[k].len + h[u] - h[to], p[to] = k;
						q.push(make_pair(-d[to], to));
					}
			}
			if (d[t] == INF)
				return;
			for (int i = 0; i < d.size(); ++i)
				if (d[i] != INF)
					h[i] += d[i], d[i] = INF;
			ll _f = INF;
			for (int u = t; u != s; u = e[p[u]].first)
				_f = min(_f, e[p[u]].cap - f[p[u]]);
			for (int u = t; u != s; u = e[p[u]].first)
				cost += _f * e[p[u]].len, f[p[u]] += _f, f[p[u] ^ 1] -= _f;
			flow += _f;
		}
	}
};

EK 求费用流

使用示例,定义一条边的费用为流量*边长,求总费用最小的最大流。BellmanFord 算法找增广路,可能被卡但是可以跑负费用流(最大费用流)。

struct EdmondKarp : Graph
{
	ll flow, cost;
	vector<ll> f;
	EdmondKarp(int n) : Graph(n) {}
	void add(Edge ed)
	{
		Graph::add(ed);
		swap(ed.first, ed.second), ed.cap = 0, ed.len *= -1;
		Graph::add(ed);
	}
	void ask(int s, int t)
	{
		vector<int> p(v.size(), e.size());
		for (f.assign(e.size(), flow = cost = 0);;)
		{
			vector<ll> d(v.size(), INF);
			vector<int> flag(v.size(), d[s] = 0);
			for (deque<int> q(flag[s] = 1, s); !q.empty(); q.pop_front())
				for (int u = q.front(), i = flag[u] = 0, k, to; i < v[u].o.size(); ++i)
					if (k = v[u].o[i], to = e[k].second,
						e[k].cap > f[k] && d[to] > d[u] + e[k].len)
					{
						d[to] = d[u] + e[k].len, p[to] = k;
						if (!flag[to])
							q.push_back(to), flag[to] = 1;
					}
			if (d[t] == INF)
				return;
			ll _f = INF;
			for (int u = t; u != s; u = e[p[u]].first)
				_f = min(_f, e[p[u]].cap - f[p[u]]);
			for (int u = t; u != s; u = e[p[u]].first)
				cost += _f * e[p[u]].len, f[p[u]] += _f, f[p[u] ^ 1] -= _f;
			flow += _f;
		}
	}
};

ZKW 求费用流

使用示例,定义一条边的费用为流量*边长,求总费用最小的最大流。

对于最终流量较大,而费用取值范围不大的图,或者是增广路径比较短的图(如二分图),zkw 算法都会比较快,原因是充分发挥优势。比如流多说明可以同一费用反复增广,费用窄说明不用改太多距离标号就会有新增广路,增广路径短可以显著改善最坏情况,因为即使每次就只增加一条边也可以很快凑成最短路。如果恰恰相反,流量不大,费用不小,增广路还较长,就不适合 zkw 算法了。

struct ZKW : Graph
{
	ll flow, cost;
	vector<ll> h, f;
	vector<int> vis;
	ZKW(int n) : Graph(n) {}
	void add(Edge ed)
	{
		Graph::add(ed);
		swap(ed.first, ed.second), ed.cap = 0, ed.len *= -1;
		Graph::add(ed);
	}
	ll dfs(int u, int t, ll r)
	{
		if (r == 0 || u == t)
			return r;
		if (vis[u])
			return 0;
		ll _f = vis[u] = 1, _r = 0;
		for (int i = 0, k; r > _r && i < v[u].o.size(); ++i)
			if (k = v[u].o[i], h[e[k].second] + e[k].len == h[u])
				_f = dfs(e[k].second, t, min(r - _r, e[k].cap - f[k])), f[k] += _f, f[k ^ 1] -= _f, _r += _f;
		return _r;
	}
	void ask(int s, int t)
	{
		h.assign(v.size(), 0);
		vis.assign(v.size(), 0);
		for (f.assign(e.size(), flow = cost = 0);;)
		{
			ll _f = dfs(s, t, INF), d = INF;
			flow += _f, cost += _f * h[s];
			for (int u = 0; u < v.size(); ++u)
				for (int i = 0, k; vis[u] && i < v[u].o.size(); ++i)
					if (k = v[u].o[i], !vis[e[k].second] && e[k].cap > f[k])
						d = min(d, e[k].len + h[e[k].second] - h[e[k].first]);
			if (d == INF)
				return;
			for (int i = 0; i < v.size(); ++i)
				if (vis[i])
					h[i] += d, vis[i] = 0;
		}
	}
};

上下界有源汇网络流

T 向 S 连容量正无穷的边,将有源汇转化为无源汇。每条边容量减去下界,设$in[i]$表示流入 i 的下界之和减去流出 i 的下界之和。新建超级源汇 SS,TT,对于$in[i]>0$的点,SS 向 i 连容量为$in[i]$的边。对于$in[i]<0$的点,i 向 TT 连容量为$−in[i]$的边。

求出以 SS,TT 为源汇的最大流,如果等于$\sum ini$则存在可行流。再求出以 S,T 为源汇的最大流即为最大流。

费用流:建完图后等价于求以 SS,TT 为源汇的的费用流。

上下界费用流:先把下界的费用加入答案。

判断边是否属于某一割集

在残余网络 (还有流量的边) 中求强连通分量,顶点不在同一 SCC 且满流的边。

判断边是否为全部最小割集的边: 在上一条的基础上,还要满足起点与 S 在同一 SCC,且终点与 T 在同一 SCC。

线性规划转费用流

首先添加松弛变量,将不等号都变为等号。分别用下一个式子减去上一个式子,如果每个变量只出现了两次且符号一正一负,那么可以转化为费用流。对于每个式子建立一个点,那么每个变量对应一条边,从一个点流出,向另一个点流入。这样,对于等式右边的常数,如果是正的,对应从源点向该点连一条流量 C,费用 0 的边;如果是负的对应从该点向汇点连一条流量 −C,费用 0 的边。对于每个变量,从它系数为正的式子向系数为负的式子连一条容量为 inf,费用为它在目标函数里系数的边。这样网络流模型就构造完毕了。

欧拉路

使用示例,给定无孤立结点图 G,若存在一条路,经过图中每边一次且仅一次,该条路称为欧拉路。

  • 无向图:当仅当该图所有顶点的度数为偶数(此时为回路),或除两个度数为奇数外(作为路径的起点和终点)、其余全为偶数。
  • 有向图:当仅当该图所有顶点出度=入度(此时为回路),或一个顶点出度=入度+1(作为起点)、另一个顶点入度=出度+1(作为终点)、其他顶点出度=入度。
struct Fleury : Graph
{
	vector<int> vis, cur, p;
	Fleury(int n) : Graph(n) {}
	void dfs(int u)
	{
		for (int &i = cur[u], k; i < v[u].i.size(); ++i) //遍历原图的反向图,这里加了一个「当前弧」优化
			if (k = v[u].i[i], !vis[k] && !vis[k ^ 1])   //无向图需要同时检查反向边未被访问过
			{
				vis[k] = 1;
				dfs(e[k].first);
				p.push_back(k);
			}
	}
	void ask() //查询欧拉回路,路径上边的序号按顺序存在p中
	{
		vis.assign(e.size(), 0), cur.assign(v.size(), 0), p.clear();
		for (int i = 0; i < v.size(); ++i)
			if (v[i].i.size() % 2)
				return dfs(i);
		dfs(0);
	}
};

混合图欧拉回路判定

首先给无向边随便定一个方向,设$\deg x$为 x 连出去的边数 − 连入 x 的边数。若存在$\deg x$为奇数,或者图不连通,则无解。否则建立源点 S,汇点 T。对于一个点 x,若$\deg x>0$,则 S 向 x 连边,容量$\frac{\deg x}{2}$;若$\deg x<0$,则 x 向 T 连边,容量$-\frac{\deg x}{2}$。 对于一条定了向的无向边$x\to y$,x 向 y 连边,容量 1,求出最大流,若与 S 和 T 连的每条边都满流,则有解。

连通性

无向图求割和双连通分量

使用示例 1

使用示例 2

  • 割边:在连通图中,删除了连通图的某条边后,图不再连通。这样的边被称为割边,也叫做桥。
  • 割点:在连通图中,删除了连通图的某个点以及与这个点相连的边后,图不再连通。这样的点被称为割点。构造 dfs 搜索树,在树上有两类节点可以成为割点:
    • 对根节点 u,若其有两棵或两棵以上的子树,则该根结点 u 为割点;
    • 对非根非叶节点 u,若其中的某棵子树的节点均没有指向 u 的祖先节点的回边,说明删除 u 之后,根结点与该棵子树的节点不再连通;则节点 u 为割点。

对于一个无向图的子图,当删除其中任意一条边后,不改变图内点的连通性,这样的子图叫做边的双连通子图。而当子图的边数达到最大时,叫做边的双连通分量。原理是图中所有割边再求一次 SCC,可直接使用求 SCC 的代码。

对于一个无向图的子图,当删除其中任意一个点后,不改变图内点的连通性,这样的子图叫做点的双连通子图。而当子图的边数达到最大时,叫做点的双连通分量。下面给出求点双连通分量的代码。

struct BiconnectedConnectedComponenet : Graph
{
	vector<int> dep, bid, stak, cutPoint, cutEdge; //bid边的端点所属双连通块
	BiconnectedConnectedComponenet(int n) : Graph(n) {}
	void ask()
	{
		dep.assign(v.size(), v.size());
		bid.assign(e.size(), e.size());
		cutPoint.assign(v.size(), 0);
		cutEdge.assign(e.size(), 0);
		for (int i = 0; i < v.size(); ++i)
			if (dep[i] == v.size())
				dfs(i, v.size());
	}
	int dfs(int u, int fa)
	{
		int low = dep[u] = fa != v.size() ? dep[fa] + 1 : 0;
		for (int i = 0, k, to, lowto; i < v[u].o.size(); ++i)
			if (k = v[u].o[i], to = e[k].second, to != fa)
			{
				if (dep[to] == v.size())
				{
					stak.push_back(k);
					low = min(low, lowto = dfs(to, u));
					if (lowto > dep[u])
						cutEdge[k] = cutEdge[k ^ 1] = 1;
					if (lowto >= dep[u])
						for (cutPoint[u] = fa != v.size() || i;;)
						{
							int x = stak.back();
							stak.pop_back();
							bid[x] = bid[x ^ 1] = k;
							if (x == k)
								break;
						}
				}
				else if (dep[to] < dep[u])
				{
					stak.push_back(k);
					low = min(low, dep[to]);
				}
			}
		return low;
	}
};

双连通图的构造

先求出所有的桥,然后删除这些桥边,剩下的每个连通块都是一个双连通子图。把每个双连通子图收缩为一个顶点,再把桥边加回来,最后的这个图一定是一棵树,边连通度为 1。统计出树中度为 1 的节点的个数,即为叶节点的个数,记为leaf。至少在树上添加(leaf+1)/2条边,就能使树达到边双连通:先把两个最近公共祖先最远的两个叶节点之间连接一条边,这样可以把这两个点到祖先的路径上所有点收缩到一起,因为一个形成的环一定是双连通的;然后再找两个最近公共祖先最远的两个叶节点,这样一对一对找完,恰好是(leaf+1)/2次,把所有点收缩到了一起。

无向图求边双连通分量&有向图求强连通分量

无向图求边双连通分量

有向图求强连通分量

struct StronglyConnectedComponenet : Graph
{
	vector<int> dep, sid, stak; //sid=点所属连通块内一点
	StronglyConnectedComponenet(int n) : Graph(n) {}
	void ask()
	{
		dep.assign(v.size(), v.size());
		sid.assign(v.size(), v.size());
		for (int i = 0; i < v.size(); ++i)
			if (dep[i] == v.size())
				dfs(i, v.size());
	}
	int dfs(int u, int fa)
	{
		int low = dep[u] = fa != v.size() ? dep[fa] + 1 : 0;
		stak.push_back(u);
		for (int i = 0, k, to; i < v[u].o.size(); ++i)
			if (k = v[u].o[i], to = e[k].second, to != fa /*, 1*/) //求强连通分量把注释去掉,即允许走回边
			{
				if (dep[to] == v.size())
					low = min(low, dfs(to, u));
				else if (sid[to] == v.size())
					low = min(low, dep[to]);
			}
		if (low == dep[u])
			for (;;)
			{
				int x = stak.back();
				stak.pop_back();
				sid[x] = u;
				if (x == u)
					break;
			}
		return low;
	}
};

2-SAT

使用示例,n 个布尔变量$x_0\dots x_{n-1}$,逻辑表达式$Y=(A_0+B_0)(A_1+B_1)\dots(A_{m-1}+B_{m-1})$,其中$A_i,B_i\in\lbrace x_j,\overline{x_j}\rbrace $,判断是否存在$x_0\dots x_{n-1}$的取值使得 Y 值为 1。因为$A+B=(\overline A\to B)(\overline B\to A)$,所以对于一个要求$A+B$,我们连$\overline A\to B,\overline B\to A$两条边。如果有一条边$A\to B$,意味着如果 A 成立那么 B 必然成立。若$\exists i,x_i,\overline{x_i}\in$同一 SCC,则不存在。

模板要求顶点数v.size()是偶数,ii ^ 1互为相反变量且奇数下标代表布尔变量为真的情况。最后输出结果在ok中。不满足对任意 iok[i] ^ ok[i ^ 1] == 1时,ask函数返回false

struct TwoSat : StronglyConnectedComponenet
{
	vector<int> ok;
	TwoSat(int n) : StronglyConnectedComponenet(n) {}
	void addOR(int a, int b) { add({a ^ 1, b}), add({b ^ 1, a}); }
	void addXOR(int a, int b) { addOR(a, b), addOR(a ^ 1, b ^ 1); }
	void addXNOR(int a, int b) { addOR(a, b ^ 1), addOR(a ^ 1, b); }
	int ask()
	{
		StronglyConnectedComponenet::ask();
		for (int i = 0; i < v.size(); i += 2)
			if (sid[i] == sid[i ^ 1])
				return 0;
		vector<vector<int>> g(v.size());
		vector<int> ind(v.size(), 0), cf(v.size(), 0), stak;
		for (int i = 0; i < v.size(); ++i)
		{
			cf[sid[i]] = sid[i ^ 1];
			for (int j = 0, k; j < v[i].o.size(); ++j)
				if (sid[e[k = v[i].o[j]].second] != sid[i])
					g[sid[e[k].second]].push_back(sid[i]), ++ind[sid[i]];
		}
		for (int i = 0; i < v.size(); ++i)
			if (sid[i] == i && !ind[i])
				stak.push_back(i);
		for (ok.assign(v.size(), -1); !stak.empty();)
		{
			int u = stak.back();
			stak.pop_back();
			if (ok[u] < 0)
				ok[u] = 1, ok[cf[u]] = 0;
			for (int i = 0; i < g[u].size(); ++i)
				if (!--ind[g[u][i]])
					stak.push_back(g[u][i]);
		}
		for (int i = 0; i < v.size(); ++i)
			if (i != sid[i])
				ok[i] = ok[sid[i]];
		return 1;
	}
};
2-SAT 暴力搜索求字典序最小的解

使用示例,时间复杂度最坏为为$O(VE)$。

struct TwoSat : Graph
{
	vector<int> ok;
	TwoSat(int n) : Graph(n) {}
	void addOR(int a, int b) { add({a ^ 1, b}), add({b ^ 1, a}); }
	void addXOR(int a, int b) { addOR(a, b), addOR(a ^ 1, b ^ 1); }
	void addXNOR(int a, int b) { addOR(a, b ^ 1), addOR(a ^ 1, b); }
	int dfs(int u, vector<int> &stak)
	{
		if (ok[u ^ 1])
			return 0;
		if (ok[u])
			return 1;
		ok[u] = 1;
		stak.push_back(u);
		for (int i = 0; i < v[u].o.size(); ++i)
			if (!dfs(e[v[u].o[i]].second, stak))
				return 0;
		return 1;
	}
	int ask()
	{
		ok.assign(v.size(), 0);
		for (int i = 0; i < v.size(); i += 2)
			if (!ok[i] && !ok[i ^ 1])
			{
				vector<int> stak;
				if (!dfs(i, stak))
				{
					for (int j = 0; j < stak.size(); ++j)
						ok[stak[j]] = 0;
					if (!dfs(i ^ 1, stak))
						return 0;
				}
			}
		return 1;
	}
};

二分图

二分图的一个等价定义是:不含有含奇数条边的环的图。

完美匹配图中所有的顶点都是匹配点。

二分图的最小点覆盖(最小割)是指最少的顶点数使得二分图 G 中的每条边都至少与其中一个点相关联。二分图中,最小割=最大匹配。

二分图的最小边覆盖(最大独立集)是指用尽量少的不相交简单路径覆盖二分图中的所有顶点。二分图中,最小点覆盖+最小边覆盖=总点数。

Hall 定理:二分图中的两部分顶点组成的集合分别为 X,Y ,则有一组无公共点的边,一端恰好为组成 X 的点的充分必要条件是:X 中的任意 k 个点至少与 Y 中的 k 个点相邻。对于区间图只需要考虑极端情况,线段树维护。

关键点:一定在最大匹配中的点。 求出任意一个最大匹配,那么只需要考虑哪些匹配点不一定在。 假设是考虑左侧的点,右侧类似: 将匹配边从右往左,非匹配边从左到右,从左侧每个未匹配点开始 DFS 到的匹配点都不是关键点。

关键边:求出任意一个最大匹配,将匹配边从右到左,剩余边从左到右,求出 SCC。 对于一条边:若它位于当前匹配中,那么若两端点位于同一 SCC,则是可能在,否则必定在;若它不位于当前匹配中,那么若两端点位于同一 SCC,则是可能在,否则必定不在。

Hungary 求最大匹配

使用示例

左边 nl 个点$0\dots nl-1$,右边 nr 个点$0\dots nr-1$,取n=max(nl,nr),左 i 右 j 间相连则g[i].push_back(j);

生成fl[i]表示左边第 i 个匹配右边第fl[i]个,fr同理。时间复杂度$O(n^3)$。未匹配的值为fl[i] == N。匹配是一个边集,其中任意两条边都没有公共顶点;扫一遍flfr判断有多少不等于N即为最大匹配数。

struct Hungary
{
	vector<int> g[N];
	int n, fl[N], fr[N], vr[N];
	int dfs(int x, int rt)
	{
		for (auto y : g[x])
			if (vr[y] != rt)
				if (vr[y] = rt, fr[y] == N || dfs(fr[y], rt))
					return fl[fr[y] = x] = y, 1;
		return 0;
	}
	void ask()
	{
		fill(fl, fl + n, N), fill(fr, fr + n, N), fill(vr, vr + n, N);
		for (int i = 0; i < n; ++i)
			if (fl[i] == N)
				dfs(i, i); //找完全匹配时如果返回值为0可直接return
	}
};

HopcroftKarp 求最大匹配

使用示例,时间复杂度$O(\sqrt{\vert V\vert }\vert E\vert )$,稀疏图上效果明显。

struct HopcroftKarp
{
	vector<int> g[N];
	int n, fl[N], fr[N], hl[N], hr[N], q[N];
	bool dfs(int x)
	{
		int c = hl[x] + 1, y = hl[x] = N;
		for (auto y : g[x])
			if (hr[y] == c)
				if (hr[y] = N, fr[y] == N || dfs(fr[y]))
					return fl[fr[y] = x] = y, 1;
		return 0;
	}
	void ask()
	{
		fill(fl, fl + n, N), fill(fr, fr + n, N);
		for (int x = 0; x < n; ++x)
			for (auto y : g[x])
				if (fr[y] == N)
				{
					fl[fr[y] = x] = y;
					break;
				}
		for (int ql, qr, ok;;)
		{
			fill(hl, hl + n, N), fill(hr, hr + n, N);
			for (int x = ql = qr = ok = 0; x < n; ++x)
				if (fl[x] == N)
					hl[q[qr++] = x] = 0;
			while (ql < qr)
				for (auto y : g[q[ql++]])
					if (hr[y] == N)
					{
						hr[y] = hl[q[ql - 1]] + 1;
						if (fr[y] == N)
							ok = 1;
						else if (hl[fr[y]] == N)
							hl[q[qr++] = fr[y]] = hr[y] + 1;
					}
			if (!ok)
				return;
			for (int x = 0; x < n; ++x)
				if (fl[x] == N)
					dfs(x);
		}
	}
};

KuhnMunkres 求最优完备匹配

使用示例

左 x 右 y 代价a[x][y]。最大费用流时,a 初始化 0;最大费用最大流时,a 初始化-N*INF

struct KuhnMunkres
{
	ll a[N][N], hl[N], hr[N], slk[N];
	int n, fl[N], fr[N], vl[N], vr[N], pre[N], q[N], ql, qr;
	int check(int i)
	{
		if (vl[i] = 1, fl[i] != N)
			return vr[q[qr++] = fl[i]] = 1;
		while (i != N)
			swap(i, fr[fl[i] = pre[i]]);
		return 0;
	}
	void bfs(int s)
	{
		fill(slk, slk + n, INF), fill(vl, vl + n, 0), fill(vr, vr + n, 0);
		for (vr[q[ql = 0] = s] = qr = 1;;)
		{
			for (ll d; ql < qr;)
				for (int i = 0, j = q[ql++]; i < n; ++i)
					if (!vl[i] && slk[i] >= (d = hl[i] + hr[j] - a[i][j]))
						if (pre[i] = j, d)
							slk[i] = d;
						else if (!check(i))
							return;
			ll d = INF;
			for (int i = 0; i < n; ++i)
				if (!vl[i] && d > slk[i])
					d = slk[i];
			for (int i = 0; i < n; ++i)
			{
				if (vl[i])
					hl[i] += d;
				else
					slk[i] -= d;
				if (vr[i])
					hr[i] -= d;
			}
			for (int i = 0; i < n; ++i)
				if (!vl[i] && !slk[i] && !check(i))
					return;
		}
	}
	void ask()
	{
		fill(fl, fl + n, N), fill(fr, fr + n, N), fill(hr, hr + n, 0);
		for (int i = 0; i < n; ++i)
			hl[i] = *max_element(a[i], a[i] + n);
		for (int j = 0; j < n; ++j)
			bfs(j);
	}
};

带花树

一般图最大匹配

使用示例

struct Blossom : Graph
{
	vector<int> f;
	Blossom(int n) : Graph(n) {}
	void ask()
	{
		f.assign(v.size(), v.size());
		vector<int> vis(v.size());
		for (int s = 0, cnt = vis.back(); s < v.size(); ++s)
			if (f[s] == v.size())
			{
				UnionfindSet ufs(v.size());
				vector<int> pre(v.size(), v.size()), flag(v.size(), 2);
				for (deque<int> q(flag[s] = 1, s); !q.empty() && f[s] == v.size(); q.pop_front())
					for (int i = 0, x = q.front(), y, a, b; i < v[x].o.size(); ++i)
						if (y = e[v[x].o[i]].second, y != f[x] && flag[y] && ufs.ask(x) != ufs.ask(y))
						{
							if (flag[y] == 1)
							{
								for (a = x, b = y, ++cnt;; swap(a, b))
									if (a != v.size())
									{
										if (vis[a = ufs.ask(a)] == cnt)
											break;
										vis[a] = cnt, a = f[a] != v.size() ? pre[f[a]] : v.size();
									}
								if (ufs.ask(x) != a)
									pre[x] = y;
								if (ufs.ask(y) != a)
									pre[y] = x;
								for (int p[2] = {x, y}, j = 0; j < 2; ++j)
									for (int x = p[j], y, z; x != a; ufs.merge(y, x), ufs.merge(x = z, y))
									{
										if (ufs.ask(z = pre[y = f[x]]) != a)
											pre[z] = y;
										if (!flag[y])
											flag[y] = 1, q.push_back(y);
										if (!flag[z])
											flag[z] = 1, q.push_back(z);
									}
							}
							else if (f[y] != v.size())
								pre[y] = x, flag[y] = 0, flag[f[y]] = 1, q.push_back(f[y]);
							else
							{
								for (pre[y] = x; y != v.size();)
									swap(y, f[f[y] = pre[y]]);
								break;
							}
						}
			}
	}
};

树形图

最小生成树

无向图

同时给出 Prim 算法(生成新树)、Kruskal 算法(消耗小)。

struct Prim : Graph
{
	struct DisGreater
	{
		bool operator()(const Edge &e1, const Edge &e2)
		{
			return e1.len > e2.len;
		}
	};
	ll ans;
	vector<int> vis;
	priority_queue<Edge, vector<Edge>, DisGreater> q;
	Prim(const Graph &g, int root) : Tree(n), ans(0), vis(g.v.size(), 0) //生成新树,每条边都要有等长反向边
	{
		for (insert(root, g); !q.empty();)
		{
			Edge ed = q.top();
			if (q.pop(), !vis[ed.second])
				insert(ed.second, g), ans += ed.len, add(ed);
		}
	}
	void insert(int u, const Graph &g) //把点和对应的相连的边加入集合
	{
		vis[u] = 1;
		for (int i = 0, k; i < g.v[u].o.size(); ++i)
			if (k = g.v[u].o[i], !vis[g.e[k].second])
				q.push(g.e[k]);
	}
};
ll kruskal(vector<Edge> &e, int n) //会清空边集e,每条边被认作无向边
{
	ll ret = 0;
	UnionFindSet ufs(n);
	for (sort(e.begin(), e.end(), DisGreater()); !e.empty(); e.pop_back())
		if (ufs.fa(e.back().from) != ufs.fa(e.back().to))
		{
			ufs.merge(e.back().from, e.back().to);
			ret += e.back().len;
		}
	return /*ufs.siz>1?INF:*/ ret; //视情况选择去注释
}

有向图

使用示例

指定以 root 为根,如果没有限定根那么新建一个虚拟点作为根,向所有边连边长最大边长+1的边,在最后生成的图中去掉此边。时间复杂度$O(VE)$。

ll zhuLiu(vector<Edge> &e, int root, int n) //不存在返回INF
{
	for (ll ret = 0;;)
	{
		vector<ll> in(n, INF);
		vector<int> pre(n, n);
		for (int i = 0, to; i < e.size(); ++i)
		{
			if (e[i].first == (to = e[i].second))
				swap(e[i--], e.back()), e.pop_back();
			else if (in[to] > e[i].len)
				in[to] = e[i].len, pre[to] = e[i].first;
		}
		for (int i = in [root] = 0; i < n; ++i)
			if (in[i] == INF)
				return INF;
		vector<int> id(n, n), vis(n, n);
		int tn = 0;
		for (int i = 0, v; i < n; ++i)
		{
			for (ret += in [v = i]; vis[v] != i && id[v] == n && v != root; v = pre[v])
				vis[v] = i;
			if (v != root && id[v] == n)
			{
				for (int u = pre[v]; u != v; u = pre[u])
					id[u] = tn;
				id[v] = tn++;
			}
		}
		if (!tn)
			return ret;
		for (int i = 0; i < n; ++i)
			if (id[i] == n)
				id[i] = tn++;
		for (int i = 0, v; i < e.size(); ++i)
			if ((e[i].first = id[e[i].first]) != (e[i].second = id[v = e[i].second]))
				e[i].len -= in[v];
		n = tn, root = id[root];
	}
}

树链剖分与 LCA

使用示例

struct Diagram : Graph
{
	Fenwick data; //暂用树状数组作为默认数据结构
	Diagram(const Graph &g, int root) : Graph(g.v.size()), data(g.v.size())
	{
		build(root, g);
		int cnt = v[root].dfn = v[root].dep = 1;
		dfs(v[root].top = root, cnt);
	}
	void build(int u, const Graph &g) //无向图dfs建树,且重边在最前,u为根节点
	{
		v[u].siz = 1;
		for (int i = 0, k, to; i < g.v[u].o.size(); ++i)
			if (k = g.v[u].o[i], to = g.e[k].second, !v[to].siz) //没访问过的点siz默认0
			{
				build(to, g);
				v[u].siz += v[to].siz;
				Graph::add(g.e[k]);
				if (v[ch(u)].siz < v[to].siz) //重边移到最前
					swap(v[u].o.front(), v[u].o.back());
			}
	}
	void dfs(int u, int &cnt)
	{
		for (int i = 0, to; i < v[u].o.size(); ++i)
		{
			v[to = ch(u, i)].dfn = ++cnt;
			v[to].top = i ? to : v[u].top;
			v[to].dep = v[u].dep + 1;
			dfs(to, cnt);
		}
	}
	int lca(int x, int y)
	{
		for (; v[x].top != v[y].top; x = fa(v[x].top))
			if (v[v[x].top].dep < v[v[y].top].dep)
				swap(x, y);
		if (v[x].dep < v[y].dep)
			swap(x, y);
		return y;
	}
	ll ask(int x, int y)
	{
		ll ans = 0;
		for (; v[x].top != v[y].top; x = fa(v[x].top))
		{
			if (v[v[x].top].dep < v[v[y].top].dep)
				swap(x, y);
			ans += data.ask(v[v[x].top].dfn, v[x].dfn);
		}
		if (v[x].dep < v[y].dep)
			swap(x, y);
		return ans += data.ask(v[y].dfn, v[x].dfn);
	}
	void add(int x, int y, ll pv)
	{
		for (; v[x].top != v[y].top; x = fa(v[x].top))
		{
			if (v[v[x].top].dep < v[v[y].top].dep)
				swap(x, y);
			data.add(v[v[x].top].dfn, v[x].dfn, pv);
		}
		if (v[x].dep < v[y].dep)
			swap(x, y);
		data.add(v[y].dfn, v[x].dfn, pv);
	}
};

点剖(点分治)

使用示例,零号点为虚节点。

struct TreeDiv : Graph
{
	int root;
	vector<int> mx, siz, vis;
	TreeDiv(int n) : Graph(n), mx(n, n), siz(n), vis(n) {}
	void dfsRoot(int u, int fa)
	{
		for (int i = mx[u] = siz[u] = 0, k, to; i < v[u].o.size(); ++i)
			if (k = v[u].o[i], to = e[k].second, to != fa && !vis[to])
				if (dfsRoot(to, u), siz[u] += siz[to], mx[u] < siz[to])
					mx[u] = siz[to];
		if (mx[u] < mx[0] - ++siz[u])
			mx[u] = mx[0] - siz[u];
		if (mx[root] >= mx[u])
			root = u;
	}
	void dfsDis(int u, int fa, ll d)
	{
		//用d更新答案
		for (int i = 0, k, to; i < v[u].o.size(); ++i)
			if (k = v[u].o[i], to = e[k].second, to != fa && !vis[to])
				dfsDis(to, u, d + e[k].len);
	}
	int cal(int u, ll d) //返回符合要求的点对数
	{
		return dfsDis(u, 0, d), /*得到答案*/;
	}
	void dfs(int u = 1)
	{
		dfsRoot(u, root = 0), ans += cal(u = root, 0), vis[u] = 1;
		for (int i = 0, k, to; i < v[u].o.size(); ++i)
			if (k = v[u].o[i], to = e[k].second, !vis[to])
				ans -= cal(to, e[k].len), mx[0] = siz[to], dfs(to);
	}
};

数论

辗转相除法

ll gcd(ll a, ll b) { return b ? gcd(b, a % b) : a; } //a、b的最大公约数
ll lcm(ll a, ll b) { return a / gcd(a, b) * b; }	 //a、b的最小公倍数
ll gcd(ll a, ll b, ll &x, ll &y)					 //扩展欧几里得,引用返回a*x+b*y=gcd(a,b)绝对值之和最小的解
{
	if (!a)
		return x = 0, y = 1, b;
	ll d = gcd(b % a, a, y, x);
	return x -= b / a * y, d;
}
ll gcd(ll a, ll b) //无取模求gcd
{
	for (ll t = 1, c, d;;)
	{
		if (a == b)
			return t * a;
		if (a < b)
			swap(a, b);
		if (a & 1)
			c = 0;
		else
			a >>= 1, c = 1;
		if (b & 1)
			d = 0;
		else
			b >>= 1, d = 1;
		if (c && d)
			t <<= 1;
		else if (!c && !d)
			a -= b;
	}
}

裴蜀定理

对任何$a,b\in Z$和它们的最大公约数$d$,关于未知数$x,y$的线性不定方程(称为裴蜀等式):$ax+by=c$当仅当$d\vert c$,可知有无穷多解。特别地,$ax+by=d$一定有解。

推论

$a,b$互质的充要条件是$ax+by=1$有整数解。

同余系运算

三种乘法的比较,结论是如果可以用__int128就用,否则就用第一种:

求乘法逆元的另外一种方法是用欧拉定理$x^{\phi(m)}\equiv1\pmod m$,x 的逆是$x^{\phi(m)-1}$。特别地,m 为素数时$\phi(m)=m-1$,此时 x 的逆就是pow(x,m-2,m)

log 函数:m 为素数时求解模方程$a^x\equiv b\pmod m$。设 P 为质数,G 为 P 的原根,则$x^y\equiv b\pmod P$等价于$y\ ind\ x\equiv b\pmod{P-1}$,其中$G\ ind\ x\equiv x\pmod P$。

struct Mod
{
	const ll M, SM;
	Mod(ll M) : M(M), SM(sqrt(M) + 0.5) {}
	ll qadd(ll &a, ll b) const { return a += b, a >= M ? a -= M : a; } //假如a+b<2*M,就不必取模了,取模运算耗时很高
	ll add(ll a, ll b) const { return qadd(a = (a + b) % M, M); }	  //考虑a和b不在同余系内甚至为负数的情况
	ll mul(ll a, ll b) const { return add(a * b, M); } //ll类型是int的时候小心 a*b 爆精度
	ll inv(ll a) const { return pow(a, M - 2); } //要求M为素数,否则return pow(a, phi(M) - 1);
	ll pow(ll a, ll b) const
	{
		ll r = 1;
		/*
		if (b < 0)
			b = -b, a = inv(a);
		*/
		for (a = add(a, M); b; b >>= 1, a = mul(a, a))
			if (b & 1)
				r = mul(r, a);
		return r;
	}
	/*
	ll mul(ll a, ll b) const { return add(a * b, -M * (ll)((long double)a / M * b)); } //long double 有效位数16~18位,模数过大时慎用!
	ll mul(ll a, ll b) const //Head算法,无循环快速计算同余乘法,根据a*b是否爆ll替换a*b%M,需要a<M且b<M,可以调用时手动取模
	{
		ll c = a / SM, d = b / SM;
		a %= SM, b %= SM;
		ll e = add(add(a * d, b * c), c * d / SM * (SM * SM - M));
		return add(add(a * b, e % SM * SM), add(c * d % SM, e / SM) * (SM * SM - M));
	}
	ll mul(ll a, ll b) const //龟速乘
	{
		ll r = 0;
		for (a %= M; b; b >>= 1, qadd(a, a))
			if (b & 1)
				qadd(r, a);
		return r;
	}
	ll inv(ll a) const //模m下a的乘法逆元,不存在返回-1(m为素数时a不为0必有逆元)
	{
		ll x, y, d = gcd(a, M, x, y);
		return d == 1 ? add(x, M) : -1;
	}
	vector<ll> sol(ll a, ll b) const //解同余方程,返回ax=b(mod M)循环节内所有解
	{
		vector<ll> ans;
		ll x, y, d = gcd(a, M, x, y);
		if (b % d)
			return ans;
		ans.push_back(mul((b / d) % (M / d), x));
		for (ll i = 1; i < d; ++i)
			ans.push_back(add(ans.back(), M / d));
		return ans;
	}
	ll log(ll a, ll b) const
	{
		unordered_map<ll, ll> x;
		for (ll i = 0, e = 1; i <= SM; ++i, e = mul(e, a))
			if (!x.count(e))
				x[e] = i;
		for (ll i = 0, v = inv(pow(a, SM)); i <= SM; ++i, b = mul(b, v))
			if (x.count(b))
				return i * SM + x[b];
		return -1;
	}
	*/
};

中国剩余定理解同余方程组

使用示例

ll crt(const vector<pair<ll, ll>> &v) //同余方程组,x%v[i].first==v[i].second,不存在返回-1
{
	ll m = v[0].first, r = v[0].second, c, d, x, y, z;
	for (int i = 1; i < v.size(); ++i)
	{
		if (c = v[i].second - r, d = gcd(m, v[i].first, x, y), c % d)
			return -1;
		gcd(m / d, z = v[i].first / d, x, y), r += c / d * x % z * m, r %= m *= z;
	}
	return r < 0 ? r + m : r;
}
struct Excrt
{
	vector<ll> b2, a2;
	int n;
	void init()
	{
		scanf("%d", &n);
		ll q, qq;
		for (int i = 0; i < n; i++)
			scanf("%lld%lld", &q, &qq), b2.push_back(q), a2.push_back(qq);
	}
	ll get()
	{
		ll M, ans, x = 0, y = 0;
		M = 1, ans = 0;
		for (int i = 0; i < n; i++)
		{
			ll a = M, b = b2[i], c = (a2[i] - ans % b + b) % b; //ax=c(mod b)
			ll gcd = exgcd(a, b, x, y), bg = b / gcd;
			if (c % gcd != 0)
			{
				return -1;
			} //返回-1表示无解
			x = mul(x, c / gcd, bg);
			y = M / gcd * b;
			ans = (x * M + ans) % y; //???k???????
			if (ans < 0)
				ans += y;
			M = y; //M??k?m?lcm
		}
		return ans;
	}
};

二次剩余

struct QuadraticRes
{
	struct numb
	{
		ll x, y;
	};
	numb mul(const numb &a, const numb &b, ll p)
	{
		return (numb){(a.x * b.x % p + a.y * b.y % p * w % p) % p, (a.x * b.y + a.y * b.x) % p};
	}
	ll pow2(numb a, ll b, ll p)
	{
		numb ret = {1, 0}, tem = a;
		while (b)
		{
			if (b & 1)
				ret = mul(ret, tem, p);
			tem = mul(tem, tem, p), b >>= 1;
		}
		return ret.x % p;
	}
	ll get(ll n, ll p)
	{
		if (Mod(p).pow(n, p / 2) == p - 1)
		{
			return -1;
		} //此时无解
		while (1)
		{
			ll a = rand();
			w = (a * a - n + p) % p;
			if (Mod(p).pow(w, p / 2) == p - 1)
			{
				return pow2((numb){a, 1}, p / 2 + 1, p); //还有一个解是p减去这个值
			}
		}
	}
};

欧拉筛

欧拉函数$\phi(n)$是小于 n 的正整数中与 n 互素的数的数目。特别地,规定$\phi(1)=1$,易知 n>2 时都为偶数。

欧拉函数是积性函数,即对任意素数$p,q$满足下列关系:$\phi(pq)=\phi(p)\phi(q)=(p-1)(q-1)$对任何两个互质的正整数$x, m(m\geq2)$有欧拉定理:$x^{\phi(m)}\equiv1\pmod m$当 m 为素数 p 时,此式变为费马小定理:$x^{p-1}\equiv1\pmod p$利用欧拉函数和它本身不同质因数的关系,用筛法$O(N)$预处理某个范围内所有数的欧拉函数值,并求出素数表。同时,利用计算欧拉函数过程中求出的最小素因子 m,可以实现$O(log N)$的素因数分解。

同时求莫比乌斯函数$\mu(n)$,存在mu中。

struct EulerSieve
{
	vector<int> p, m, phi, mu; //素数序列,最小素因子,欧拉函数,莫比乌斯函数
	EulerSieve(int N) : m(N, 0), phi(N, 0), mu(N, 0)
	{
		phi[1] = mu[1] = 1;					 //m[1]=0,m[i]==i可判断i是素数
		for (long long i = 2, k; i < N; ++i) //防i*p[j]爆int
		{
			if (!m[i])
				p.push_back(m[i] = i), phi[i] = i - 1, mu[i] = -1; //i是素数
			for (int j = 0; j < p.size() && (k = i * p[j]) < N; ++j)
			{
				phi[k] = phi[i] * p[j];
				if ((m[k] = p[j]) == m[i])
				{
					mu[k] = 0;
					break;
				}
				phi[k] -= phi[i];
				mu[k] = -mu[i];
			}
		}
	}
};

直接求欧拉函数

ll phi(ll n)
{
	ll phi = n;
	for (ll i = 2; i * i <= n; ++i)
		if (!(n % i))
			for (phi = phi / i * (i - 1); !(n % i);)
				n /= i;
	if (n > 1)
		phi = phi / n * (n - 1);
	return phi;
}

常见数论函数变换

$\sum_{d\vert n}\mu(d)=[n=1]$

$\phi(n)=\sum_{i=1}^n[\gcd(i,n)=1]=\sum_{i=1}^n\sum_{k\mid i,k\mid n}\mu(k)=\sum_{k\mid n}\frac nk\mu(k)$

前缀和

欧拉函数前缀和$S_\phi(n)=\frac{(n+1)n}2-\sum_{d=1}^nS_\phi(\frac{n}{d})$

莫比乌斯函数前缀和$S_\mu(n)=1-\sum_{d=1}^nS_\mu(\frac{n}{d})$

莫比乌斯反演

若$f(n)=\sum_{d\vert n}g(d)$,则$g(n)=\sum_{d\vert n}\mu(d)f(\frac{n}{d})$

若$f(n)=\sum_{i=1}^nt(i)g(\frac{n}{i})$,则$g(n)=\sum_{i=1}^n\mu(i)t(i)f(\frac{n}{i})$(此时代$t(i)=[gcd(n,i)>1]$可得上式)

举例(其中$T=kd$):

[\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)
=\sum_d d\sum_{i=1}^n\sum_{j=1}^m[\gcd (i,j)=d]
=\sum_{d}d\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}[\gcd (i,j)=1]
=\sum_{d}d\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}\sum_{k\mid i,k\mid j}\mu(k)
=\sum_d d\sum_k\mu(k)\sum_{k\mid i}^{\lfloor\frac nd\rfloor}\sum_{k\mid j}^{\lfloor\frac md\rfloor}
=\sum_{d}d\sum_k\mu(k)\lfloor\frac n{kd}\rfloor\lfloor\frac m{kd}\rfloor
=\sum_{T}\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor\sum_{k\mid T}\frac Tk\mu(k)
=\sum_{T}\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor\varphi(T)]

$\varphi(T)$可以使用线性筛预处理处理,我们就可以枚举$T$求上式了,时间复杂度$O(n)$。多组数据$n,m$询问上式,时间复杂度就变成了$O(Tn)$。事实上,$\lfloor\frac{n}{T}\rfloor$是不会轻易变化的,是过了连续的一段后才发生变化的,那么我们就可以计算出这一段的结束位置,对$\varphi$函数作前缀和,就可以直接分块了,这样的时间复杂度是$O(T\sqrt{n})$的。

PollardRho 大数素因子分解

时间复杂度$O(N^{1/4})$,数据多的时候可考虑欧拉筛优化。

注意这里模乘法很容易爆 long long,看情况选用快速乘法。

struct PollardRho
{
	bool isPrime(ll n, int S = 12) //MillerRabin素数测试,S为测试次数,用前S个素数测一遍,S=12可保证unsigned long long范围内无错;n<2请特判
	{
		static ll d, u, t, p[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
		for (d = n - 1; !(d & 1);)
			d >>= 1; //未对0,1做特判!
		Mod mo(n);
		for (int i = 0; i < S; ++i)
		{
			if (!(n % p[i]))
				return n == p[i];
			for (t = mo.pow(p[i], u = d); t != n - 1 && t != 1 && u != n - 1;)
				t = mo.mul(t, t), u <<= 1;
			if (t != n - 1 && !(u & 1))
				return 0;
		}
		return 1;
	}
	void fac(ll n, vector<ll> &factor)
	{
		/*
		if (n < e.m.size()) //欧拉筛预处理优化,可以防止很多小素数因子的情况
		{
			for (; n > 1; n /= e.m[n])
				factor.push_back(e.m[n]);
			return;
		}
		*/
		if (isPrime(n))
			return factor.push_back(n);
		Mod mo(n);
		for (ll c = 1;; ++c)
			for (ll i = 0, k = 1, x = rand() % (n - 1) + 1, y, p;;)
			{
				if (++i == k)
					y = x, k <<= 1;
				if (x = mo.add(mo.mul(x, x), c), p = __gcd(abs(x - y), n), p == n)
					break;
				if (p > 1)
					return fac(p, factor), fac(n / p, factor);
			}
	}
};

快速变换

蝴蝶变换(雷德变换)

保存 FTT 和 FNTT 时交换的对应位置(即保存的是置换)。

struct Rader : vector<int>
{
	Rader(int n) : vector<int>(1 << int(ceil(log2(n))))
	{
		for (int i = at(0) = 0; i < size(); ++i)
			if (at(i) = at(i >> 1) >> 1, i & 1)
				at(i) += size() >> 1;
	}
};

快速傅里叶变换

使用示例,高精度乘法。

任意模数 FFT

struct FFT : Rader
{
	vector<complex<lf>> w;
	FFT(int n) : Rader(n), w(size())
	{
		for (int i = 0; i < size(); ++i)
			w[i] = polar(1.0, 2 * M_PI * i / size());
	}
	void fft(vector<complex<lf>> &x) const
	{
		for (int i = 0; i < x.size(); ++i)
			if (i < at(i))
				std::swap(x[i], x[at(i)]);
		for (int i = 1; i < size(); i <<= 1)
			for (int j = 0; j < i; ++j)
				for (int k = j; k < size(); k += i << 1)
				{
					complex<lf> t = w[size() / (i << 1) * j] * x[k + i];
					x[k + i] = x[k] - t, x[k] += t;
				}
	}
	vector<ll> ask(const vector<ll> &a, const vector<ll> &b) const
	{
		vector<complex<lf>> xa(a.begin(), a.end()), xb(b.begin(), b.end());
		xa.resize(size()), xb.resize(size()), fft(xa), fft(xb);
		for (int i = 0; i < size(); ++i)
			xa[i] *= xb[i];
		fft(xa);
		vector<ll> ans(size());
		for (int i = 0; i < size(); ++i)
			ans[i] = xa[(size() - i) & (size() - 1)].real() / size() + 0.5;
		return ans;
	}
	vector<ll> askMod(const vector<ll> &a, const vector<ll> &b, ll M) const //任意模数
	{
		vector<complex<lf>> e(size()), c(size());
		for (int i = 0; i < a.size(); ++i)
			e[i].real(a[i] & 0x7fff), e[i].imag(a[i] >> 15);
		for (int i = 0; i < b.size(); ++i)
			c[i].real(b[i] & 0x7fff), c[i].imag(b[i] >> 15);
		fft(e), fft(c);
		vector<complex<lf>> d(c);
		for (int i = 0; i < size(); ++i)
		{
			int fr = (size() - i) & (size() - 1);
			c[i] *= 0.5 * complex<lf>(e[i].X + e[fr].X, e[i].Y - e[fr].Y);
			d[i] *= 0.5 * complex<lf>(e[i].Y + e[fr].Y, e[fr].X - e[i].X);
		}
		fft(c), fft(d);
		vector<ll> ans(size());
		for (int i = 0; i < size(); ++i)
		{
			ll fr = (size() - i) & (size() - 1),
			   p = c[fr].X / size() + 0.5,
			   o = c[fr].Y / size() + 0.5,
			   x = d[fr].X / size() + 0.5,
			   u = d[fr].Y / size() + 0.5;
			ans[i] = (p % M + ((o + x) % M << 15) + (u % M << 30)) % M;
		}
		return ans;
	}
};

快速数论变换

原理和 FFT 相同,解决特殊情况下 FFT 的浮点误差,并且可以在同余系进行变换。

对于形如$m=2^nc+1$的费马素数,记其原根为 g,则旋转因子为$g^{(m-1)/n}$,满足$g^{m-1}=1$且$2^n\mid m-1$。

使用示例

struct FNTT : Rader, Mod
{
	vector<ll> w;
	FNTT(int N, ll M, ll G) : Rader(N), Mod(M), w(size(), pow(G, (M - 1) / size()))
	{
		for (int i = w[0] = 1; i < size(); ++i)
			w[i] = mul(w[i], w[i - 1]);
	}
	void fntt(vector<ll> &x) const
	{
		for (int i = 0; i < size(); ++i)
			if (i < at(i))
				std::swap(x[i], x[at(i)]);
		for (int i = 1, j; i < size(); i <<= 1)
			for (int j = 0; j < i; ++j)
				for (int k = j; k < size(); k += i << 1)
				{
					ll t = mul(w[size() / (i << 1) * j], x[k + i]);
					qadd(x[k + i] = x[k], M - t), qadd(x[k], t);
				}
	}
	vector<ll> ask(vector<ll> a, vector<ll> b) const
	{
		a.resize(size()), b.resize(size()), fntt(a), fntt(b);
		for (int i = 0; i < size(); ++i)
			a[i] = mul(a[i], b[i]);
		fntt(a), reverse(a.begin() + 1, a.end());
		ll u = inv(size());
		for (int i = 0; i < size(); ++i)
			a[i] = mul(a[i], u);
		return a;
	}
};

常见素数的原根

r * 2 ^ k + 1 r k g,mod(r * 2 ^ k + 1)的原根
3 1 1 2
5 1 2 2
17 1 4 3
97 3 5 5
193 3 6 5
257 1 8 3
7681 15 9 17
12289 3 12 11
40961 5 13 3
65537 1 16 3
786433 3 18 10
5767169 11 19 3
7340033 7 20 3
23068673 11 21 3
104857601 25 22 3
167772161 5 25 3
469762049 7 26 3
998244353 119 23 3
1004535809 479 21 3
2013265921 15 27 31
2281701377 17 27 3
3221225473 3 30 5
75161927681 35 31 3
77309411329 9 33 7
206158430209 3 36 22
2061584302081 15 37 7
2748779069441 5 39 3
6597069766657 3 41 5
39582418599937 9 42 5
79164837199873 9 43 5
263882790666241 15 44 7
1231453023109121 35 45 3
1337006139375617 19 46 3
3799912185593857 27 47 5
4222124650659841 15 48 19
7881299347898369 7 50 6
31525197391593473 7 52 3
180143985094819841 5 55 6
1945555039024054273 27 56 5
4179340454199820289 29 57 3

快速沃尔什变换

如果要在同余系中进行运算,则下面代码需要修改。

void fwt(vector<ll> &x, void f(ll &l, ll &r))
{
	for (int i = 1; i < x.size(); i <<= 1)
		for (int j = 0; j < i; ++j)
			for (int k = j; k < x.size(); k += i << 1)
				f(x[k], x[k + i]);
}
void fwt(ll *b, ll *e, void f(ll &l, ll &r)) //再给一个递归二分的代码便于理解
{
	if (e - b < 2)
		return;
	ll *m = b + (e - b) / 2;
	fwt(b, m, f), fwt(m, e, f);
	while (m < e)
		f(*(b++), *(m++));
}

AND

void tf(ll &l, ll &r) { l += r; }
void utf(ll &l, ll &r) { l -= r; }

OR

void tf(ll &l, ll &r) { r += l; }
void utf(ll &l, ll &r) { r -= l; }

XOR

void tf(ll &l, ll &r)
{
	ll tl = l + r, tr = l - r;
	l = tl, r = tr;
}
void utf(ll &l, ll &r) { tf(l, r), l >>= 1, r >>= 1; }

XNOR、NAND、NOR

直接用异或运算、与运算、或运算的方法求出来,然后将互反的两位交换即可。

Pell 方程

形如$x^2-Dy^2=1$(D 为任意正整数)的方程称为佩尔方程,必有最小正整数解$(x_0,y_0)$,用

$x_n=x_0x_{n-1}+Dy_0y_{n-1},y_n=y_0x_{n-1}+x_0y_{n-1}$

可递推方程的第 n 小整数解(可用矩阵快速幂求),同时还有

$2x_0x_n=x_{n-1}+x_{n+1},2x_0y_n=y_{n-1}+y_{n+1}$

Bertrand 猜想

$\forall n>3,\exist n<p<n\times 2$其中 n 为整数,p 为质数。

威尔逊定理

当且仅当 p 为素数时:$(p-1)!\equiv -1\pmod p$,

Jacobi’s Four Square Theorem

设$a^2+b^2+c^2+d^2=n$的自然数解个数为$r4(n)$,$d(n)$为$n$的约数和,由 Jacobi’s Four Square Theorem 可知,若$n$是奇数,则$r4(n)=8d(n)$,否则$r4(n)=24d(k)$,$k$是$n$去除所有$2$后的结果。

组合数学

组合数取模

使用示例

为方便,记$C(n,m)=C_n^m=\binom{n}{m}$。

struct Factorial : Mod
{
	vector<ll> fac, ifac;
	Factorial(int N, ll M) : fac(N, 1), ifac(N, 1), Mod(M)
	{
		for (int i = 2; i < N; ++i)
			fac[i] = mul(fac[i - 1], i), ifac[i] = mul(M - M / i, ifac[M % i]);
		for (int i = 2; i < N; ++i)
			ifac[i] = mul(ifac[i], ifac[i - 1]);
	}
	ll c(int n, int m) { return mul(mul(fac[n], ifac[m]), ifac[n - m]); }
	ll lucas(ll n, ll m) //卢卡斯定理求C(n,m)%M,适用于模数M小于N的情况,或者m较小的时候也可以暴力求
	{
		if (!m)
			return 1;
		if (n < m || n % M < m % M)
			return 0;
		if (n < M && m < M)
			return c(n, m);
		return mul(lucas(n / M, m / M), lucas(n % M, m % M));
	}
};

组合数 LCM

$(n + 1)lcm(C(n,0),C(n,1),\dots,C(n,k))=lcm(n+1,n,n−1,n−2,\dots,n−k+1)$

区间 lcm 的维护:对于一个数,将其分解质因数,若有因子$p^k$,那么拆分出 k 个数 $p^1,p^2,\dots,p^k$,权值都为 p,那么区间$[l,r]$内所有数的 lcm 的答案=所有在该区间中出现过的数的权值之积,可持久化线段树维护之。

Stirling 数

第一类斯特林数

第一类斯特林数$S(p,k)$的一个的组合学解释是:将 p 个物体排成 k 个非空循环排列的方法数。

递推公式:$S(p,k)=(p−1)S(p−1,k)+S(p−1,k−1),1\leq k\leq p−1;S(p,0)=0,p\ge1;S(p,p)=1,p\ge0$

第二类斯特林数

第二类斯特林数$S(p,k)$的一个的组合学解释是:将 p 个物体划分成 k 个非空不可辨别的(可以理解为盒子没有编号)集合的方法数。

递推公式:$S(p,k)=kS(p−1,k)+S(p−1,k−1),1\leq k\leq p−1;S(p,0)=0,p\ge 1;S(p,p)=1,p\ge0$

卷积形式:$S(n,m)=\frac{1}{m!}\sum_{k=0}^m(-1)^kC(m,k)(m-k)^n=\sum_{k=0}^m\frac{(-1)^k}{k!}\frac{(m-k)^n}{(m-k)!}$

同时有转化:$x^k=\sum_{i=1}^ki!C(x,i)S(k,i)$

斯特林近似公式

$n!\approx\sqrt{2\pi n}(\frac{n}{e})^n$

小球入盒模型通解

k 个球 m 个盒子 空盒子 方案数
各不相同 各不相同 允许 $m^k$
各不相同 各不相同 $m!Stirling2(k,m)$
各不相同 完全相同 允许 $\sum_{i=1}^mStirling2(k,i)$
各不相同 完全相同 $Stirling2(k,m)$
完全相同 各不相同 允许 $C(m+k−1,k)$
完全相同 各不相同 $C(k−1,m−1)$
完全相同 完全相同 允许 $\frac{1}{(1−x)(1−x^2)…(1−x^m)}的x^k项的系数$
完全相同 完全相同 $\frac{x^m}{(1−x)(1−x^2)…(1−x^m)}的x^k项的系数$

置换

使用示例

struct Permutation : vector<int>
{
	Permutation(int n = 0) : vector<int>(n) {}
	friend Permutation operator*(const Permutation &f, const Permutation &g)
	{
		Permutation ans(f.size());
		for (int i = 0; i < f.size(); ++i)
			ans[i] = g[f[i]];
		return ans;
	}
	friend Permutation inv(const Permutation &f)
	{
		Permutation ans(f.size());
		for (int i = 0; i < f.size(); ++i)
			ans[f[i]] = i;
		return ans;
	}
	friend vector<vector<int>> cycle(const Permutation &f)
	{
		vector<int> vis(f.size(), 0);
		vector<vector<int>> ans;
		for (int i = 0; i < f.size(); ++i)
			if (!vis[i])
			{
				ans.push_back(vector<int>());
				for (int j = i; !vis[j]; j = f[j])
					vis[j] = 1, ans.back().push_back(j);
			}
		return ans;
	}
};

生成字典序

下一排列

对给定的排列$a_1a_2\dots a_n$,找到$a_j$使得$a_j<a_{j+1},a_{j+1}>a_{j+2}>\dots>a_n$即这列数中最后一个相邻递增数对,然后把$a_{j+1},a_{j+2},\dots,a_n$中大于$a_j$的最小数放到位置 j,然后$a_j\dots a_n$中剩余的数从小到大排序放到$[j+1,n]$中。

bool nextPermutation(ll *b, ll *e) //标准库有这个函数next_permutation
{
	ll *i = e - 1, *j = e - 2;
	while (j >= b && *j >= *(j + 1))
		--j;
	if (j < b)
		return 0;
	while (*i <= *j)
		--i;
	return swap(*i, *j), reverse(j + 1, e), 1;
}

二项式反演

$f(n)=\sum_{k=0}^nC(n,k)g(k),g(n)=\sum_{k=0}^n(−1)^{n−k}C(n,k)f(k)$

第 k 小期望

$f(n,k)$表示有 n 个变量,和为 1,第 k 小的期望。

$f(n,k)=\frac{1}{n^2}+(1-\frac{1}{n})f(n-1,k-1),f(n,0)=0$

错排数

考虑一个有 n 个元素的排列,若一个排列中所有的元素都不在自己原来的位置上,那么这样的排列就称为原排 列的一个错排。

n 个元素的错排数$D_n$满足递推公式:$D_1=0,D_2=1,D_n=(n−1)(D_{n−2}+D_{n−1})$

通项:$D(n)=n![\frac{(-1)^2}{2!}+\dots+\frac{(-1)^{n-1} }{(n-1)!}+\frac{(-1)^n}{n!}]=\lfloor\frac{n!}{e}+\frac{1}{2}\rfloor$

Bonuli 数

使用示例

$B_n = -\frac{1}{C(n+1,n)}(C(n+1,0)B_0+C(n+1,1)B_1+\dots+C(n+1,n-1)B_{n-1})=-\frac{1}{n+1}(C(n+1,0)B_0+C(n+1,1)B_1+\dots+C(n+1,n-1)B_{n-1})$

可用于计算任意正整数次数的幂和:$\sum_{i=1}^ni^k=\frac{1}{k+1}\sum_{j=0}^kC(k+1,j)B_jn^{k+1-j}$

struct Bonuli : Factorial
{
	vector<ll> b;
	Bonuli(int N, ll M) : Factorial(N, M), b(N, 0)
	{
		for (int i = b[0] = 1; i < N; ++i)
		{
			for (int j = 0; j < i; ++j)
				b[i] = qadd(b[i], mul(b[j], c(i + 1, j)));
			b[i] = qadd(M, -mul(b[i], mul(fac[i], ifac[i + 1])));
		}
	}
	ll ask(ll n, int k)
	{
		ll r = 0, w = 1, u = add(n, 1);
		for (int i = 1; i <= k + 1; ++i)
			r = qadd(r, mul(mul(b[k + 1 - i], c(k + 1, i)), w = mul(w, u)));
		return mul(r, mul(fac[k], ifac[k + 1]));
	}
};

Catalan 数

$h_1=1,h_n=\frac{4n−2}{n+1}h_{n−1}=\frac{C(2n,n)}{n+1}=C(2n,n)−C(2n,n−1)$。

在一个格点阵列中,从$(0,0)$点走到$(n,m)$点且不经过对角线$x=y$的方法数:$C(n+m−1,m)−C(n+m−1,m−1),x>y;C(n+m,m)−C(n+m,m−1),x\ge y$。

常见的 Catalan 数:括号序的个数、凸多边形三角剖分的方案数等。

Bell 数

把 n 个带标号的物品划分为若干不相交集合的方案数称为贝尔数,其递推公式:$B_n=\sum_{i=0}^{N-1}C_{n-1}^iB_i$

前几项贝尔数:

1,2,5,15,52,203,877,4140,21147,115975,678570,4213597,27644437,190899322,1382958545,...

等价类容斥

考虑容斥,Bell(p)枚举所有等价情况。对于一种情况,强制了一个等价类里面的数都要相同,其它的可以相同也可以不同。

容斥系数为:$(−1)^{p−等价类个数}(每个等价类大小−1)!之积$。

Grey 码

格雷序列第 i 个是i^(i>>1)。长为 n 的 01 序列共$2^n$个,下标从$0\dots 2^n-1$。

扩展 Cayley 公式

对于 n 个点,m 个连通块的图,假设每个连通块有 a[i]个点,那么用 s−1 条边把它连通的方案数为$n^{s−2}a[1]a[2]\dots a[m]$。

超立方体

n 维超立方体有$2^{n−i}C(n,i)$个 i 维元素。

枚举位集 I 的非空子集 J

for(J=I; J; J=I&J1) {}

线性代数

矩阵

typedef array<array<ll, N>, N> Mat;

乘法和快速幂

Mat operator*(const Mat &a, const Mat &b)
{
	Mat r;
	for (int i = 0; i < r.size(); ++i)
		for (int j = 0; j < r.size(); ++j)
			for (int k = r[i][j] = 0; k < r.size(); ++k)
				M.qadd(r[i][j], M.mul(a[i][k], b[k][j]));
	return r;
}
Mat pow(Mat a, ll b)
{
	Mat r;
	for (int i = 0; i < r.size(); ++i)
		for (int j = 0; j < r[i].size(); ++j)
			r[i][j] = i == j;
	for (; b; b >>= 1, a = a * a)
		if (b & 1)
			r = r * a;
	return r;
}

行列式

ll det(Mat a, int n)
{
	ll ans = 1;
	for (int i = 0; i < n; ++i)
	{
		for (int j = i + 1; j < n; ++j)
			while (fabs(a[j][i]) > EPS)
			{
				ll t = a[i][i] / a[j][i];
				for (int k = i; k < n; ++k)
					a[i][k] -= t * a[j][k], swap(a[i][k], a[j][k]);
			}
		if (fabs(ans *= a[i][i]) < EPS)
			return 0;
	}
	return ans;
}

高斯消元

void GaussElimination(Mat &a, int n) //a为增广矩阵,要求n*n的系数矩阵可逆,运行结束后a[i][n]为第i个未知数的值
{
	for (int i = 0, r; i < n; ++i)
	{
		for (int j = r = i; j < n; ++j)
			if (fabs(a[r][i]) < fabs(a[j][i]))
				r = j;
		if (r != i)
			swap_ranges(a[r].begin(), a[r].begin() + n + 1, a[i]);
		for (int j = n; j >= i; --j)
			for (int k = i + 1; k < n; ++k)
				a[k][j] -= a[k][i] * a[i][j] / a[i][i];
	}
	for (int i = n - 1; ~i; --i)
	{
		for (int j = i + 1; j < n; ++j)
			a[i][n] -= a[j][n] * a[i][j];
		a[i][n] /= a[i][i];
	}
}

线性基

向量线性基

add 返回要插入的向量 z 是否与已插入的线性无关。

struct Base
{
	vector<vector<double>> v;
	Base(int N) : v(N, vector<double>(N, 0)) {} //R^N的子空间
	bool add(vector<double> x)
	{
		for (int i = 0; i < x.size(); ++i)
			if (fabs(x[i]) > EPS)
			{
				if (fabs(v[i][i]) < EPS)
					return v[i] = x, 1;
				double t = x[i] / v[i][i];
				for (int j = 0; j < x.size(); ++j)
					x[j] -= t * v[i][j];
			}
		return 0;
	}
};

异或线性基

若要查询第 k 小子集异或和,则把 k 写成二进制,对于是 1 的第 i 位,把从低位到高位第 i 个不为 0 的数异或进答案。若要判断是否有非空子集的异或和为 0,如果不存在自由基,那么说明只有空集的异或值为 0,需要高斯消元来判断。

struct BaseXOR
{
	vector<ll> a;
	BaseXOR() : a(64, 0) {}
	ll ask() //查询最大子集异或和
	{
		ll t = 0;
		for (int i = a.size() - 1; ~i; --i)
			t = max(t, t ^ a[i]);
		return t;
	}
	bool add(ll x)
	{
		for (int i = a.size() - 1; ~i; --i)
			if (x >> i & 1)
			{
				if (a[i])
					x ^= a[i];
				else
					return a[i] = x, 1;
			}
		return 0;
	}
	bool check(ll x) //判断一个数是否能够被异或出,0根据需要特判
	{
		for (int i = a.size() - 1; ~i; --i)
			if (x >> i & 1)
				if (x ^= a[i], !x)
					return 1;
		return 0;
	}
};

离散数学

归并排序求逆序对

使用示例

ll merge_sort(ll *b, ll *e) //int存答案会爆
{
	if (e - b < 2)
		return 0;
	ll *m = b + (e - b) / 2, ans = merge_sort(b, m) + merge_sort(m, e);
	for (ll *i = b, *j = m; i < m && j < e; ++i)
		for (; j < e && *j < *i; ++j)
			ans += m - i;
	return inplace_merge(b, m, e), ans;
}

约瑟夫问题

ll josephus(ll n, ll k) //编号0~n-1,每k个出列,时间复杂度O(min(n,k))
{
	if (n < 3)
		return k % n;
	if (n < k)
		return (Josephus(n - 1, k) + k) % n;
	ll ret = Josephus(n - n / k, k) - n % k;
	return ret < 0 ? ret + n : ret + ret / (k - 1);
}

堆的操作

和 STL 的函数push_heappop_heap用法完全相同,实现的是大根堆。

ll pushHeap(ll *b, ll *e)
{
	ll tmp = *--e;
	for (int i = e - b, p;; b[i] = b[p], i = p)
		if (i <= 0 || b[p = (i - 1) / 2] >= tmp)
			return b[i] = tmp;
}
ll popHeap(ll *b, ll *e)
{
	ll tmp = b[0];
	for (int i = 0, siz = --e - b, ch;; b[i] = b[ch], i = ch)
	{
		if (ch = i * 2 + 1, ch + 1 < siz && b[ch] < b[ch + 1])
			++ch;
		if (i >= siz / 2 || b[siz] >= b[ch])
			return b[i] = b[siz], b[siz] = tmp;
	}
}

蔡勒公式

$w=(\lfloor\frac{c}{4}\rfloor-2c+y+\lfloor\frac{y}{4}\rfloor+\lfloor\frac{13(m+1)}{5}\rfloor+d-1)\mod7$

w:$0,1,\dots,6$对应周日,周一,$\dots$,周六

c:世纪减 1(即年份前两位数)。

y:年份后两位数。

m:月($3\leq m\leq14$,即在蔡勒公式中,1、2 月要看作上一年的 13、14 月来计算)。

d:日。

曼哈顿距离的变换

$\mid x_1−x_2\mid +\mid y_1−y_2\mid=\max (\mid (x_1 + y_1)−(x_2 + y_2)\mid ,\mid (x_1 −y_1)−(x_2 −y_2)\mid )$

皮克定理

顶点坐标均是整点(或说正方形格点)的简单多边形中,面积 S 和内部格点数目 n、边上格点数目 m 的满足关系$S=n+\frac{m}{2}-1$。

矩形面积并、矩形面积交、矩形周长并(线段树、扫描线总结)

详见这篇博文

数学分析

增长趋势

[n\to+\infty,\forall p,q>0,a>1,{(\ln n)}^q\ll n^p\ll a^n\ll n!\ll n^n]

积分表

反读可得导数表,此处略。

[\int k\,\mathrm{d}x=kx+C
\int x^a\,dx=\frac{x^{a+1}}{a+1}+C
\int\frac{1}{x}\,dx=\ln\mid x\mid +C
\int e^x\,dx=e^x + C
\int a^x\,dx=\frac{a^x}{\ln a}+C
\int\cos x\,dx=\sin x+C
\int\sin x\,dx=-\cos x+C
\int\frac{1}{cos^2x}\,dx=\int\sec^2 x\,dx=\tan x+C
\int\frac{1}{sin^2x}\,dx=\int\csc^2 x\,dx=-\cot x+C
\int\frac{1}{\sqrt{1-x^2}}\,dx=\arcsin x+C=-\arccos x+C
\int\frac{1}{1+x^2}\,dx=\arctan x+C=-arccot\,x+C
\int\sec x\tan x\,dx=\sec x+C
\int\csc x\cot x\,dx=-\csc x+C
\int\tan x\,dx=-\ln\mid \cos x\mid +C
\int\cot x\,dx=\ln\mid \sin x\mid +C
\int\sec x\,dx=\ln\mid \sec x+\tan x\mid +C
\int\csc x\,dx=\ln\mid \csc x-\cot x\mid +C
\int sh\,x\,dx=ch\,x+C
\int ch\,x\,dx=sh\,x+C
\int\frac{1}{x^2+a^2}\,dx=\frac{1}{a}\arctan\frac{x}{a}+C
\int\frac{1}{x^2-a^2}\,dx=\frac{1}{2a}\ln\mid \frac{x-a}{x+a}\mid +C
\int\frac{1}{\sqrt{a^2-x^2}}\,dx=\arcsin\frac{x}{a}+C
\int\frac{1}{\sqrt{x^2-a^2}}\,dx=\ln\mid x+\sqrt{x^2-a^2}\mid +C
\int\frac{1}{\sqrt{x^2+a^2}}\,dx=\ln\mid x+\sqrt{x^2+a^2}\mid +C\]

积分求几何量

弧长

若简单闭曲线

[\begin{cases} x=x(t),
y=y(t), \end{cases} t\in[\alpha,\beta]]

端点处重合($x(\alpha)=x(\beta),y(\alpha)=y(\beta)$)且其他地方不自交,$x(t),y(t)$连续且满足

[[x’(t)]^2+[y’(t)]^2\ne0,\forall t\in[\alpha,\beta]]

此时称曲线光滑,其长度

[s=\int_\alpha^\beta\sqrt{[x’(t)]^2+[y’(t)]^2}\,dt]

此式可对称推广到高维空间曲线。 极坐标下,

[r=r(\theta),\theta\in[\alpha,\beta]]

的长度为

[s=\int_\alpha^\beta\sqrt{[r(\theta)]^2+[r’(\theta)]^2}\,d\theta]

面积

若简单闭曲线

[\begin{cases} x=x(t),
y=y(t), \end{cases} t\in[\alpha,\beta]]

端点处连续($x(\alpha)=x(\beta),y(\alpha)=y(\beta)$)且其他地方不自交,$x(t),y(t)$都逐段有连续微商,则此闭合曲线围起来的有界区域面积

[S=-\int_\alpha^\beta x’(t)y(t)\,dt=-\int_\alpha^\beta y(t)\,dx(t)=-\oint_\Gamma y\,dx=\oint_\Gamma x\,dy]

等式右边称为曲线$\Gamma$上的积分,其计算方法是带入参数方程到定积分计算式中,积分上下限为始点与终点对应的参数值。下限并不总是小于上限,参数从下限到上限变化时对应曲线的正向(沿正向观察时,曲线所围的区域永远在左侧)。 极坐标下,连续非负曲线$r=r(\theta)$与向径$\theta=\alpha,\theta=\beta$,其中$0\leq\beta-\alpha\leq2\pi$所围成的平面图形面积

[S=\frac{1}{2}\int_\alpha^\beta r^2(\theta)\,d\theta]

体积

记立体过 x 点且垂直于 x 轴的截面面积为$S(x)$,则其体积

[V=\int_a^bS(x)\,dx]

连续曲线$y=f(x)\ge 0,x\in[a,b]$绕 x 轴旋转一周产生的旋转体体积

[V=\pi\int_a^by^2\,dx]

旋转体侧面积

若曲线由参数方程

[\begin{cases} x=x(t),
y=y(t), \end{cases} t\in[\alpha,\beta]]

给出,则其绕 x 轴旋转体的侧面积

[s=2\pi\int_\alpha^\beta y(t)\sqrt{[x’(t)]^2+[y’(t)]^2}\,dt]

方向导数

设三元函数$u=f(x,y,z)$在点$P_0(x_0,y_0,z_0)$的某邻域内有定义,任意给定始于点$P_0$的射线$l$,$P(x,y,z)$为 l 上且含于定义域内的点。若极限

[\lim_{r(p,p_0)\to0^+}\frac{f(P)-f(P_0)}{r(P,P_0)}=\lim_{r(p,p_0)\to0^+}\frac{\Delta_lf(P_0)}{r(P,P_0)}]

存在,则称该极限值为函数$f$在点$P_0$沿方向$l$的方向导数,记为 $\frac{\partial f}{\partial l}\mid _{P_0}$或$\frac{\partial f(P_0)}{\partial l}$,$\frac{\Delta_lf(P_0)}{r(P,P_0)}$称为函数在$P_0$点沿$l$方向的增量。

特别地,$\frac{\partial f(P_0)}{\partial x}$就是函数在$P_0$点沿$x$轴正向的方向导数,$y,z$轴上的方向导数同理。若函数在$P_0$点可微,则其在$P_0$沿任何方向$l$的方向导数都存在,则有以下公式

[\frac{\partial f(P_0)}{\partial l}=(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\frac{\partial f}{\partial z})\mid _{P_0}\cdot\vec{l_0}]

其中$\vec{l_0}=(\cos\alpha,\cos\beta,cos\gamma)=\frac{1}{\rho}(\Delta x,\Delta y,\Delta z)$为$l$的方向余弦。

曲率

若曲线由参数方程

[\begin{cases} x=x(t),
y=y(t), \end{cases} t\in[\alpha,\beta]]

给出且有二阶微商,则其在一点的曲率

[K=\frac{\mid y’‘x’-y’x’‘\mid }{[x’^2+y’^2]^{\frac{3}{2}}}]

若$y=f(x)$,则

[K=\frac{\mid y’‘\mid }{(1+y’^2)^\frac{3}{2}}]

同时记$\frac{1}{K}$为曲率半径。

空间曲线的切线与法平面

若已知曲线上一点$P(x_0,y_0,z_0)$处的切向量为$\tau(x_0,y_0,z_0)=(A,B,C)$则曲线在该点的切线方程为

[\frac{x-x_0}A=\frac{y-y_0}B=\frac{z-z_0}C]

法平面方程为

[A(x-x_0)+B(y-y_0)+C(z-z_0)=0]

当曲线由参数方程

[\begin{cases} x=x(t),
y=y(t),
z=z(t), \end{cases} t\in[\alpha,\beta]]

给出时,曲线在 P 点的切向量为

[\tau=\pm(x’(t_0),y’(t_0),z’(t_0))]

更一般地,若曲线用两曲面的交线给出

[\begin{cases} F(x,y,z)=0,
G(x,y,z)=0, \end{cases}]

且在 P 点的某邻域能确定函数组$y=y(x),z=z(x)$满足$y_0=y(x_0),z_0=z(x_0)$,且$y’(x),z’(x)$存在,则曲线在 P 点的切向量

[\tau=\pm(\frac{\partial(F,G)}{\partial(y,z)},\frac{\partial(F,G)}{\partial(z,x)},\frac{\partial(F,G)}{\partial(x,y)})]

空间曲面的切平面与法线

若已知曲面上一点$P(x_0,y_0,z_0)$处的切平面的法向量为$\vec n=(A’,B’,C’)$则曲线在该点的法线方程为

[\frac{x-x_0}{A’}=\frac{y-y_0}{B’}=\frac{z-z_0}{C’}]

切平面方程为

[A’(x-x_0)+B’(y-y_0)+C’(z-z_0)=0]

当曲面方程为$\pi:F(x,y,z)=0$在曲面上任取一条过 P 的曲线,设其方程为

[\begin{cases} x=x(t),
y=y(t),
z=z(t), \end{cases} t\in[\alpha,\beta]]

此时有$F(x(t),y(t),z(t))=0$令$t=t_0$两边对 t 求导,并写成向量的内积式,得

[(F_x,F_y,F_z)_P\cdot(x’(t_0),y’(t_0),z’(t_0))=0]

则曲线在 P 点的法向量为

[\vec{n}=\pm(F_x,F_y,F_z)_P]

若曲线由参数方程给出

[\begin{cases} x=x(u,v),
y=y(u,v),
z=z(u,v), \end{cases}]

则曲线在 P 点的法向量

[\vec{n}=\pm(\frac{\partial(y,z)}{\partial(u,v)},\frac{\partial(z,x)}{\partial(u,v)},\frac{\partial(x,y)}{\partial(u,v)})]

高阶导数与泰勒公式

用$f^{(n)}(x)$表示 f(x)的 n 阶导数,只要让余项<EPS即可计算指定函数到任意精确度,特别取 a=0 时称为麦克劳林公式。

[f(x)=f(a)+f^{(1)}(a)(x-a)+\frac{f^{(2)}(a)}{2!}(x-a)^2+\dots+\frac{f^{(n)}(a)}{n!}(x-a)^n+R_n(x)]

佩亚诺余项

[R_n(x)=o((x-a)^n)]

积分余项

[R_n(x)=\frac{1}{n!}\int_a^x(x-t)^nf^{(n+1)}(t)\,dt]

拉格朗日余项

[R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}(x-a)^{n+1},a<\xi<x]

柯西余项

[R_n(x)=\frac{(x-a)^{n+1}}{n!}(1-\theta)^nf^{(n+1)}(a+\theta(x-a)),0<\theta<1]

对数函数

[[\ln(1+x)]^{(n)}=(-1)^{n-1}(n-1)!(1+x)^{-n} \ln(1+x)=x-\frac{x^2}{2}+\frac{x^3}{3}-\frac{x^4}{4}+\dots+(-1)^{n-1}\frac{x^n}{n}+R_n(x)]

幂函数

[[(1+x)^a]^{(n)}=a(a-1)\dots(a-n+1)(1+x)^{a-n}
(1+x)^a=1+ax+\frac{a(a-1)}{2!}x^2+\dots+\frac{a(a-1)\dots(a-n+1)}{n!}x^n+R_n(x)]

三角函数

[(\sin x)^{(n)}=\sin(x+\frac{n\pi}{2})
\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\dots+(-1)^{k-1}\frac{x^{2k-1}}{(2k-1)!}+R{2k}(x) \R{2k}(x)=(-1)^k\frac{\cos\theta x}{(2k+1)!}x^{2k+1}
(\cos x)^{(n)}=\cos(x+\frac{n\pi}{2})
\cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\dots+(-1)^{k-1}\frac{x^{2k-2}}{(2k-2)!}+R{2k-1}(x)
R
{2k-1}(x)=(-1)^k\frac{\cos\theta x}{(2k)!}x^{2k}]

指数函数

[(e^x)^{(n)}=e^x
e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\dots+\frac{x^n}{n!}+R_n(x)
R_n(x)=\frac{e^{\theta x}}{(n+1)!}x^{n+1},\xi=\theta x,0<\theta<1]

二元函数

设$f(x,y)$在$P_0(x_0,y_0)$的某邻域$O(P_0)$内有直到$n+1$阶连续偏导数,则对$O(P_0)$内$\forall(x_0+\Delta x,y_0+\Delta y),\exists\theta\in(0,1)$,使得

[f(x0+\Delta x,y_0+\Delta y)=\sum{k=0}^n\frac{1}{k!}(\frac{\partial}{\partial x}\Delta x+\frac{\partial}{\partial y}\Delta y)^kf(x_0,y_0)+R_n]

其中

[R_n=\frac{1}{(n+1)!}(\frac{\partial}{\partial x}\Delta x+\frac{\partial}{\partial y}\Delta y)^{n+1}f(x_0+\theta\Delta x,y_0+\theta\Delta y)]

级数部分和

幂级数

快速计算幂级数的部分和$\sum_{i=1}^ni^k\mod M$可借助伯努利数,详见组合数学模板。

[\sum_{i=1}^ni^1=\frac 1 2n(n+1)
\sum_{i=1}^ni^2=\frac 1 6n(n+1)(2n+1)
\sum_{i=1}^ni^3=\frac 1 4[n(n+1)]^2
\sum_{i=1}^ni^4=\frac 1{30}n(n+1)(2n+1)(3n^2+3n-1)
\sum_{i=1}^ni^5=\frac 1{12}[n(n+1)]^2(2n^2+2n-1)
\sum_{i=1}^ni^6=\frac 1{42}n(n+1)(2n+1)(3n^4+6n^3-3n+1)]

调和级数

[n\to\infty,\sum_{i=1}^n\frac 1 i\to\ln n+r,r\approx0.5772156649015328\dots]

二分求零点、三分求极值点

需要$f(x)$在区间$[l,r]$上单调/凹凸性唯一。

lf bs(lf l, lf r, lf f(lf x))
{
	if (r - l < EPS)
		return l;
	lf m = (l + r) / 2;
	return sgn(f(l) * f(m)) < 0 ? bs(l, m, f) : ts(m, r, f);
}
lf ts(lf l, lf r, lf f(lf x))
{
	if (r - l < EPS)
		return l;
	lf d = (r - l) / 3, lm = l + d, rm = r - d;
	return f(lm) < f(rm) ? ts(l, rm, f) : ts(lm, r, f); //极小值
}

自适应辛普森求积分

使用示例

这篇论文论证了加一个十五分之一的偏移收敛会比较快…

struct Simpson
{
	lf simpson(lf a, lf b, lf f(lf x))
	{
		return (f(a) + 4 * f((a + b) / 2) + f(b)) * (b - a) / 6;
	}
	lf ask(lf a, lf b, lf f(lf x), lf e = EPS)
	{
		lf c = (a + b) / 2, L = simpson(a, c, f), R = simpson(c, b, f), delta = (L + R - simpson(a, b, f)) / 15;
		return fabs(delta) < e ? L + R + delta : ask(a, c, f, e / 2) + ask(c, b, f, e / 2);
	}
};

插值法

拉格朗日插值法:插值多项式和插值基函数的形式对称,容易编程。但是,增加节点时,需要重新计算每一个插值基函数。要在$\pmod p$意义下进行的话,那么 p 只能是质数。 牛顿插值法:当插值节点增加时,之前已计算的结果仍然能用,每增加一个节点,只要再增加一项即可,从而避免了重复性计算。如果要 mod 非质数的话,那么就要用牛顿插值法。

typedef complex<lf> Coord;
#define X real()
#define Y imag()
struct Lagrange
{
	lf ask(const vector<Coord> &p, lf x) //返回p确定的多项式函数在x处的值
	{
		lf ret = 0;
		for (int i = 0; i < p.size(); ++i)
		{
			lf tmp = p[i].Y;
			for (int j = 0; j < p.size(); ++j)
				if (i != j)
					tmp *= (x - p[j].X) / (p[i].X - p[j].X);
			ret += tmp;
		}
		return ret;
	}
	vector<lf> ask(vector<Coord> p) //返回p确定的多项式系数向量
	{
		vector<lf> ret(p.size()), sum(p.size());
		ret[0] = p[0].Y, sum[0] = 1;
		for (int i = 1; i < p.size(); ++i)
		{
			for (int j = p.size() - 1; j >= i; --j)
				p[j].Y = (p[j].Y - p[j - 1].Y) / (p[j].X - p[j - i].X);
			for (int j = i; ~j; --j)
				sum[j] = (j ? sum[j - 1] : 0) - sum[j] * p[i - 1].X,
				ret[j] += sum[j] * p[i].Y;
		}
		return ret;
	}
};
struct Newton
{
	lf differenceQuotient(const vector<Coord> &p, int k) //计算差商
	{
		lf ret = 0;
		for (int i = 0; i <= k; ++i)
		{
			lf tmp = p[i].Y;
			for (int j = 0; j <= k; ++j)
				if (i != j)
					tmp /= p[i].X - p[j].X;
			ret += tmp;
		}
		return ret;
	}
	lf ask(const vector<Coord> &p, lf x)
	{
		lf ret = p[0].Y;
		for (int i = 1; i < p.size(); ++i)
		{
			lf tmp = differenceQuotient(p, i); //多次求,可O(n^3)预处理优化
			for (int j = 0; j < i; ++j)
				tmp *= x - p[j].X;
			ret += tmp;
		}
		return ret;
	}
};

[]

计算几何

二维

点和向量化为坐标 Coord 进行运算,使用 stl 中的 complex 实现。

复数相乘的几何意义为长度相乘,极角相加。

用直线上的一点 p 和方向向量 v 表示一条经过 p 的直线,直线上的所有点 q 满足 q=p+t*v,其中 t 是参数;当限制 t≥0 时,该参数方程表示射线;限制 0≤t≤1 时,该参数方程表示线段。

此外,如果已知线段端点 a1 和 a2,可以通过 Line(a1,a2-a1)来得到对应的参数形式。

Morley 定理:三角形每个内角的三等分线相交成等边三角形。

欧拉定理:平面图的点数 V、边数 E 和面数 F 满足 V+F-E=2。

typedef double lf;
typedef complex<lf> Coord;
const lf EPS = 1e-9;
#define X real()
#define Y imag()
struct Line
{
	Coord p, v;
	Line(Coord p = Coord(), Coord v = Coord()) : p(p), v(v) {}
	Coord point(lf t) const { return p + v * t; }
};
struct Circle
{
	Coord c;
	lf r;
	Circle(Coord c = Coord(), lf r = 0) : c(c), r(r) {}
	Coord point(lf t) const { return c + polar(r, t); } //t为参数,幅角
};
/*
Coord(lf x=0,lf y=0);//构造函数
lf real(Coord a);//a的实部(复平面的横坐标),也可写作a.real()
lf imag(Coord a);//a的虚部(复平面的纵坐标),也可写作a.imag()

lf abs(Coord a);//向量a的模长,或是点a到原点的距离
lf norm(Coord a);//abs的平方,比abs快,但是要注意浮点数精度溢出

lf arg(Coord a);//a的幅角,与atan2(a.real(),a.imag())等价
Coord polar(lf r,lf t);//极坐标生成方式,r为幅值,t为幅角

//运算符重载+、-、*、/(以及对应的赋值运算,但是赋值运算不能写在表达式中,详见参考地址)、<<、>>(输出括号形式的坐标)
*/
int sgn(lf d) { return (d > EPS) - (d < -EPS); }

bool operator!=(const Coord &A, const Coord &B) { return sgn(A.X - B.X) || sgn(A.Y - B.Y); } //不等运算符,涉及到浮点数比较要重写

bool operator==(const Coord &A, const Coord &B) { return !(A != B); }

bool cmpCoord(const Coord &A, const Coord &B) { return sgn(A.X - B.X) ? sgn(A.X - B.X) < 0 : sgn(A.Y - B.Y) < 0; } //复数没有小于运算,只能这样定义一个比较函数

bool cmpLine(const Line &A, const Line &B) { return sgn(arg(A.v) - arg(B.v)) < 0; } //按极角排序,求凸包中使用

lf Dot(const Coord &A, const Coord &B) { return A.X * B.X + A.Y * B.Y; }

lf Cross(const Coord &A, const Coord &B) { return A.X * B.Y - B.X * A.Y; }

lf Angle(const Coord &A, const Coord &B) { return acos(Dot(A, B) / abs(A) / abs(B)); }

lf Area2(const Coord &A, const Coord &B, const Coord &C) { return Cross(B - A, C - A); } //三角形ABC有向面积的两倍

Coord Rotate(const Coord &A, lf rad) { return A * polar(1.0, rad);} //向量A逆时针旋转rad弧度

Coord Normal(const Coord &A) //A的法向量,把A逆时针旋转九十度并长度化为1
{
	lf L = abs(A);
	return Coord(-A.Y / L, A.X / L);
}

bool onLeft(const Coord &P, const Line &L) { return sgn(Cross(L.v, P - L.p)) > 0; } //p是否在有向直线L左侧,不含线上

lf DistanceToLine(const Coord &P, const Line &L) { return Cross(L.v, P - L.p) / abs(L.v); } //点到直线距离(有向)

lf DistanceToLine(const Coord &P, const Coord &A, const Coord &B) { return DistanceToLine(P, Line(A, B - A)); }

lf DistanceToSegment(const Coord &P, const Coord &A, const Coord &B) //点到线段的距离(无向)
{
	if (A == B)
		return abs(P - A);
	Coord v1 = B - A, v2 = P - A, v3 = P - B;
	if (sgn(Dot(v1, v2)) < 0)
		return abs(v2);
	if (sgn(Dot(v1, v3)) > 0)
		return abs(v3);
	return fabs(DistanceToLine(P, Line(A, B - A)));
}

Coord getLineProjection(const Coord &P, const Line &L) { return L.point(Dot(L.v, P - L.p) / norm(L.v)); } //点在直线上的投影

Coord getLineProjection(const Coord &P, const Coord &A, const Coord &B) { return getLineProjection(P, Line(A, B - A)); }

Coord getSymmetry(const Coord &P, const Coord &O) { return O + O - P; } //P关于O的对称点

Coord getSymmetry(const Coord &P, const Line &L) { return getSymmetry(P, getLineProjection(P, L)); } //P关于L的对称点

Coord getLineIntersection(const Line &L1, const Line &L2) { return L1.point(Cross(L2.v, L1.p - L2.p) / Cross(L1.v, L2.v)); } //直线交点,须确保两直线相交

Coord getLineIntersection(const Coord &A1, const Coord &A2, const Coord &B1, const Coord &B2) { return getLineIntersection(Line(A1, A2 - A1), Line(B1, B2 - B1)); }

bool SegmentProperIntersection(const Coord &A1, const Coord &A2, const Coord &B1, const Coord &B2) //线段相交判定,交点不在一条线段的端点
{
	lf C1 = Cross(A2 - A1, B1 - A1), C2 = Cross(A2 - A1, B2 - A1),
	   C3 = Cross(B2 - B1, A1 - B1), C4 = Cross(B2 - B1, A2 - B1);
	return sgn(C1) * sgn(C2) < 0 && sgn(C3) * sgn(C4) < 0;
}

bool onSegment(const Coord &P, const Coord &A1, const Coord &A2) { return sgn(Dot(A1 - P, A2 - P)) < 0 && !sgn(Cross(A1 - P, A2 - P)); } //判断点是否在线段上,不包含端点

lf PolygonArea(const vector<Coord> &p) //计算多边形的有向面积,凸多边形即为面积
{
	lf s = 0;
	for (int i = 2; i < p.size(); ++i)
		s += Area2(p[0], p[i - 1], p[i]);
	return s / 2;
}

int inPolygon(const Coord &p, const vector<Coord> &poly) //点在多边形内的判定,转角法,正值为内部,0为外部,-1在边界上
{
	int ans = 0;
	for (int i = 0, k, d1, d2, n = poly.size(); i != n; ++i)
	{
		if (onSegment(p, poly[i], poly[(i + 1) % n]))
			return -1; //在边界上
		k = sgn(Cross(poly[(i + 1) % n] - poly[i], p - poly[i]));
		d1 = sgn(poly[i].Y - p.Y);
		d2 = sgn(poly[(i + 1) % n].Y - p.Y);
		if (k > 0 && d1 <= 0 && d2 > 0)
			++ans;
		if (k < 0 && d2 <= 0 && d1 > 0)
			--ans;
	}
	return ans;
}

vector<Coord> ConvexHull(vector<Coord> p, int collineation = 1) //获得凸包,不希望凸包的边上有输入点第二个参数传0
{
	vector<Coord> ans;
	sort(p.begin(), p.end(), cmpCoord); //先比横坐标再比纵坐标
	for (int i = 0; i < p.size(); ++i)  //求出下凸包
	{
		while (ans.size() > 1 && sgn(Area2(ans[ans.size() - 2], ans[ans.size() - 1], p[i])) < collineation)
			ans.pop_back();
		ans.push_back(p[i]);
	}
	for (int i = p.size() - 2, k = ans.size(); i >= 0; --i) //求出上凸包
	{
		while (ans.size() > k && -sgn(Area2(ans[ans.size() - 1], ans[ans.size() - 2], p[i])) < collineation)
			ans.pop_back();
		ans.push_back(p[i]);
	}
	if (p.size() > 1)
		ans.pop_back();
	return ans;
}

vector<Coord> cutPolygon(const vector<Coord> &poly, const Coord &A, const Coord &B) //用有向直线A->B切割多边形poly, 返回「左侧」。 如果退化,可能会返回一个单点或者线段,复杂度O(n^2)
{
	vector<Coord> newpoly;
	for (int i = 0, n = poly.size(); i != n; ++i)
	{
		Coord C = poly[i], D = poly[(i + 1) % n];
		if (sgn(Cross(B - A, C - A)) >= 0)
			newpoly.push_back(C);
		if (!sgn(Cross(B - A, C - D)))
		{
			Coord ip = getLineIntersection(Line(A, B - A), Line(C, D - C));
			if (onSegment(ip, C, D))
				newpoly.push_back(ip);
		}
	}
	return newpoly;
}

vector<Coord> getHalfPlaneIntersection(vector<Line> L) //半平面交
{
	sort(L.begin(), L.end(), cmpLine);  //按极角排序
	vector<Coord> p(L.size(), Coord()); //p[i]为q[i]和q[i+1]的交点
	int first = 0, last = 0;			//双端队列的第一个元素和最后一个元素
	vector<Line> q(L.size(), Line());   //双端队列
	q[0] = L[0];						//队列初始化为只有一个半平面L[0]

	for (int i = 0, n = L.size(); i != n; ++i)
	{
		while (first < last && !onLeft(p[last - 1], L[i]))
			--last;
		while (first < last && !onLeft(p[first], L[i]))
			++first;

		q[++last] = L[i];
		if (!sgn(Cross(q[last].v, q[last - 1].v)))
		{
			--last;
			if (onLeft(L[i].p, q[last]))
				q[last] = L[i];
		}
		if (first < last)
			p[last - 1] = getLineIntersection(q[last - 1], q[last]);
	}
	while (first < last && !onLeft(p[last - 1], q[first]))
		--last; //删除无用平面

	if (last - first <= 1)
		return vector<Coord>(); //空集
	p[last] = getLineIntersection(q[last], q[first]);
	return vector<Coord>(p.begin() + first, p.begin() + last + 1); //从deque复制到输出中
}

int getLineCircleIntersection(const Line &L, const Circle &C, vector<Coord> &sol)
{
	lf a = L.v.X,
	   b = L.p.X - C.c.X,
	   c = L.v.Y,
	   d = L.p.Y - C.c.Y,
	   e = a * a + c * c,
	   f = 2 * (a * b + c * d),
	   g = b * b + d * d - C.r * C.r,
	   delta = f * f - 4 * e * g;

	if (sgn(delta) < 0)
		return 0;
	if (!sgn(delta))
		return sol.push_back(L.point(-f / (2 * e))), 1;
	sol.push_back(L.point((-f - sqrt(delta)) / (2 * e)));
	sol.push_back(L.point((-f + sqrt(delta)) / (2 * e)));
	return 2;
}

int getCircleIntersection(const Circle &C1, const Circle &C2, vector<Coord> &sol)
{
	lf d = abs(C1.c - C2.c);

	if (!sgn(d))
		return sgn(C1.r - C2.r) ? 0 : -1; //重合返回-1

	if (sgn(C1.r + C2.r - d) < 0 || sgn(fabs(C1.r - C2.r) - d) > 0) //外离或内含
		return 0;

	lf a = arg(C2.c - C1.c),
	   da = acos((C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d));

	Coord p1 = C1.point(a - da), p2 = C1.point(a + da);

	sol.push_back(p1);

	if (p1 == p2)
		return 1; //相切
	return sol.push_back(p2), 2;
}

Line getTangent(const Coord &C, const Coord &P) { return Line(P, Normal(C - P)); } //圆心C,圆上一点P处切线

int getTangents(const Coord &p, const Circle &C, vector<Coord> &sol) //点到圆的切点,返回个数
{
	Coord u = p - C.c;
	lf d = abs(u);
	if (d < C.r)
		return 0;	  //点在圆内
	if (!sgn(d - C.r)) //点在圆上
		return sol.push_back(p), 1;
	lf base = arg(u), ang = acos(C.r / d);
	sol.push_back(C.point(base + ang));
	sol.push_back(C.point(base - ang));
	return 2;
}

int getTangents(Circle A, Circle &B, vector<Coord> &a, vector<Coord> &b) //公共切线的切点
{
	int cnt = 0;

	if (A.r < B.r)
		swap(A, B), swap(a, b); //有时需标记交换

	lf d = abs(A.c - B.c),
	   rdiff = A.r - B.r,
	   rsum = A.r + B.r;

	if (sgn(d - rdiff) < 0)
		return 0; //内含

	lf base = arg(B.c - A.c);

	if (!sgn(d) && !sgn(rdiff))
		return -1; //重合,无穷多条切线

	if (!sgn(d - rdiff)) //内切,外公切线
	{
		a.push_back(A.point(base));
		b.push_back(B.point(base));
		return 1;
	}

	//有外公切线的情形
	lf ang = acos(rdiff / d);
	a.push_back(A.point(base + ang));
	b.push_back(B.point(base + ang));
	a.push_back(A.point(base - ang));
	b.push_back(B.point(base - ang));
	cnt += 2;

	if (!sgn(d - rsum))
	{
		a.push_back(A.point(base));
		b.push_back(B.point(base + M_PI));
		++cnt;
	}
	else if (sgn(d - rsum) > 0)
	{
		lf ang_in = acos(rsum / d);
		a.push_back(A.point(base + ang_in));
		b.push_back(B.point(base + ang_in + M_PI));
		a.push_back(A.point(base - ang_in));
		b.push_back(B.point(base - ang_in + M_PI));
		cnt += 2;
	}
	return cnt;
}

lf AreaCircleWithTriangle(const Circle &C, Coord A, Coord B) //C和三角形OAB的相交面积,如果三角形顶点不在O上则把圆和三角形同时平移,直到有一个顶点在O上
{
	int sg = sgn(Cross(A, B));
	if (!sg || A == C.c || B == C.c)
		return 0;

	lf OA = abs(A - C.c), OB = abs(B - C.c), angle = Angle(A, B),
	   d = DistanceToLine(Coord(), A, B);

	if (sgn(OA - C.r) <= 0 && sgn(OB - C.r) <= 0)
		return Cross(A, B) / 2;

	if (sgn(OA - C.r) >= 0 && sgn(OB - C.r) >= 0 && sgn(d - C.r) >= 0)
		return sg * C.r * C.r * angle / 2;
	if (sgn(OA - C.r) >= 0 && sgn(OB - C.r) >= 0 && sgn(d - C.r) < 0)
	{
		Coord prj = getLineProjection(Coord(), A, B);
		if (!onSegment(prj, A, B))
			return sg * C.r * C.r * angle / 2;

		vector<Coord> p;
		Line L = Line(A, B - A);
		getLineCircleIntersection(L, C, p);

		lf s1 = C.r * C.r * angle / 2,
		   s2 = C.r * C.r * Angle(p[0], p[1]) / 2;
		s2 -= fabs(Cross(p[0], p[1]) / 2);
		s1 = s1 - s2;

		return sg * s1;
	}
	if (sgn(OB - C.r) < 0)
		swap(A, B);

	Line L = Line(A, B - A);
	vector<Coord> inter;
	getLineCircleIntersection(L, C, inter);
	Coord inter_point = inter[!onSegment(inter[0], A, B)];

	lf s = fabs(Cross(inter_point, A) / 2);
	s += C.r * C.r * Angle(inter_point, B) / 2;

	return s * sg;
}

lf AreaCircleWithPolygon(const Circle &C, const vector<Coord> &p)
{
	lf ans = 0;
	for (int i = 0; i < p.size(); ++i)
		ans += AreaCircleWithTriangle(C, p[i], p[(i + 1) % p.size()]);
	return fabs(ans);
}

Coord getGravityCenter(const vector<Coord> &p) //多边形重心
{
	Coord a(0, 0);
	lf am = 0, mj;
	for (int i = 0; i < p.size(); ++i)
	{
		mj = Cross(p[i], p[(i + 1) % p.size()]);
		a += mj * (p[i] + p[(i + 1) % p.size()]);
		am += mj;
	}
	return a / am / 3.0;
}

三维

typedef double lf;
const lf EPS = 1e-9, INF = 1e9;
int sgn(lf d) { return (d > EPS) - (d < -EPS); }
struct Coord3
{
	lf X, Y, Z;
	friend bool operator!=(const Coord3 &a, const Coord3 &b) { return sgn(a.X - b.X) || sgn(a.Y - b.Y) || sgn(a.Z - b.Z); }
	friend bool operator==(const Coord3 &a, const Coord3 &b) { return !(a != b); }
	Coord3 &operator+=(const Coord3 &b) { return X += b.X, Y += b.Y, Z += b.Z, *this; }
	friend Coord3 operator+(Coord3 a, const Coord3 &b) { return a += b; }
	Coord3 &operator-=(const Coord3 &b) { return X -= b.X, Y -= b.Y, Z -= b.Z, *this; }
	friend Coord3 operator-(Coord3 a, const Coord3 &b) { return a -= b; }
	Coord3 &operator*=(lf d) { return X *= d, Y *= d, Z *= d, *this; }
	friend Coord3 operator*(Coord3 a, lf d) { return a *= d; }
	friend Coord3 operator*(lf d, Coord3 a) { return a *= d; }
	Coord3 &operator/=(lf d) { return X /= d, Y /= d, Z /= d, *this; }
	friend Coord3 operator/(Coord3 a, lf d) { return a /= d; }
	friend lf Dot(const Coord3 &A, const Coord3 &B) { return A.X * B.X + A.Y * B.Y + A.Z * B.Z; }
	friend Coord3 Cross(const Coord3 &A, const Coord3 &B) { return {A.Y * B.Z - A.Z * B.Y, A.Z * B.X - A.X * B.Z, A.X * B.Y - A.Y * B.X}; }
	friend lf norm(const Coord3 &A) { return Dot(A, A); }
	friend lf abs(const Coord3 &A) { return sqrt(norm(A)); }
	friend lf Angle(const Coord3 &A, const Coord3 &B) { return acos(Dot(A, B) / abs(A) / abs(B)); }
	friend lf Area2(Coord3 A, Coord3 B, Coord3 C) { return abs(Cross(B - A, C - A)); }
	friend lf Volume6(Coord3 A, Coord3 B, Coord3 C, Coord3 D) { return Dot(D - A, Cross(B - A, C - A)); }	 //四面体体积
	friend Coord3 Centroid(Coord3 A, Coord3 B, Coord3 C, Coord3 D) { return (A + B + C + D) / 4.0; }		  //四面体的重心
	friend lf DistanceToPlane(Coord3 p, Coord3 p0, const Coord3 &n) { return Dot(p - p0, n) / abs(n); }		  //点p到平面p0-n的有向距离
	friend Coord3 getPlaneProjection(Coord3 p, Coord3 p0, const Coord3 &n) { return p - n * Dot(p - p0, n); } //点p在平面p0-n上的投影。n必须为单位向量
	friend Coord3 LinePlaneIntersection(Coord3 p1, Coord3 p2, Coord3 p0, Coord3 n)							  //直线p1-p2 与平面p0-n的交点,假设交点唯一存在
	{
		Coord3 v = p2 - p1;
		lf t = Dot(n, p0 - p1) / Dot(n, p2 - p1); //分母为0,直线与平面平行或在平面上
		return p1 + v * t;						  //如果是线段 判断t是否在0~1之间
	}
	friend lf DistanceToLine(Coord3 P, Coord3 A, Coord3 B) //点P到直线AB的距离
	{
		Coord3 v1 = B - A, v2 = P - A;
		return abs(Cross(v1, v2)) / abs(v1);
	}
	friend lf DistanceToSeg(Coord3 P, Coord3 A, Coord3 B) //点到线段的距离
	{
		if (A == B)
			return abs(P - A);
		Coord3 v1 = B - A, v2 = P - A, v3 = P - B;
		if (sgn(Dot(v1, v2)) < 0)
			return abs(v2);
		if (sgn(Dot(v1, v3)) > 0)
			return abs(v3);
		return fabs(DistanceToLine(P, A, B));
	}
	friend bool LineDistance3D(Coord3 p1, Coord3 u, Coord3 p2, Coord3 v, lf &s) //求异面直线 p1+s*u与p2+t*v的公垂线对应的s,如果平行|重合,返回0
	{
		lf b = Dot(u, u) * Dot(v, v) - Dot(u, v) * Dot(u, v);
		if (!sgn(b))
			return 0;
		lf a = Dot(u, v) * Dot(v, p1 - p2) - Dot(v, v) * Dot(u, p1 - p2);
		return s = a / b, 1;
	}
	friend bool SameSide(Coord3 p1, Coord3 p2, Coord3 a, Coord3 b) { return sgn(Dot(Cross(b - a, p1 - a), Cross(b - a, p2 - a))) >= 0; } //p1和p2是否在线段a-b的同侧
	friend bool PointInTri(Coord3 PP, Coord3 P[3])																						 //点P在三角形P0,P1,p中
	{
		return SameSide(PP, P[0], P[1], P[2]) &&
			   SameSide(PP, P[1], P[0], P[2]) &&
			   SameSide(PP, P[2], P[0], P[1]);
	}
	friend bool TriSegIntersection(Coord3 P[3], Coord3 A, Coord3 B, Coord3 &PP) //三角形P0P1p是否和线段AB相交,如有则为PP
	{
		Coord3 n = Cross(P[1] - P[0], P[2] - P[0]);
		if (sgn(Dot(n, B - A)) == 0)
			return false;						 //线段A-B和平面P0P1p平行或共面
		lf t = Dot(n, P[0] - A) / Dot(n, B - A); //平面A和直线P1-p有惟一交点
		if (sgn(t) < 0 || sgn(t - 1) > 0)
			return false; //不在线段AB上
		return PointInTri(PP = A + (B - A) * t, P);
	}
	friend bool TriTriIntersection(Coord3 T1[3], Coord3 T2[3]) //空间两三角形是否相交
	{
		Coord3 P;
		for (int i = 0; i < 3; ++i)
			if (TriSegIntersection(T1, T2[i], T2[(i + 1) % 3], P) ||
				TriSegIntersection(T2, T1[i], T1[(i + 1) % 3], P))
				return 1;
		return 0;
	}
	friend lf SegSegDistance(Coord3 a1, Coord3 b1, Coord3 a2, Coord3 b2, Coord3 &ans1, Coord3 &ans2) //空间两直线上最近点对 返回最近距离 点对保存在ans1 ans2中
	{
		Coord3 v1 = (a1 - b1), v2 = (a2 - b2);
		Coord3 N = Cross(v1, v2);
		Coord3 ab = (a1 - a2);
		lf ans = Dot(N, ab) / abs(N);
		Coord3 d1 = b1 - a1, d2 = b2 - a2, cd = Cross(d1, d2);
		lf nd = norm(cd), t1 = Dot(Cross(a2 - a1, d2), cd) / nd, t2 = Dot(Cross(a2 - a1, d1), cd) / nd;
		return ans1 = a1 + (b1 - a1) * t1, ans2 = a2 + (b2 - a2) * t2, fabs(ans);
	}
	friend bool InsideWithMinDistance(Coord3 PP, Coord3 *P, lf dist) //判断PP是否在三角形P中,并且到三条边的距离都至少为dist。保证P,A,B,C共面
	{
		return PointInTri(PP, P) &&
				   DistanceToLine(PP, P[0], P[1]) >= dist ||
			   DistanceToLine(PP, P[1], P[2]) >= dist ||
			   DistanceToLine(PP, P[2], P[0]) >= dist;
	}
	friend lf minBall(const vector<Coord3> &data, const lf eps = EPS * 1e-3) //模拟退火求最小球覆盖,EPS玄学调整;返回半径
	{
		lf step = 1, ans = INF;
		for (Coord3 z{0, 0, 0}; step > eps; step *= 0.99)
		{
			int s = 0;
			for (int i = 0; i < data.size(); ++i)
				if (abs(z - data[s]) < abs(z - data[i]))
					s = i;
			ans = min(ans, abs(z - data[s]));
			z -= (z - data[s]) * step;
		}
		return ans;
	}
};
struct Sphere
{
	Coord3 o;
	lf r;
	Coord3 point(lf lat, lf lng) const //经纬度确定球面上一点
	{
		lat *= M_PI / 180;
		lng *= M_PI / 180;
		return o + r * Coord3{cos(lat) * cos(lng), cos(lat) * sin(lng), sin(lat)};
	}
};
struct ConvexPolyhedron //空间多边形和凸包问题
{
	struct Face
	{
		int v[3];
		Face(int a, int b, int c) { v[0] = a, v[1] = b, v[2] = c; }
		Coord3 Normal(const vector<Coord3> &P) const { return Cross(P[v[1]] - P[v[0]], P[v[2]] - P[v[0]]); }
		bool CanSee(const vector<Coord3> &P, int i) const { return Dot(P[i] - P[v[0]], Normal(P)) > 0; } //f是否能看见P[i]
	};
	vector<Face> faces;
	vector<Coord3> p;
	ConvexPolyhedron(vector<Coord3> P) : p(P)
	{
		for (int i = 0; i < p.size(); ++i)
			P[i] += Coord3{randEPS(), randEPS(), randEPS()};
		vector<vector<int>> vis(P.size(), vector<int>(P.size()));
		faces.push_back(Face(0, 1, 2)); //由于已经进行扰动,前三个点不共线
		faces.push_back(Face(2, 1, 0));
		for (int i = 3; i < P.size(); ++i)
		{
			vector<Face> next;
			for (int j = 0; j < faces.size(); ++j) //计算每条边的「左面」的可见性
			{
				Face &f = faces[j];
				int res = f.CanSee(P, i);
				if (!res)
					next.push_back(f);
				for (int k = 0; k < 3; ++k)
					vis[f.v[k]][f.v[(k + 1) % 3]] = res;
			}
			for (int j = 0; j < faces.size(); ++j)
				for (int k = 0; k < 3; ++k)
				{
					int a = faces[j].v[k], b = faces[j].v[(k + 1) % 3];
					if (vis[a][b] != vis[b][a] && vis[a][b]) //(a,b)是分界线,左边对P[i]可见
						next.push_back(Face(a, b, i));
				}
			swap(faces, next);
		}
	}
	lf randEPS() { return (rand() / lf(RAND_MAX) - 0.5) * EPS; }
	Coord3 centroid() //三维凸包重心
	{
		Coord3 C = p[0], tot{0, 0, 0};
		lf totv = 0;
		for (int i = 0; i < faces.size(); ++i)
		{
			Coord3 p1 = p[faces[i].v[0]], p2 = p[faces[i].v[1]], p3 = p[faces[i].v[2]];
			lf v = -Volume6(p1, p2, p3, C);
			totv += v;
			tot += v * Centroid(p1, p2, p3, C);
		}
		return tot / totv;
	}
	lf dist(Coord3 C) //凸包内一点到表面最近距离
	{
		lf ans = INF;
		for (int i = 0; i < faces.size(); ++i)
		{
			Coord3 p1 = p[faces[i].v[0]], p2 = p[faces[i].v[1]], p3 = p[faces[i].v[2]];
			ans = min(ans, fabs(-Volume6(p1, p2, p3, C) / Area2(p1, p2, p3)));
		}
		return ans;
	}
};

高精度

无符号整数

使用范例

vector 自带大小比较为字典序比较, !=== 运算,其余需要时一定记得重载!

减法,当被减数小于减数时减为 0。

struct Wint : vector<int> //继承vector
{
	static const int width = 9, base = 1e9;
	Wint(unsigned long long n = 0) //普通初始化,当整型数和Wint同时运算时会提升至Wint
	{
		for (; n; n /= base)
			push_back(n % base);
	}
	explicit Wint(const string &s) //字符串初始化函数,未判断字符串合法情况
	{
		for (int len = int(s.size() - 1) / width + 1, b, e, i = 0; i < len; ++i)
			for (e = s.size() - i * width, b = max(0, e - width), push_back(0); b != e; ++b)
				back() = back() * 10 + s[b] - '0';
		trim(0);
	}
	Wint &trim(bool up = 1) //去前导0,是否需要进位,很常用的小函数,为方便返回自身
	{
		for (int i = 1; up && i < size(); ++i)
		{
			if ((*this)[i - 1] < 0)
				--(*this)[i], (*this)[i - 1] += base;
			if ((*this)[i - 1] >= base)
				(*this)[i] += (*this)[i - 1] / base, (*this)[i - 1] %= base;
		}
		while (!empty() && back() <= 0)
			pop_back();
		for (; up && !empty() && back() >= base; (*this)[size() - 2] %= base)
			push_back(back() / base);
		return *this;
	}
	friend istream &operator>>(istream &is, Wint &n)
	{
		string s; //懒
		return is >> s, n = Wint(s), is;
	}
	friend ostream &operator<<(ostream &os, const Wint &n)
	{
		if (n.empty())
			return os.put('0');
		os << n.back();
		char ch = os.fill('0');
		for (int i = n.size() - 2; ~i; --i)
			os.width(n.width), os << n[i];
		return os.fill(ch), os;
	}
	friend bool operator<(const Wint &a, const Wint &b)
	{
		if (a.size() != b.size())
			return a.size() < b.size();
		for (int i = a.size() - 1; ~i; --i)
			if (a[i] != b[i])
				return a[i] < b[i];
		return 0;
	}
	friend bool operator>(const Wint &a, const Wint &b) { return b < a; }
	friend bool operator<=(const Wint &a, const Wint &b) { return !(a > b); }
	friend bool operator>=(const Wint &a, const Wint &b) { return !(a < b); }
	Wint &operator+=(const Wint &b)
	{
		if (size() < b.size())
			resize(b.size()); //保证有足够的位数
		for (int i = 0; i < b.size(); ++i)
			(*this)[i] += b[i];
		return trim(); //单独进位防自运算
	}
	friend Wint operator+(Wint a, const Wint &b) { return a += b; }
	Wint &operator++() { return *this += 1; } //前置版本,懒
	Wint operator++(int)					  //后置版本
	{
		Wint b(*this);
		return ++*this, b;
	}
	Wint &operator-=(const Wint &b) //a<b会使a变为0
	{
		if (size() < b.size())
			resize(b.size()); //保证有足够的位数
		for (int i = 0; i < b.size(); ++i)
			(*this)[i] -= b[i];
		return trim(); //单独进位防自运算
	}
	friend Wint operator-(Wint a, const Wint &b) { return a -= b; }
	Wint &operator--() { return *this -= 1; } //前置版本,懒
	Wint operator--(int)					  //后置版本
	{
		Wint b(*this);
		return --*this, b;
	}
	Wint &operator*=(const Wint &b) //高精度乘法,常规写法
	{
		Wint c;
		c.assign(size() + b.size(), 0);
		for (int j = 0, k, l; j < b.size(); ++j)
			if (b[j]) //稀疏优化,特殊情况很有效
				for (int i = 0; i < size(); ++i)
				{
					unsigned long long n = (*this)[i];
					for (n *= b[j], k = i + j; n; n /= base)
						c[k++] += n % base;
					for (l = i + j; c[l] >= base || l + 1 < k; c[l++] %= base)
						c[l + 1] += c[l] / base;
				}
		return swap(c), trim(0);
	}
	/*
	Wint &operator*=(const Wint &b) //一种效率略高但对位宽有限制的写法
	{
		vector<unsigned long long> n(size() + b.size(), 0); //防爆int
		//乘法算完后统一进位效率高,防止乘法溢出(unsigned long long范围0~1.8e19)
		//位宽为9时size()不能超过18(十进制162位),位宽为8时size()不能超过1800(十进制14400位)等等。
		for (int j = 0; j != b.size(); ++j)
			if (b[j]) //稀疏优化,特殊情况很有效
				for (int i = 0; i != size(); ++i)
					n[i + j] += (unsigned long long)(*this)[i] * b[j];
		for (int i = 1; i < n.size(); ++i) //这里用<防止位数0,单独进位防自运算
			n[i] += n[i - 1] / base, n[i - 1] %= base;
		return assign(n.begin(), n.end()), trim(0);
	}
	Wint &operator*=(const Wint &b) //fft优化乘法,注意double仅15位有效数字,调小Wint::width不超过2,计算自2*log2(base)+2*log2(len)<53
	{
		vector<ll> ax(begin(), end()), bx(b.begin(), b.end());
		ax = FFT(size() + b.size()).ask(ax, bx);
		for (int i = 1; i < ax.size(); ++i)
			ax[i] += ax[i - 1] / base, ax[i - 1] %= base;
		return assign(ax.begin(), ax.end()), trim(0);
	}
	Wint &operator*=(const Wint &b) //ntt优化,Wint::width不超过2
	{
		vector<ll> ax(begin(), end()), bx(b.begin(), b.end());
		ax = FNTT(size() + b.size(), (7 << 26) + 1, 3).ask(ax, bx);
		for (int i = 1; i < ax.size(); ++i)
			ax[i] += ax[i - 1] / base, ax[i - 1] %= base;
		return assign(ax.begin(), ax.end()), trim(0);
	}
	*/
	friend Wint operator*(Wint a, const Wint &b) { return a *= b; }
	Wint &operator/=(Wint b)
	{
		Wint r, c, d = b.base / (b.back() + 1);
		*this *= d, b *= d, c.assign(size(), 0);
		for (int i = size() - 1; ~i; --i)
		{
			r.insert(r.begin(), (*this)[i]);
			unsigned long long s = 0;
			for (int j = b.size(); j + 1 >= b.size(); --j) //b.size()==0肯定第一行就出问题的
				s = s * b.base + (j < r.size() ? r[j] : 0);
			for (d = c[i] = s / b.back(), d *= b; r < d; r += b)
				--c[i];
			r -= d;
		}
		return swap(c), trim(0); //r为加倍后的余数,可通过高精度除低精度得到真正余数,此处略
	}
	friend Wint operator/(Wint a, const Wint &b) { return a /= b; }
	Wint &operator%=(const Wint &b) { return *this -= *this / b * b; }
	friend Wint operator%(Wint a, const Wint &b) { return a %= b; }
	//开平方,改自ZJU模板
	bool cmp(long long c, int d, const Wint &b) const
	{
		if ((int)b.size() - (int)size() < d + 1 && c)
			return 1;
		long long t = 0;
		for (int i = b.size() - 1, lo = -(base << 1); lo <= t && t <= 0 && ~i; --i)
			if (t = t * base - b[i], 0 <= i - d - 1 && i - d - 1 < size())
				t += (*this)[i - d - 1] * c;
		return t > 0;
	}
	Wint &sub(const Wint &b, long long k, int d)
	{
		int l = b.size() + d;
		for (int i = d + 1; i <= l; ++i)
		{
			long long tmp = (*this)[i] - k * b[i - d - 1];
			if (tmp < 0)
			{
				(*this)[i + 1] += (tmp - base + 1) / base;
				(*this)[i] = tmp - (tmp - base + 1) / base * base;
			}
			else
				(*this)[i] = tmp;
		}
		for (int i = l + 1; i < size() && (*this)[i] < 0; ++i)
		{
			(*this)[i + 1] += ((*this)[i] - base + 1) / base;
			(*this)[i] -= ((*this)[i] - base + 1) / base * base;
		}
		return trim(0);
	}
	friend Wint sqrt(Wint a)
	{
		Wint n;
		n.assign(a.size() + 1 >> 1, 0);
		for (int i = n.size() - 1, l, r; ~i; --i)
		{
			for (l = 0, r = a.base, n[i] = l + r >> 1; r - l > 1; n[i] = l + r >> 1)
			{
				if (n.cmp(n[i], i - 1, a))
					r = n[i];
				else
					l = n[i];
			}
			a.sub(n, n[i], i - 1), n[i] += l + r >> 1;
		}
		for (int i = 0; i < n.size(); ++i)
			n[i] >>= 1;
		return n.trim(0);
	}
	/*
	friend Wint sqrt(const Wint &a) //常规牛顿迭代实现的开平方算法,慢但是好敲
	{
		Wint b = a, c = (b + 1) / 2;
		while (b != c)
			swap(b, c), c = (b + a / b) / 2;
		return c;
	}
	friend Wint sqrt(const Wint &a)
	{
		Wint ret, t;
		ret.assign((a.size() + 1) >> 1, 0);
		for (int i = ret.size() - 1, l, r; ~i; --i)
		{
			for (l = 0, r = a.base; r - l > 1;)
			{
				ret[i] = l + (r - l) / 2;
				t = ret * ret;
				if (a < t)
					r = ret[i];
				else
					l = ret[i];
			}
			if (!l && i == ret.size() - 1)
				ret.pop_back();
			else
				ret[i] = l;
		}
		return ret;
	}
	*/
};

分数

使用示范

struct Fraction
{
	ll num, den;
	Fraction(ll n = 0, ll d = 1) : num(n), den(d)
	{
		d = __gcd(num, den), num /= d, den /= d;
		if (den < 0)
			num = -num, den = -den;
	}
	friend Fraction operator+(const Fraction &A, const Fraction &B)
	{
		ll d = __gcd(A.den, B.den);
		return Fraction(B.den / d * A.num + A.den / d * B.num, A.den / d * B.den);
	}
	Fraction &operator+=(const Fraction &c) { return *this = *this + c; }
	Fraction operator-() const
	{
		Fraction r(*this);
		return r.num = -r.num, r;
	}
	friend Fraction operator-(const Fraction &a, const Fraction &c) { return -c + a; }
	Fraction &operator-=(const Fraction &c) { return *this = *this - c; }
	friend Fraction operator*(const Fraction &A, const Fraction &B) { return Fraction(A.num * B.num, A.den * B.den); }
	Fraction &operator*=(const Fraction &c) { return *this = *this * c; }
	friend Fraction operator/(const Fraction &A, const Fraction &B) { return Fraction(A.num * B.den, A.den * B.num); }
	Fraction &operator/=(const Fraction &c) { return *this = *this / c; }
	friend Fraction operator%(const Fraction &a, const Fraction &c) { return Fraction(a.num * c.den % (c.num * a.den), a.den * c.den); }
	Fraction &operator%=(const Fraction &c) { return *this = *this % c; }
	friend bool operator==(const Fraction &a, const Fraction &b) { return a.num * b.den == a.den * b.num; }
	friend bool operator!=(const Fraction &a, const Fraction &b) { return !(a == b); }
	friend bool operator<(const Fraction &a, const Fraction &b) { return a.num * b.den < a.den * b.num; }
	friend bool operator>(const Fraction &a, const Fraction &b) { return b < a; }
	friend bool operator<=(const Fraction &a, const Fraction &b) { return !(a > b); }
	friend bool operator>=(const Fraction &a, const Fraction &b) { return !(a < b); }
	friend Fraction abs(Fraction f)
	{
		if (f.num < 0)
			f.num = -f.num;
		return f;
	}
	friend ostream &operator<<(ostream &os, const Fraction &f) { return !f.num ? os << 0 : f.den == 1 ? os << f.num : os << f.num << '/' << f.den; }
};

bitset 高精度

代替整型进行位运算,更方便并且可以处理超过最大整形范围大小的位集合。 你可以把 bitset 看作可以位运算的 bool 数组,换言之,bitset 的大小是固定的。因此,用 bitset 做状态压缩是很方便的,也可以方便的实现集合的交并补操作。 bitset 仅重载了相等不等和位运算符,原生不支持四则运算和大小比较,所以很少代替高精度数。

typedef bitset<127> Bint;
/*
Bint b(unsigned long long u=0);//用u的低N位初始化b,N是一个可转成ULL类型的常量表达式,高位补0
Bint bs(string s,int pos,int m=string::npos,char zero='0',char one='1');//用s从pos位开始的m位初始化b,s中只含zero和one

b.size();//b的大小,即N
b.count();//b中1的个数

b[pos];//b中pos位的引用

b.set();//b全赋1
b.reset();//b全赋0
b.flip();//b全反转

b.to_ull();//b转成unsigned long long,b.size()>64时出错
b.to_string(char zero='0',char one='1');//按参数输出字符串

os<<b;//按'0'和'1'打印b
is>>b;//按'0'和'1'读入b,当下一个字符不是'0'或'1'或读到b.size()个数后停止

==、!=、<<、>>、&、|、^//保持它们在整型运算中的含义
*/
//大小比较,其他运算符类似
bool operator<(const Bint &a, const Bint &b)
{
	for (int i = a.size() - 1; ~i; --i)
		if (a[i] != b[i])
			return a[i] < b[i];
	return 0;
}
bool operator!=(const Bint &a, const Bint &b)
{
	for (int i = a.size() - 1; ~i; --i)
		if (a[i] != b[i])
			return 1;
	return 0;
}
//加法
Bint operator+(const Bint &a, const Bint &b) { return b.any() ? (a ^ b) + ((a & b) << 1) : a; }
Bint &operator+=(Bint &a, const Bint &b) { return a = a + b; }
//减法
Bint operator-(const Bint &a) { return Bint(1) + ~a; }
Bint &operator-=(Bint &a, const Bint &b) { return a += -b; }
Bint operator-(Bint a, const Bint &b) { return a -= b; }
//乘法
Bint operator*(Bint a, Bint b)
{
	Bint r(0);
	for (; b.any(); b >>= 1, a += a)
		if (b[0])
			r += a;
	return r;
}
Bint &operator*=(Bint &a, const Bint &b) { return a = a * b; }
//整除,取模
Bint operator%=(Bint &dividend, const Bint &divisor)
{
	if (dividend < divisor || divisor.none())
		return dividend;
	Bint c, res = 0;
	for (int k; divisor < dividend;)
	{
		for (k = 0, c = divisor; !(dividend < c); c <<= 1, ++k)
			if (dividend < divisor + c)
			{
				res += 1 << k;
				break;
			}
		if (dividend < divisor + c)
			break;
		res += 1 << (k - 1);
		dividend -= c >> 1;
	}
	return dividend; //res是商
}
//输入输出,bitset已经原生重载了输入输出运算符,避免歧义。
istream &getb(istream &is, Bint &val)
{
	int sign = 1, ch = is.get();
	for (; !isdigit(ch) && ch != EOF; ch = is.get())
		if (ch == '-')
			sign = -sign;
	for (val = 0; isdigit(ch); ch = is.get())
		val = (val << 3) + (val << 1) + (ch ^ '0');
	if (sign == -1)
		val = -val;
	return is.putback(ch);
}
ostream &putb(ostream &os, const Bint &val)
{
	if (Bint(9) < val)
		putb(os, val / 10);
	return os.put(val.to_ulong() % 10 + '0');
}

C++语言相关

GCC 内置位运算

int __builtin_clz(unsigned int x);      //求前缀0的个数
int __builtin_ctz(unsigned int x);      //求后缀0的个数
int __builtin_ffs(unsigned int x);      //x的二进制末尾最后一个1的位置,从1开始
int __builtin_popcount(unsigned int x); //x二进制中1的个数,相当于bitset::count()
int __builtin_parity(unsigned int x);   //判断x的二进制中1的个数的奇(1)偶(0)性,这些函数都有相应的usigned long和usigned long long版本,只需要在函数名后面加上l或ll就可以了,比如__builtin_clzll

开栈

For C++

#pragma comment(linker, "/STACK:102400000,102400000") //For C++

For G++

int main() //For G++
{
	int size = 256 << 20; //256MB
	char *p = (char *)malloc(size) + size;
	__asm__ __volatile__("movq %0, %%rsp\n" ::"r"(p)); //64bit,一定要最后写一句`exit(0);`退出程序,否则会得到非零退出的错误,可能RE。
}

读入优化

C 文件指针版

ll getll(FILE *in = stdin)
{
	ll val = 0, sgn = 1, ch = getc(in);
	for (; !isdigit(ch) && ch != EOF; ch = getc(in))
		if (ch == '-')
			sgn = -sgn;
	for (; isdigit(ch); ch = getc(in))
		val = val * 10 + ch - '0';
	return ungetc(ch, in), sgn * val;
}
lf getlf(FILE *in = stdin)
{
	lf val = getll(in), p = val < 0 ? -1 : 1;
	ll ch = getc(in);
	if (ch == '.')
		for (ch = getc(in); isdigit(ch); ch = getc(in))
			val += (p /= 10) * (ch - '0');
	return ungetc(ch, in), val;
}

仿 C++IO 流沙雕版

使用前

使用后

#define cin kin
struct Istream
{
	char b[20 << 20], *i, *e; //20MB
	Istream(FILE *in) : i(b), e(b + fread(b, sizeof(*b), sizeof(b), in)) {}
	bool eof() const
	{
		return i == e;
	}
	Istream &operator>>(long long &val)
	{
		return val = strtoll(i, &i, 10 /*进制,取值2~36*/), *this;
	}
	Istream &operator>>(ll &val) //极限快
	{
		while (*i < '0')
			++i; //无符号
		for (val = 0; *i >= '0'; ++i)
			val = (val << 3) + (val << 1) + *i - '0';
		return *this;
	}
	Istream &operator>>(double &val)
	{
		return val = strtod(i, &i), *this;
	}
	Istream &operator>>(string &s)
	{
		while (!eof() && isspace(*i))
			++i;
		for (s.clear(); !eof() && !isspace(*i); ++i)
			s += *i;
		return *this;
	}
} kin(stdin);