ACM Templates

.ACM-Templates/.Not classified/full.cpp

@@ -0,0 +1,45 @@
+#if __cplusplus >= 201103L
+#include <bits/stdc++.h>
+#endif
+#include <algorithm>
+#include <bitset>
+#include <cassert>
+#include <cctype>
+#include <climits>
+#include <cmath>
+#include <complex>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <ctime>
+#include <functional>
+#include <iomanip>
+#include <iostream>
+#include <list>
+#include <map>
+#include <numeric>
+#include <queue>
+#include <set>
+#include <stack>
+#include <string>
+#include <vector>
+using namespace std;
+
+typedef long long ll;
+const int N = 100005;
+const int M = 1000000007;
+
+int n, m;
+int a[N];
+
+int main() {
+  int C = 0, T;
+  scanf("%d", &T);
+  while (++C <= T) {
+    scanf("%d", &n);
+    for (int i = 0; i < n; i++) {
+      scanf("%d", &a[i]);
+    }
+
+  }
+}

.ACM-Templates/.Not classified/minimized.cpp

@@ -0,0 +1,18 @@
+#include <bits/stdc++.h>
+using namespace std;
+
+typedef long long ll;
+const int N = 100005;
+const int M = 1000000007;
+
+int n, m;
+int a[N];
+
+int main() {
+  while (~scanf("%d", &n)) {
+    for (int i = 0; i < n; i++) {
+      scanf("%d", &a[i]);
+    }
+
+  }
+}

.ACM-Templates/Data-Structures/link-cut-tree/main.cpp

@@ -0,0 +1,80 @@
+const int MAXN = 100010;
+struct Node {
+	Node *ch[2], *p; int size, value;
+	bool rev;
+	Node(int t = 0);
+	inline bool dir(void) {return p->ch[1] == this;}
+	inline void SetC(Node *x, bool d) {
+		ch[d] = x; x->p = this;
+	}
+	inline void Rev(void) {
+		swap(ch[0], ch[1]); rev ^= 1;
+	}
+	inline void Push(void) { 
+		if (rev) {
+			ch[0]->Rev();
+			ch[1]->Rev();
+			rev = 0;
+		}
+	}
+	inline void Update(void) { 
+		size = ch[0]->size + ch[1]->size + 1;
+	}
+}Tnull, *null = &Tnull, *fim[MAXN];
+// 要记得额外更新null的信息
+Node::Node(int _value){ch[0] = ch[1] = p = null; rev = 0;}
+inline bool isRoot(Node *x) {return x->p == null || (x != x->p->ch[0] && x != x->p->ch[1]);}
+inline void rotate(Node *x) {
+	Node *p = x->p; bool d = x->dir();
+	p->Push(); x->Push();
+	if (!isRoot(p)) p->p->SetC(x, p->dir()); else x->p = p->p;
+	p->SetC(x->ch[!d], d);
+	x->SetC(p, !d);
+	p->Update();
+}
+inline void splay(Node *x) {
+	x->Push();
+	while (!isRoot(x)) {
+		if (isRoot(x->p)) rotate(x);
+		else {
+			if (x->dir() == x->p->dir()) {rotate(x->p); rotate(x);} 
+			else {rotate(x); rotate(x);}
+		}
+	}
+	x->Update();
+}
+inline Node* Access(Node *x) {
+	Node *t = x, *q = null;
+	for (; x != null; x = x->p) {
+		splay(x); x->ch[1] = q; q = x;
+	}
+	splay(t); //info will be updated in the splay;
+	return q;
+}
+inline void Evert(Node *x) {
+	Access(x); x->Rev();
+}
+inline void link(Node *x, Node *y) {
+	Evert(x); x->p = y;
+}
+inline Node* getRoot(Node *x) {
+	Node *tmp = x;
+	Access(x);
+	while (tmp->Push(), tmp->ch[0] != null) tmp = tmp->ch[0];
+	splay(tmp);
+	return tmp;
+}
+// 一定要确定x和y之间有边
+inline void cut(Node *x, Node *y) { 	
+	Access(x); splay(y);
+	if (y->p != x) swap(x, y);
+	Access(x); splay(y);
+	y->p = null;
+}
+inline Node* getPath(Node *x, Node *y) {
+	Evert(x); Access(y);
+	return y;
+}
+inline void clear(void) {
+	null->rev = 0; null->sie = 0; null->value = 0;
+}

.ACM-Templates/Data-Structures/link-cut-tree/remark.tex

@@ -0,0 +1 @@
+动态树模板，如果维护的权值在边上，需要把边拆成点。

.ACM-Templates/Data-Structures/splay/main.cpp

@@ -0,0 +1,158 @@
+// 最大可能的节点数
+const int MAXN = 100010;
+// 每个打标记的操作就是更新这个节点的信息，然后对子节点打标记
+struct Node {
+	Node *ch[2], *p; int size, value;
+	bool rev;
+	inline bool dir(void) {return p->ch[1] == this;}
+	inline void SetC(Node *x, bool d) {
+		ch[d] = x; x->p = this;
+	}
+	inline void Rev(void) {
+		swap(ch[0], ch[1]); rev ^= 1;
+	}
+	// null永远不会push
+	inline void Push(void) { 
+		if (rev) {
+			ch[0]->Rev();
+			ch[1]->Rev();
+			rev = 0;
+		}
+	}
+	// null永远不会update
+	inline void Update(void) { 
+		size = ch[0]->size + ch[1]->size + 1;
+	}
+	inline void initInfo(void) {
+		rev = 0;
+	}
+}Tnull, *null = &Tnull, *data, POOL[MAXN];
+class Splay {public:
+	Node *root;
+	inline void rotate(Node *x) {
+		Node *p = x->p; bool d = x->dir();
+		p->Push(); x->Push();
+		p->p->SetC(x, p->dir());
+		p->SetC(x->ch[!d], d);
+		x->SetC(p, !d);
+		p->Update();
+	}
+	inline void splay(Node *x, Node *G) {
+		if (G == null) root = x;
+		while (x->p != G) {
+			if (x->p->p == G) rotate(x); else {
+				if (x->dir() == x->p->dir()) {rotate(x->p); rotate(x);}
+				else {rotate(x); rotate(x);}
+			}
+		}
+		x->Push();
+		x->Update();
+	}
+	inline Node* Renew(int value) {
+		Node *ret = data++;
+		ret->ch[0] = ret->ch[1] =ret->p = null;
+		ret->size = 1; ret->value = value;
+		ret->initInfo();
+		return ret;
+	}
+	inline Node* getMin(Node *x) {Node *tmp = x; while (tmp->ch[0] != null) tmp = tmp->ch[0]; return tmp;}
+	inline Node* getMax(Node *x) {Node *tmp = x; while (tmp->ch[1] != null) tmp = tmp->ch[1]; return tmp;}
+	// 查询第k大
+	inline Node* getKth(int k) {
+		Node *tmp = root; 
+		assert(k > 0 && k <= root->size);
+		while (true) {
+			tmp->Push();
+			if (tmp->ch[0]->size + 1 == k) return tmp;
+			if (tmp->ch[0]->size >= k) tmp = tmp->ch[0]; 
+			else k -= tmp->ch[0]->size + 1, tmp = tmp->ch[1];
+		}
+	}
+	// 以下为splay当作平衡树使用
+	// 查找树中value = v的元素, 返回之后splay
+	inline Node* find(int v) {
+		Node *tmp = root;
+		while (tmp != null) {
+			tmp->Push();
+			if (tmp->value == v) return tmp;
+			if (v < tmp->value) tmp = tmp->ch[0]; else tmp = tmp->ch[1];
+		}
+		return null;
+	}
+	// 统计有多少元素小于等于v, 当flag = 1时，统计多少元素严格小于v, 一定要记得splay最后的那个tmp
+	inline int Count(int v, bool flag = 0) {
+		Node *tmp = root, *last = null; int ret = 0;
+		while (tmp != null) {
+			tmp->Push();
+			last = tmp;
+			if ((!flag && tmp->value > v) || (flag && tmp->value >= v)) {
+				tmp = tmp->ch[0];
+			} else ret += tmp->ch[0]->size + 1, tmp = tmp->ch[1];
+		}
+		if (last != null) splay(last, null);
+		return ret;
+	}
+	// 删除x这个结点
+	inline void erase(Node* x) {
+		splay(x, null);
+		if (x->ch[0] == null || x->ch[1] == null) {
+			int d = x->ch[1] != null;
+			root = x->ch[d]; root->p = null;
+			return;
+		}
+		Node *L = getMax(x->ch[0]), *R = getMax(x->ch[1]);
+		splay(L, x); splay(R, x);
+		L->SetC(R, 1); L->p = null; root = L;
+		L->Update();
+	}
+	// 插入一个值为value的节点，初始要以Insert(root, null, value)来调用, 返回之后splay
+	inline Node* Insert(Node *&now, Node* father, int value) {
+		if (now == null) {
+			now = Renew(value); now->p = father;
+			return now;
+		}
+		Node *ret;
+		now->Push();
+		if (value <= now->value) ret = Insert(now->ch[0], now, value);
+		else ret = Insert(now->ch[1], now, value);
+		now->Update();
+		return ret;
+	}
+	// 以下为splay维护序列, 初始要在原序列中放入一个-inf和inf来防止边界条件
+	// 得到原数列中[l,r]区间对应的结点，如果l == r + 1则表示是一个空区间
+	inline Node* getInterval(int l, int r) {
+		assert(l <= r + 1);
+		Node *L = getKth(l), *R = getKth(r + 2);
+		splay(L, null); splay(R, L);
+		return R->ch[0];
+	}
+	// 删除一段区间[l,r]
+	inline void eraseInterval(int l, int r) {
+		getInterval(l, r);
+		root->ch[1]->ch[0] = null;
+		root->ch[1]->Update();
+		root->Update();
+	}
+	// 在位置l的后面插入一段区间x (0 <= l <= n)
+	inline void insertInterval(int l, Node *x) {
+		Node *L = getKth(l + 1), *R = getKth(l + 2);
+		splay(L, null); splay(R, L);
+		R->SetC(x, 0);
+		R->Update(); L->Update();
+	}
+	// 把数列a的[l,r]构建为一个splay
+	inline Node* Build(int l, int r, int a[]) {
+		if (l > r) return null;
+		int mid = (l + r) >> 1;
+		Node *ret = Renew(a[mid]);
+		if (l == r) return ret;
+		ret->SetC(Build(l, mid - 1, a), 0);
+		ret->SetC(Build(mid + 1, r, a), 1);
+		ret->Update();
+		return ret;
+	}
+}T;
+void clear(void) {
+	data = POOL;
+	T.root = null;
+}

.ACM-Templates/Data-Structures/splay/remark.tex

@@ -0,0 +1 @@
+splay树模板，分为两个部分，第一个部分为splay当作平衡树使用，第二个部分为splay当作线段树维护序列使用。注意在每次Insert之后都要splay其返回值到根

.ACM-Templates/Data-Structures/suffix-array/main.cpp

@@ -0,0 +1,37 @@
+// string is 1-base, sa is 1-base
+int w[MAXM];
+inline void Sort(int a[], int ret[], int n, int m = MAXM - 1) {
+	for (int i = 0; i <= m; i++) w[i] = 0;
+	for (int i = 1; i <= n; i++) w[a[i]]++;
+	for (int i = 1; i <= m; i++) w[i] += w[i - 1];
+	for (int i = n; i >= 1; i--) ret[w[a[i]]--] = i;
+}
+int wa[MAXN], wb[MAXN], tmp[MAXN];
+inline void getSA(int ch[], int sa[], int n) {
+	int *x = wa, *y = wb;
+	for (int i = 1; i <= n; i++) x[i] = ch[i];
+	Sort(ch, sa, n);
+	for (int j = 1, p = 1, m = MAXN - 1; p < n; m = p, j <<= 1) {
+		p = 0;
+		for (int i = n - j + 1; i <= n; i++) y[++p] = i;
+		for (int i = 1; i <= n; i++) if (sa[i] > j) y[++p] = sa[i] - j;
+		for (int i = 1; i <= n; i++) tmp[i] = x[y[i]];
+		Sort(tmp, sa, n, m);
+		for (int i = 1; i <= n; i++) sa[i] = y[sa[i]];
+		swap(x, y); x[sa[1]] = p = 1;
+		for (int i = 2; i <= n; i++) {
+			if (y[sa[i]] == y[sa[i - 1]] && y[sa[i] + j] == y[sa[i - 1] + j]) x[sa[i]] = p;
+			else x[sa[i]] = ++p;
+		}
+	}
+	sa[0] = n + 1; // for calculate height.
+}
+int rank[MAXN];
+inline void getHeight(int ch[], int sa[], int height[], int n) {
+	for (int i = 1; i <= n; i++) rank[sa[i]] = i;
+	for (int i = 1, t = 0; i <= n; i++) {
+		if (t > 0) t--;
+		while (ch[i + t] == ch[sa[rank[i] - 1] + t]) t++;
+		height[rank[i]] = t;
+	}
+}

.ACM-Templates/Data-Structures/suffix-array/remark.tex

@@ -0,0 +1,2 @@
+倍增算法构建后缀数组，时间复杂度为 $O(n \lg n)$. 注意字符串和最终的后缀数组均从 1 开始编号. 要保证字符串中都为大于0的字符, 而且字符串的第n+1位应该为0.
+

.ACM-Templates/Data-Structures/suffix-automaton/main.cpp

@@ -0,0 +1,29 @@
+struct Node {
+	Node *next[26], *par; int val, end; // 26 is volatile
+}POOL[MAXN << 1], *data, *root, *last; //Note that the size of POOL should be doubled.
+inline void Add(int x) {
+	Node *p = last, *np = data++;
+	np->val = p->val + 1; np->end = true;
+	while (p && !p->next[x]) 
+		p->next[x] = np, p = p->par;
+	if (p == 0) {
+		np->par = root;
+	} else {
+		Node *q = p->next[x];
+		if (q->val == p->val + 1) {
+			np->par = q;
+		} else {
+			Node *nq = data++;
+			nq->val = p->val + 1;
+			memcpy(nq->next, q->next, sizeof q->next);
+			nq->par = q->par;
+			np->par = q->par = nq;
+			while (p && p->next[x] == q)
+				p->next[x] = nq, p = p->par;
+		}
+	}
+	last = np;
+}
+void Clear(void) {
+	data = POOL; last = root = data++;
+}

.ACM-Templates/Data-Structures/suffix-automaton/remark.tex

@@ -0,0 +1 @@
+构建后缀自动机。统计一个串出现的次数时，只需统计其节点所对应的子树中，end为true的节点的个数即可。将所有节点按val值从小到大排序后即可得到parent树由根开始的BFS序。

.ACM-Templates/Graph/Hopcroft-Karp/main.cpp

@@ -0,0 +1,43 @@
+const int MAXN = 10000;
+int match[MAXN];  
+int levelx[MAXN], levely[MAXN], link[MAXN]; 
+int d[MAXN]; 
+inline bool Bfs(void) { 
+	int head = 1, tail = 0; 
+	memset(levely, 0, sizeof levely); 
+	for (int i = 1; i <= n; i++) { 
+		if (!match[i]) d[++tail] = i; 
+		levelx[i] = 0; 
+	} 
+	bool ret = false; 
+	while (head <= tail) { 
+		int now = d[head++]; 
+		for (Edge *p = a[now]; p; p = p->next) if (levely[p->y] == 0) { 
+			levely[p->y] = levelx[now] + 1; 
+			if (link[p->y] == 0) ret = true; 
+			else levelx[link[p->y]] = levely[p->y] + 1, d[++tail] = link[p->y]; 
+		} 
+	} 
+	return ret; 
+} 
+bool Find(int u) { 
+	for (Edge *p = a[u]; p; p = p->next) if (levely[p->y] == levelx[u] + 1) { 
+		levely[p->y] = 0; 
+		if (link[p->y] == 0 || Find(link[p->y])) { 
+			match[u] = p->y; link[p->y] = u; 
+			return true; 
+		} 
+	} 
+	return false; 
+} 
+inline void Match(void) { 
+	while (Bfs()) 
+		for (int i = 1; i <= n; i++) 
+			if (!match[i]) Find(i); 
+} 
+inline void clear(void) {
+	memset(match, 0, sizeof match);
+	memset(link, 0, sizeof link);
+	memset(levelx, 0, sizeof levelx);
+	memset(levely, 0, sizeof levely);
+}

.ACM-Templates/Graph/Hopcroft-Karp/remark.tex

@@ -0,0 +1 @@
+Hopcroft 算法，其中二分图左边点的编号为1-n，右边点编号为n+1 - N，match[i]表示左边的i匹配上的右边的点，0无匹配，link[i]表示右边的i匹配上的左边的点，0无匹配

.ACM-Templates/Graph/KM-algorithm/main.cpp

@@ -0,0 +1,56 @@
+const int MAXN = 210;
+int nx, ny; // 左边的点标号从1到nx，右边点标号从1到ny
+long long inf, cost[MAXN][MAXN], fx[MAXN], fy[MAXN], dist[MAXN]; //权值若为long long的话，只需改动此行即可
+int used[MAXN], maty[MAXN], which[MAXN];
+inline void AddEdge(int x, int y, int z) {
+	cost[x][y] = z;
+}
+pair<int, long long> KM(void) {
+	for (int x = 1; x <= nx; x++) {
+		int y0 = 0; maty[0] = x;
+		for (int y = 0; y <= ny; y++) { dist[y] = inf + 1; used[y] = false; }
+		do {
+			used[y0] = true;
+			int x0 = maty[y0], y1;
+			long long delta = inf + 1;
+			for (int y = 1; y <= ny; y++) if (!used[y]) {
+				long long curdist = cost[x0][y] - fx[x0] - fy[y];
+				if (curdist < dist[y]) {
+					dist[y] = curdist;
+					which[y] = y0;
+				}
+				if (dist[y] < delta) {
+					delta = dist[y];
+					y1 = y;
+				}
+			}
+			for (int y = 0; y <= ny; y++) if (used[y]) {
+				fx[maty[y]] += delta;
+				fy[y] -= delta;
+			} else dist[y] -= delta;
+			y0 = y1;
+		} while (maty[y0] != 0);
+		do {
+			int y1 = which[y0];
+			maty[y0] = maty[y1];
+			y0 = y1;
+		} while (y0);
+	}
+	long long ret = 0;
+	int npair = 0;
+	for (int y = 1; y <= ny; y++) {
+		int x = maty[y];
+		if (cost[x][y] < inf) {
+			ret += cost[x][y];
+			npair++;
+		}
+	}
+	return make_pair(npair, ret);
+}
+inline void clear(void) {
+	memset(fx, 0x9f, sizeof fx);
+	memset(fy, 0x9f, sizeof fy);
+	memset(cost, 0x3f, sizeof cost);
+	memset(maty, 0, sizeof maty);
+	inf = cost[0][0];
+}

.ACM-Templates/Graph/KM-algorithm/remark.tex

@@ -0,0 +1,3 @@
+KM算法求最小权完美匹配（注意不是全局最大权匹配），如果两边点数不相等或者非完美二分图，则可以添加虚拟点和权值为inf的边
+
+若求最大权匹配，则将边权取反即可

.ACM-Templates/Graph/Maximum-Flow/main.cpp

@@ -0,0 +1,42 @@
+struct Edge {
+	int y, f; Edge *next, *opt;
+}*a[MAXN], DATA[MAXM << 1], *data = DATA;
+inline void Add(int x, int y, int c) {
+	Edge *tmp = data++;
+	tmp->y = y; tmp->f = c, tmp->next = a[x]; a[x] = tmp;
+	tmp = data++; tmp->y = x; tmp->f = 0; tmp->next = a[y]; a[y] = tmp;
+	a[x]->opt = a[y]; a[y]->opt = a[x];
+}
+int n, m, vs, vt, L;
+int level[MAXN], d[MAXN];
+inline bool Bfs(void) {
+	memset(level, -1, sizeof level);
+	d[1] = vs; level[vs] = 0;
+	int head = 1, tail = 1;
+	while (head <= tail) {
+		int now = d[head++];
+		e[now] = a[now];
+		for (Edge *p = a[now]; p; p = p->next) if (p->f > 0&& level[p->y] == -1)
+			level[d[++tail] = p->y] = level[now] + 1;
+	}
+	return level[vt] != -1;
+}
+inline int Extend(int u, int sum) {
+	if (u == vt) return sum;
+	int r = 0, t;
+	for (Edge *p = e[u]; p && r < sum; p = p->next) if (level[p->y] == level[u] + 1 && p->f > 0) {
+		t = std::min(sum - r, p->f);
+		t = Extend(p->y, t);
+		p->f -= t, p->opt->f += t, r += t;
+		e[u] = p;
+	}
+	if (!r) level[u] = -1;
+	return r;
+}
+inline int Dinic(void) {
+	int r = 0, t;
+	while (Bfs()) {
+		while ((t = Extend(vs, inf))) r += t;
+	}
+	return r;
+}

.ACM-Templates/Graph/Maximum-Flow/remark.tex

@@ -0,0 +1 @@
+Dinic求解最大流，有当前弧优化

.ACM-Templates/Graph/Minimum-cost-Flow/main.cpp

@@ -0,0 +1,47 @@
+struct Edge {
+    int y, f, c; Edge *next, *opt;
+    Edge(int y, int f, int c, Edge *next):y(y), f(f), c(c), next(next){}
+}*a[MAXN];
+inline void AddEdge(int x, int y, int f, int c) {
+    a[x] = new Edge(y, f, c, a[x]);
+    a[y] = new Edge(x, 0, -c, a[y]);
+    a[x]->opt = a[y]; a[y]->opt = a[x];
+}
+int d[MAXN], vis[MAXN], dis[MAXN]; Edge *path[MAXN];
+inline pair<int, int> Spfa(void) {
+    int Flow = 0, Cost = 0;
+    while (true) {
+        memset(vis, 0, sizeof vis);
+        memset(dis, 0x7f, sizeof dis);
+        memset(path, 0, sizeof path);
+        int head = 1, tail = 1, sum = 1; d[1] = vs; vis[vs] = true; dis[vs] = 0;
+        while (sum) {
+            int now = d[head++]; if (head == MAXN) head = 1; sum--;
+            for (Edge *p = a[now]; p; p = p->next) if (p->f > 0 && dis[p->y] > dis[now] + p->c) {
+                dis[p->y] = dis[now] + p->c;
+                path[p->y] = p;
+                if (!vis[p->y]) {
+                    ++tail; if (tail == MAXN) tail = 1; sum++;
+                    d[tail] = p->y;
+                    vis[p->y] = true;
+                }
+            }
+            vis[now] = false;
+        }
+        if (dis[vt] == dis[0]) return make_pair(Flow, Cost);
+        int tmp = vt, Min = ~0U>>1;
+        while (path[tmp]) {
+            Min = min(Min, path[tmp]->f);
+            tmp = path[tmp]->opt->y;
+        }
+        Flow += Min;
+        tmp = vt;
+        while (path[tmp]) {
+            path[tmp]->f -= Min;
+            path[tmp]->opt->f += Min;
+            Cost += Min * path[tmp]->c;
+            tmp = path[tmp]->opt->y;
+        }
+    }
+    return make_pair(Flow, Cost);
+}

.ACM-Templates/Graph/Minimum-cost-Flow/remark.tex

@@ -0,0 +1 @@
+spfa求解最小费用最大流

.ACM-Templates/Graph/get-bridge/main.cpp

@@ -1,6 +1,5 @@
 int dfs[MAXN], low[MAXN], tim = 0;
-inline void Dfs(int u) 
-{
+inline void Dfs(int u) {
 	dfs[u] = low[u] = ++tim;
 	for (Edge *p = a[u]; p; p = p->next) if (p->flag) {
 		if (!dfs[p->y]) {
@@ -13,4 +12,3 @@ inline void Dfs(int u)
 		p->bridge = p->opt->bridge = true;
 	}
 }
-

.ACM-Templates/Graph/get-bridge/remark.tex

@@ -0,0 +1 @@
+Tarjan 算法求一个无向图的桥

.ACM-Templates/Graph/最小树形图/main.cpp

@@ -0,0 +1,55 @@
+typedef pair<int, int> PII;
+PII tmp[MAXN];
+int dis[MAXN][MAXN]; int map[MAXN][MAXN];
+int vis[MAXN], q[MAXN], inp[MAXN];
+int Dfs(int u) {
+	vis[u] = true; int ret = 1;
+	for (int i = 1; i <= n; i++) if (map[u][i] && !vis[i]) ret += Dfs(i);
+	return ret;
+}
+int Zhuliu(void) {
+	memset(map, 0, sizeof map);
+	memset(vis, 0, sizeof vis);
+	if (Dfs(1) != n) return -1;
+	int done = 0;
+	while (true) {
+		memset(vis, 0, sizeof vis); memset(inp, 0, sizeof inp);
+		int Ans = 0;
+		for (int i = 1; i <= n; i++) vis[i] = i;
+		for (int i = 2; i <= n; i++) {
+			tmp[i] = PII(1000000000, 0);
+			for (int j = 1; j <= n; j++) if (map[j][i] == 1) tmp[i] = min(tmp[i], PII(dis[j][i], j));
+			inp[tmp[i].second]++; if (tmp[i].second) Ans += tmp[i].first;
+		}
+		int head = 1, tail = 0;
+		for (int i = 1; i <= n; i++) if (!inp[i]) q[++tail] = i;
+		while (head <= tail) {
+			int now = q[head++];
+			if (!--inp[tmp[now].second]) q[++tail] = tmp[now].second;
+		}
+		bool ok = true;
+		for (int i = 1, t; i <= n; i++) if (inp[i] > 0) {
+			t = i; ok = false;
+			do {
+				inp[t] = -i;
+				t = tmp[t].second;
+				vis[t] = i;
+			} while (t != i);
+		}
+		if (ok) return Ans + done;
+		for (int i = 1; i <= n; i++) if (inp[i] < 0) {
+			done += tmp[i].first;
+			for (int j = 1; j <= n; j++) if (vis[j] != vis[i]) {
+				if (map[i][j]) {
+					map[vis[i]][vis[j]] = 1;
+					dis[vis[i]][vis[j]] = min(dis[vis[i]][vis[j]], dis[i][j]);
+				}
+				if (map[j][i]) {
+					dis[vis[j]][vis[i]] = min(dis[vis[j]][vis[i]], dis[j][i] - tmp[i].first);
+					map[vis[j]][vis[i]] = 1;
+				}
+			}
+			if (vis[i] != i) for (int j = 1; j <= n; j++) map[i][j] = map[j][i] = 0;
+		}
+	}
+}

.ACM-Templates/Graph/最小树形图/remark.tex

@@ -0,0 +1,3 @@
+最小树形图算法，使用邻接矩阵存图，返回-1表示无最小树形图。否则，tmp[i].second表示i的入边
+map[i][j] = 1 表示表示存在一条从i到j的有向边
+dis[i][j] 表示从i到j边的长度

.ACM-Templates/Math/geometry/main.cpp

@@ -0,0 +1,424 @@
+#include <cstdio>
+#include <cmath>
+#include <vector>
+#include <deque>
+#include <algorithm>
+const double eps = 1e-13;
+const double pi = 3.14159265358979324;
+struct point2 {
+	double x, y;
+	point2& operator += (point2 a) { x+=a.x, y+=a.y; return *this; }
+	point2& operator -= (point2 a) { x-=a.x, y-=a.y; return *this; }
+	point2& operator *= (double a) { x*=a, y*=a; return *this; }
+	point2& operator /= (double a) { x/=a, y/=a; return *this; }
+};
+point2 operator + (point2 a, point2 b) { point2 c(a); c += b; return c; }
+point2 operator - (point2 a, point2 b) { point2 c(a); c -= b; return c; }
+point2 operator * (point2 a, double b) { point2 c(a); c *= b; return c; }
+point2 operator * (double a, point2 b) { point2 c(b); c *= a; return c; }
+point2 operator / (point2 a, double b) { point2 c(a); c /= b; return c; }
+double operator * (point2 a, point2 b) { return a.x*b.x+a.y*b.y; }
+double operator % (point2 a, point2 b) { return a.x*b.y-a.y*b.x; }
+double dis(point2 a) { return sqrt(a.x * a.x + a.y * a.y); }
+double arg(point2 a) { return atan2(a.y, a.x); }
+point2 rotate(point2 a, double th) {
+	point2 b;
+	b.x = a.x * cos(th) - a.y * sin(th);
+	b.y = a.x * sin(th) + a.y * cos(th);
+	return b;
+}
+int parallel(point2 a, point2 b) {
+	return a * a < eps * eps || b * b < eps * eps
+		|| (a % b) * (a % b) / ((a * a) * (b * b)) < eps * eps;
+}
+int perpend(point2 a, point2 b) {
+	return a * a < eps * eps || b * b < eps * eps
+		|| (a * b) * (a * b) / ((a * a) * (b * b)) < eps * eps;
+}
+struct line2 {
+	point2 a, s;
+};
+struct circle2 {
+	point2 a;
+	double r;
+};
+double point_line_dis(point2 a, line2 b, point2 *res = NULL) {
+	point2 p;
+	p = b.a + ((a - b.a) * b.s) / (b.s * b.s) * b.s;
+	if (res != NULL) *res = p;
+	return dis(a - p);
+}
+int line_line_cross(line2 a, line2 b, point2 *res = NULL) {
+	if (parallel(a.s, b.s))
+		if (parallel(b.a - a.a, a.s))
+			return -1;
+		else
+			return 0;
+	double k1 = (b.a - a.a) % b.s / (a.s % b.s);
+	if (res != NULL) *res = a.a + k1 * a.s;
+	return 1;
+}
+int segment_segment_cross(line2 a, line2 b, point2 *res = NULL) {
+	if (a.s * a.s < eps * eps && b.s * b.s < eps * eps)
+		if ((b.a - a.a) * (b.a - a.a) < eps * eps) {
+			if (res != NULL) *res = a.a;
+			return 1;
+		} else
+			return 0;
+	if (parallel(a.s, b.s) && parallel(b.a - a.a, a.s) && parallel(a.a - b.a, b.s)) {
+		double y1, y2, y3, y4;
+		point2 y1p = a.a, y2p = a.a + a.s, y3p = b.a, y4p = b.a + b.s;
+		if (std::abs(a.s.x) < std::abs(a.s.y) || std::abs(b.s.x) < std::abs(b.s.y))
+			y1 = y1p.y, y2 = y2p.y, y3 = y3p.y, y4 = y4p.y;
+		else
+			y1 = y1p.x, y2 = y2p.x, y3 = y3p.x, y4 = y4p.x;
+		if (y1 > y2) std::swap(y1, y2), std::swap(y1p, y2p);
+		if (y3 > y4) std::swap(y3, y4), std::swap(y3p, y4p);
+		if (y2 - y1 < y4 - y3) std::swap(y1, y3), std::swap(y1p, y3p), std::swap(y2, y4), std::swap(y2p, y4p);
+		if (y3 > y2 + (y2 - y1) * eps || y4 < y1 - (y2 - y1) * eps)
+			return 0;
+		else if (fabs(y3 - y2) < (y2 - y1) * eps || fabs(y3 - y4) < eps) {
+			if (res != NULL) *res = y3p;
+			return 1;
+		} else if (fabs(y4 - y1) < (y2 - y1) * eps || fabs(y1 - y2) < eps) {
+			if (res != NULL) *res = y1p;
+			return 1;
+		} else
+			return -1;
+	} else {
+		double k1 = (b.a - a.a) % a.s, k2 = (b.a + b.s - a.a) % a.s;
+		k1 /= a.s * a.s, k2 /= a.s * a.s;
+		double k3 = (a.a - b.a) % b.s, k4 = (a.a + a.s - b.a) % b.s;
+		k3 /= b.s * b.s, k4 /= b.s * b.s;
+		int ret = (k1 < eps && k2 > -eps || k1 > -eps && k2 < eps)
+			&& (k3 < eps && k4 > -eps || k3 > -eps && k4 < eps);
+		if (ret) line_line_cross(a, b, res);
+		return ret;
+	}
+}
+int line_circle_cross(line2 a, circle2 b, point2 *res1 = NULL, point2 *res2 = NULL) {
+	point2 p;
+	double d = point_line_dis(b.a, a, &p);
+	if (d / b.r > 1 + eps)
+		return 0;
+	else if (d / b.r > 1 - eps) {
+		if (res1 != NULL) *res1 = p;
+		return 1;
+	} else {
+		d = sqrt(b.r * b.r - d * d) / dis(a.s);
+		if (res1 != NULL) *res1 = p + d * a.s;
+		if (res2 != NULL) *res2 = p - d * a.s;
+		return 2;
+	}
+}
+int circle_circle_cross(circle2 a, circle2 b, point2 *res1 = NULL, point2 *res2 = NULL) {
+	double d = dis(a.a - b.a);
+	point2 u = (b.a - a.a) / d;
+	if (d / (a.r + b.r) > 1 + eps)
+		return 0;
+	else if (d / (a.r + b.r) > 1 - eps) {
+		if (res1 != NULL) *res1 = a.a + u * a.r;
+		return 1;
+	} else if ((d - fabs(a.r - b.r)) / (a.r + b.r) > eps) {
+		double th = acos((a.r * a.r + d * d - b.r * b.r) / (2 * a.r * d));
+		if (res1 != NULL) *res1 = a.a + rotate(u * a.r, th);
+		if (res2 != NULL) *res2 = a.a + rotate(u * a.r, -th);
+		return 2;
+	} else if ((d - fabs(a.r - b.r)) / (a.r + b.r) > -eps) {
+		if (a.r / b.r < 1 - eps) {
+			if (res1 != NULL) *res1 = b.a - u * b.r;
+			return 1;
+		} else if (a.r / b.r > 1 + eps) {
+			if (res1 != NULL) *res1 = a.a + u * a.r;
+			return 1;
+		} else return -1;
+	} else
+		return 0;
+}
+int point_circle_tangent(point2 a, circle2 b, point2 *res1 = NULL, point2 *res2 = NULL) {
+	double d = dis(a - b.a);
+	point2 u = (a - b.a) / d;
+	if (d / b.r > 1 + eps) {
+		double th = acos(b.r / d);
+		if (res1 != NULL) *res1 = b.a + rotate(u * b.r, th);
+		if (res2 != NULL) *res2 = b.a + rotate(u * b.r, -th);
+		return 2;
+	} else if (d / b.r > 1 - eps) {
+		if (res1 != NULL) *res1 = a;
+		return 1;
+	} else
+		return 0;
+}
+int circle_circle_tangent(circle2 a, circle2 b, line2 *reso1 = NULL, line2 *reso2 = NULL, line2 *resi1 = NULL, line2 *resi2 = NULL) {
+	double d = dis(a.a - b.a);
+	point2 u = (b.a - a.a) / d;
+	int cnt = 0;
+	if ((d - fabs(a.r - b.r)) / (a.r + b.r) > eps) {
+		double th = acos((a.r - b.r) / d);
+		if (reso1 != NULL) {
+			reso1->a = a.a + rotate(u * a.r, th);
+			reso1->s = b.a + rotate(u * b.r, th) - reso1->a;
+		}
+		if (reso2 != NULL) {
+			reso2->a = a.a + rotate(u * a.r, -th);
+			reso2->s = b.a + rotate(u * b.r, -th) - reso2->a;
+		}
+		cnt += 2;
+	} else if ((d - fabs(a.r - b.r)) / (a.r + b.r) > -eps) {
+		if (a.r / b.r < 1 - eps) {
+			if (reso1 != NULL) {
+				reso1->a = b.a - u * b.r;
+				reso1->s = rotate(u, pi / 2);
+			}
+			cnt++;
+		} else if (a.r / b.r > 1 + eps) {
+			if (reso1 != NULL) {
+				reso1->a = a.a + u * a.r;
+				reso1->s = rotate(u, pi / 2);
+			}
+			cnt++;
+		} else return -1;
+	}
+	if (d / (a.r + b.r) > 1 + eps) {
+		double th = acos((a.r + b.r) / d);
+		if (resi1 != NULL) {
+			resi1->a = a.a + rotate(u * a.r, th);
+			resi1->s = b.a - rotate(u * b.r, th) - resi1->a;
+		}
+		if (resi2 != NULL) {
+			resi2->a = a.a + rotate(u * a.r, -th);
+			resi2->s = b.a - rotate(u * b.r, -th) - resi2->a;
+		}
+		cnt += 2;
+	} else if (d / (a.r + b.r) > 1 - eps) {
+		if (resi1 != NULL) {
+			resi1->a = a.a + u * a.r;
+			resi1->s = rotate(u, pi / 2);
+		}
+		cnt++;
+	}
+	return cnt;
+}
+typedef std::vector<point2> convex2;
+double area(const convex2 &a) {
+	double s = 0;
+	for (int i = 0; i < a.size(); i++)
+		s += a[i] % a[(i+1)%a.size()];
+	return fabs(s)/2;
+}
+int point_convex_inside(point2 a, const convex2 &b) {
+	double sum = 0;
+	for (int i = 0; i < b.size(); i++)
+		sum += fabs((b[i] - a) % (b[(i+1)%b.size()] - a));
+	return fabs(sum / (2*area(b)) - 1) < eps;
+}
+int line_convex_cross(line2 a, const convex2 &b, point2 *res1 = NULL, point2 *res2 = NULL) {
+	int cnt = 0;
+	for (int i = 0; i < b.size(); i++) {
+		line2 ltmp; point2 ptmp;
+		ltmp.a = b[i], ltmp.s = b[(i+1)%b.size()] - b[i];
+		int flag = line_line_cross(a, ltmp, &ptmp);
+		if (flag == -1) return -1;
+		if (flag == 0) continue;
+		double k = (ptmp - ltmp.a) * ltmp.s / (ltmp.s * ltmp.s);
+		if (k < eps || k > 1+eps) continue;
+		if (cnt == 0 && res1 != NULL) *res1 = ptmp;
+		else if (cnt == 1 && res2 != NULL) *res2 = ptmp;
+		cnt++;
+	}
+	return cnt;
+}
+int convex_gen_cmp(point2 a, point2 b) {
+	return a.y < b.y - eps || fabs(a.y - b.y) < eps && a.x < b.x - eps;
+}
+int convex_gen(const convex2 &a, convex2 &b) {
+	std::deque<point2> q;
+	convex2 t(a);
+	std::sort(t.begin(), t.end(), convex_gen_cmp);
+	q.push_back(t[0]), q.push_back(t[1]);
+	for (int i = 2; i < t.size(); i++) {
+		while (q.size() > 1) {
+			point2 p1 = t[i]-q[q.size()-1], p2 = q[q.size()-1]-q[q.size()-2];
+			if (p1 % p2 > 0 || parallel(p1, p2)) q.pop_back();
+			else break;
+		}
+		q.push_back(t[i]);
+	}
+	int pretop = q.size();
+	for (int i = t.size()-1; i >= 0; i--) {
+		while (q.size() > pretop) {
+			point2 p1 = t[i]-q[q.size()-1], p2 = q[q.size()-1]-q[q.size()-2];
+			if (p1 % p2 > 0 || parallel(p1, p2)) q.pop_back();
+			else break;
+		}
+		q.push_back(t[i]);
+	}
+	q.pop_back();
+	if (q.size() < 3) {
+		b.clear();
+		return 0;
+	}
+	b.clear();
+	for (int i = 0; i < q.size(); i++) b.push_back(q[i]);
+	return 1;
+}
+int halfplane_cross_cmp(line2 a, line2 b) {
+	double c1 = arg(a.s), c2 = arg(b.s);
+	return c1 < c2-eps || fabs(c1-c2) < eps && b.s % (a.a - b.a) / dis(b.s) > eps;
+}
+int halfplane_cross(const std::vector<line2> &a, convex2 &b) {
+	std::vector<line2> t(a);
+	std::sort(t.begin(), t.end(), halfplane_cross_cmp);
+	int j = 0;
+	for (int i = 0; i < t.size(); i++)
+		if (!i || arg(t[i].s) > arg(t[i-1].s) + eps) t[j++] = t[i];
+	if (j > 0 && arg(t[j].s) > arg(t[0].s) + 2*pi - eps) j--;
+	t.resize(j);
+	std::deque<line2> q;
+	q.push_back(t[0]), q.push_back(t[1]);
+	point2 p;
+	for (int i = 2, k = 0; i < t.size(); i++) {
+		for (; k < q.size() && t[i].s % q[k].s > 0; k++);
+		point2 s1 = q[q.size()-1].s, s2 = q[0].s;
+		if (k > 0 && k < q.size() && s1 % s2 > 0 && !parallel(s1, s2)) {
+			line_line_cross(q[k], q[k-1], &p);
+			double r1 = t[i].s % (p - t[i].a) / dis(t[i].s);
+			line_line_cross(q[0], q[q.size()-1], &p);
+			double r2 = t[i].s % (p - t[i].a) / dis(t[i].s);
+			if (r1 < eps && r2 < eps) {
+				b.clear();
+				return 0;
+			} else if (r1 > -eps && r2 > -eps)
+				continue;
+		}
+		while (q.size() > 1) {
+			line_line_cross(q[q.size()-1], q[q.size()-2], &p);
+			if (t[i].s % (p - t[i].a) / dis(t[i].s) < eps) {
+				q.pop_back();
+				if (k == q.size()) k--;
+			} else break;
+		}
+		while (q.size() > 1) {
+			line_line_cross(q[0], q[1], &p);
+			if (t[i].s % (p - t[i].a) / dis(t[i].s) < eps) {
+				q.pop_front();
+				k--; if (k < 0) k = 0;
+			} else break;
+		}
+		q.push_back(t[i]);
+	}
+	b.clear();
+	for (int i = 0; i < q.size(); i++) {
+		line_line_cross(q[i], q[(i+1)%q.size()], &p);
+		b.push_back(p);
+	}
+	return 1;
+}
+struct point3 {
+	double x, y, z;
+	point3& operator += (point3 a) { x+=a.x,y+=a.y,z+=a.z; return *this; }
+	point3& operator -= (point3 a) { x-=a.x,y-=a.y,z-=a.z; return *this; }
+	point3& operator *= (double a) { x*=a, y*=a, z*=a; return *this; }
+	point3& operator /= (double a) { x/=a, y/=a, z/=a; return *this; }
+};
+point3 operator + (point3 a, point3 b) { point3 c(a); c += b; return c; }
+point3 operator - (point3 a, point3 b) { point3 c(a); c -= b; return c; }
+point3 operator * (point3 a, double b) { point3 c(a); c *= b; return c; }
+point3 operator * (double a, point3 b) { point3 c(b); c *= a; return c; }
+point3 operator / (point3 a, double b) { point3 c(a); c /= b; return c; }
+double operator * (point3 a, point3 b) { return a.x*b.x+a.y*b.y+a.z*b.z; }
+point3 operator % (point3 a, point3 b) {
+	point3 c;
+	c.x = a.y * b.z - a.z * b.y;
+	c.y = a.z * b.x - a.x * b.z;
+	c.z = a.x * b.y - a.y * b.x;
+	return c;
+}
+double dis(point3 a) { return sqrt(a.x * a.x + a.y * a.y + a.z * a.z); }
+int parallel(point3 a, point3 b) {
+	return a * a < eps * eps || b * b < eps * eps
+		|| (a % b) * (a % b) / ((a * a) * (b * b)) < eps * eps;
+}
+int perpend(point3 a, point3 b) {
+	return a * a < eps * eps || b * b < eps * eps
+		|| (a * b) * (a * b) / ((a * a) * (b * b)) < eps * eps;
+}
+struct line3 {
+	point3 a, s;
+};
+struct face3 {
+	point3 a, n;
+};
+struct circle3 {
+	point3 a;
+	double r;
+};
+point2 project(point3 a, face3 b, point3 xs) {
+	point3 ys; ys = b.n % xs;
+	point2 c; c.x = ((a - b.a) * xs) / dis(xs), c.y = ((a - b.a) * ys) / dis(ys);
+	return c;
+}
+line2 project(line3 a, face3 b, point3 xs) {
+	line2 c; c.a = project(a.a, b, xs), c.s = project(a.a + a.s, b, xs) - c.a;
+	return c;
+}
+point3 revproject(point2 a, face3 b, point3 xs) {
+	point3 ys; ys = b.n % xs;
+	return a.x * xs / dis(xs) + a.y * ys / dis(ys) + b.a;
+}
+line3 revproject(line2 a, face3 b, point3 xs) {
+	line3 c; c.a = revproject(a.a, b, xs), c.s = revproject(a.a + a.s, b, xs) - c.a;
+	return c;
+}
+double point_line_dis(point3 a, line3 b, point3 *res = NULL) {
+	point3 p;
+	p = b.a + ((a - b.a) * b.s) / (b.s * b.s) * b.s;
+	if (res != NULL) *res = p;
+	return dis(a - p);
+}
+double point_face_dis(point3 a, face3 b, point3 *res = NULL) {
+	point3 p;
+	p = ((a - b.a) * b.n) / (b.n * b.n) * b.n;
+	if (res != NULL) *res = a - p;
+	return dis(p);
+}
+double line_line_dis(line3 a, line3 b, point3 *resa = NULL, point3 *resb = NULL) {
+	point3 p;
+	if (parallel(a.s, b.s)) {
+		double d = point_line_dis(a.a, b, &p);
+		if (resa != NULL) *resa = a.a;
+		if (resb != NULL) *resb = p;
+		return d;
+	}
+	point3 n = a.s % b.s;
+	face3 f; f.a = b.a, f.n = n;
+	double d = point_face_dis(a.a, f, &p);
+	double k1 = ((b.a - p) % b.s) * n / (n * n);
+	if (resb != NULL) *resb = p + k1 * a.s;
+	if (resa != NULL) *resa = a.a + k1 * a.s;
+	return d;
+}
+int line_face_cross(line3 a, face3 b, point3 *res = NULL) {
+	if (perpend(a.s, b.n))
+		if (perpend(b.a - a.a, b.n))
+			return -1;
+		else
+			return 0;
+	double k = (b.a - a.a) * b.n / (a.s * b.n);
+	if (res != NULL) *res = a.a + k * a.s;
+	return 1;
+}
+int face_face_cross(face3 a, face3 b, line3 *res = NULL) {
+	if (parallel(a.n, b.n))
+		if (perpend(b.a - a.a, a.n))
+			return -1;
+		else
+			return 0;
+	point3 s = a.n % b.n;
+	point3 p;
+	line3 t; t.a = a.a, t.s = a.n % s;
+	line_face_cross(t, b, &p);
+	if (res != NULL)
+		res->a = p, res->s = s;
+	return 1;
+}

.ACM-Templates/Math/geometry/remark.tex

@@ -0,0 +1,9 @@
+计算几何模板，包含：
+\begin{itemize}
+\item point2/3 点/向量，方法：向量加+减-，数乘*，点乘*，叉乘%，模dis，辐角arg（2维），旋转rotate（2维），平行parallel，垂直perpend，三维点根据所在面的法向量和一个面上向量投影为二维点以及反投影
+\item line2/3 线，表示为起点+方向向量，方法：点线距和垂足，线线交、线段线段交（线段交可退化，2维），线线距和最近点对（3维，若平行返回任意一对），线的投影（同点）
+\item face3 面，表示为起点+法向量，方法：点面距和垂足，线面交，面面交
+\item circle2 圆，表示为圆心+半径，方法：线圆交，圆圆交，点到圆的切点，圆与圆的公切线（若公切线有两个切点，返回直线两端点为两切点，否则返回直线端点为切点）
+\item convex2 凸包，方法：面积，点不严格在凸包内，线与凸包交，水平序Graham求凸包，半平面交（半平面为传入向量左侧平面）
+\end{itemize}
+说明：对需要求点/线的，点/线作为指针传入，若为NULL表示不需要（默认）；对需要求点/线个数的，返回值为-1表示有无穷多个。

.ACM-Templates/Math/numbertheory/main.cpp

@@ -0,0 +1,96 @@
+#include <cstdlib>
+#include <cmath>
+#include <map>
+#include <vector>
+#include <algorithm>
+void ext_gcd(int a, int b, int &x, int &y) {
+	if (!b) x = 1, y = 0;
+	else if (!a) x = 1, y = -1;
+	else if (a > b) ext_gcd(a % b, b, x, y), y += a / b * x;
+	else ext_gcd(a, b % a, x, y), x += b / a * y;
+}
+long long flsum_t (long long a, long long b, long long c, long long n) {
+	if (n < 0) return 0;
+	if (c < 0) a = -a, b = -b, c = -c;
+	n++; long long res = 0;
+	if (a < 0 || a >= c) {
+		long long ra = (a % c + c) % c;
+		long long k = (a - ra) / c;
+		res += k * n * (n - 1) / 2;
+		a = ra;
+	}
+	if (b < 0 || b >= c) {
+		long long rb = (b % c + c) % c;
+		long long k = (b - rb) / c;
+		res += k * n;
+		b = rb;
+	}
+	if (a * n + b < c) return res;
+	else return res + flsum_t(c, (a * n + b) % c, a, (a * n + b) / c - 1);
+}
+long long flsum (long long a, long long b, long long c, long long st, long long ed) {
+	return flsum_t(a, b, c, ed) - flsum_t(a, b, c, st - 1);
+}
+int power(int n, int k, int r) {
+	int t = n, s = 1;
+	for (; k; k >>= 1, t = 1LL * t * t % r)
+		if (k & 1) s = 1LL * s * t % r;
+	return s;
+}
+int millerrabin(int x, int tester) {
+	int k = x-1; for (; !(k & 1); k >>= 1);
+	int y = power(tester, k, x);
+	if (y == 1) return 1;
+	for (; k < x-1; k <<= 1, y = 1LL * y * y % x)
+		if (y == x-1) return 1;
+	return 0;
+}
+int isprime(int x) {
+	if (x == 2 || x == 7 || x == 61) return 1;
+	return millerrabin(x, 2) && millerrabin(x, 7) && millerrabin(x, 61);
+}
+int rho_f(int x, int c, int p) {
+	return (1LL * x * x + c) % p;
+}
+int rho(int n) {
+	int c = rand() % (n-1) + 1, x = 2, y = x, d = 1;
+	while (d == 1) {
+		x = rho_f(x, c, n);
+		y = rho_f(rho_f(y, c, n), c, n);
+		d = std::__gcd(y > x ? y-x : x-y, n);
+	}
+	return d;
+}
+void factor(int n, std::vector<int> &res) {
+	if (n == 1) return;
+	else if (isprime(n)) res.push_back(n);
+	else if (n == 4) res.push_back(2), res.push_back(2);
+	else {
+		int d;
+		while ((d = rho(n)) == n);
+		factor(d, res), factor(n / d, res);
+	}
+}
+int ind(int a, int b, int m) {
+	a %= m, b %= m;
+	std::map<int, int> hash;
+	int r = (int)(sqrt(m)), k = 1;
+	if (r * r < m) r++;
+	for (int i = 0; i < r; i++, k = 1LL * k * a % m)
+		if (hash.find(k) == hash.end())
+			hash.insert(std::make_pair(k, i));
+	int s = 1;
+	std::map<int, int>::iterator it;
+	for (int i = 0; i < r; i++, s = 1LL * s * k % m) {
+		int x, y, t;
+		ext_gcd(s, m, x, y);
+		t = 1LL * b * x % m;
+		if ((it = hash.find(t)) != hash.end())
+			return i * r + it->second;
+	}
+}
+void prepare_inv(int *inv, int p) {
+	inv[1] = 1;
+	for (int i = 2; i < p; i++)
+		inv[i] = 1LL * inv[p%i] * (p - p/i) % p;
+}

.ACM-Templates/Math/numbertheory/remark.tex

@@ -0,0 +1,7 @@
+数论模板，其中包含：
+\begin{itemize}
+\item ext\_gcd：扩展欧几里得方法解 $ax+by=\gcd(a,b)$，该函数保证当 $a,b>0$ 时 $x>0$。
+\item flsum：欧几里得思想解 $\sum_{x=st}^{ed}\lfloor\frac{ax+b}c\rfloor$。
+\item ind：小步大步走算法求 $a^x=m\pmod b$。
+\item prepare\_inv：$O(p)$ 求模 $p$ 域下所有非零元的逆元（$p$ 素数）。
+\end{itemize}

.ACM-Templates/Math/numerical/main.cpp

@@ -0,0 +1,43 @@
+#include <cmath>
+#include <algorithm>
+#include <complex>
+const double pi = 3.14159265358979324;
+typedef double (*__F) (double);
+double area(double x, double y, double fl, double fmid, double fr) {
+	return (fl + 4 * fmid + fr) * (y - x) / 6;
+}
+double area_simpson_solve(__F f, double x, double mid, double y, double fl, double fmid, double fr, double pre, double zero) {
+	double lmid = (x + mid) / 2, rmid = (mid + y) / 2;
+	double flmid = f(lmid), frmid = f(rmid);
+	double al = area(x, mid, fl, flmid, fmid), ar = area(mid, y, fmid, frmid, fr);
+	if (fabs(al + ar - pre) < zero) return al + ar;
+	else return area_simpson_solve(f, x, lmid, mid, fl, flmid, fmid, al, zero) + area_simpson_solve(f, mid, rmid, y, fmid, frmid, fr, ar, zero);
+}
+double area_simpson(__F f, double x, double y, double zero = 1e-10) {
+	double mid = (x + y) / 2, fl = f(x), fmid = f(mid), fr = f(y);
+	return area_simpson_solve(f, x, mid, y, fl, fmid, fr, area(x, y, fl, fmid, fr), zero);
+}
+typedef std::complex<double> complex;
+void fft_prepare(int maxn, complex *&e) {
+	e = new complex[2 * maxn - 1];
+	e += maxn - 1;
+	e[0] = complex(1, 0);
+	for (int i = 1; i < maxn; i <<= 1)
+		e[i] = complex(cos(2 * pi * i / maxn), sin(2 * pi * i / maxn));
+	for (int i = 3; i < maxn; i++)
+		if ((i & -i) != i) e[i] = e[i - (i & -i)] * e[i & -i];
+	for (int i = 1; i < maxn; i++) e[-i] = e[maxn - i];
+}
+/* f = 1: dft; f = -1: idft */
+void dft(complex *a, int N, int f, complex *e, int maxn) {
+	int d = maxn / N * f;
+	complex x;
+	for (int n = N, m; m = n / 2, m >= 1; n = m, d *= 2)
+		for (int i = 0; i < m; i++)
+			for (int j = i; j < N; j += n)
+				x = a[j] - a[j+m], a[j] += a[j+m], a[j+m] = x * e[d * i];
+	for (int i = 0, j = 1; j < N - 1; j++) {
+		for (int k = N / 2; k > (i ^= k); k /= 2);
+		if (j < i) std::swap(a[i], a[j]);
+	}
+}

.ACM-Templates/Math/numerical/remark.tex

@@ -0,0 +1,6 @@
+数值计算模板，内含：
+\begin{itemize}
+\item area\_simpson：Simpson法求定积分。
+\item fft\_prepare：求 $maxn$ 次单位根，包括负指数。
+\item dft：做长度为 $N(N|maxn)$ 的DFT，其中若 $f=-1$ 则求IDFT。
+\end{itemize}