dna049 的 C++ 模板

此处只包含通用代码块，以下为其它内容：

C++ 代码规范
C++ 图论模板（因为图论内容实在太多）
C++ 字符串模板（内容相对独立）
STL 和 C++17 学习记录

欢迎使用或转载代码块，唯一要求添加一行注释: https://dna049.com

编译器比你想象中的聪明，所以尽量写 编译器好优化 并且易读的代码，很多位运算，除法都是可以被优化的！
优质的代码本身就是一种解释，添加完全不必要的注释只会让人恶心
以后尽量使用 vector 而非数组，结合 C++17 特征，可以简化代码且便于编译器优化！
以前觉得 main 函数 return 0 只是标准写法，现在(2020/5/22)才知道能返回就能提前优雅的结束，并且可以判断是否正常结束！
全局变量数组元素自动默认初始化为 0，局部变量要加 = {} 才会初始化为 0

此篇代码不再更新，已经发现不少 bug。下面是代码存储位置（并不包含全部此篇内容）：

通用代码块

竞赛入门书籍推荐：GuideToCompetitiveProgramming

佛祖保佑，永无 bug

#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
// #define print(x) std::cout << (x) << std::endl 
// #define clog(x) std::clog << (#x) << " is " << (x) << '\n';
using LL = long long;

int main(){
	//freopen("in","r",stdin);
	std::cin.tie(nullptr)->sync_with_stdio(false);
	auto start = std::clock();
	
	std::cout << "Time used: " << (std::clock() - start) << "ms" << std::endl;
	return 0;
}
/*
------ Welcome to my blog: http://dna049.com ------
                   _ooOoo_
                  o8888888o
                  88" . "88
                  (| -_- |)
                  O\  =  /O
               ____/`---'\____
             .'  \\|     |//  `.
            /  \\|||  :  |||//  \
           /  _||||| -:- |||||-  \
           |   | \\\  -  /// |   |
           | \_|  ''\---/''  |   |
           \  .-\__  `-`  ___/-. /
         ___`. .'  /--.--\  `. . __
      ."" '<  `.___\_<|>_/___.'  >'"".
     | | :  `- \`.;`\ _ /`;.`/ - ` : | |
     \  \ `-.   \_ __\ /__ _/   .-` /  /
======`-.____`-.___\_____/___.-`____.-'======
                   `=---='
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
            佛祖保佑       永无BUG

 * ┌───┐   ┌───┬───┬───┬───┐ ┌───┬───┬───┬───┐ ┌───┬───┬───┬───┐ ┌───┬───┬───┐
 * │Esc│   │ F1│ F2│ F3│ F4│ │ F5│ F6│ F7│ F8│ │ F9│F10│F11│F12│ │P/S│S L│P/B│  ┌┐    ┌┐    ┌┐
 * └───┘   └───┴───┴───┴───┘ └───┴───┴───┴───┘ └───┴───┴───┴───┘ └───┴───┴───┘  └┘    └┘    └┘
 * ┌───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───────┐ ┌───┬───┬───┐ ┌───┬───┬───┬───┐
 * │~ `│! 1│@ 2│# 3│$ 4│% 5│^ 6│& 7│* 8│( 9│) 0│_ -│+ =│ BacSp │ │Ins│Hom│PUp│ │N L│ / │ * │ - │
 * ├───┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─────┤ ├───┼───┼───┤ ├───┼───┼───┼───┤
 * │ Tab │ Q │ W │ E │ R │ T │ Y │ U │ I │ O │ P │{ [│} ]│ | \ │ │Del│End│PDn│ │ 7 │ 8 │ 9 │   │
 * ├─────┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴─────┤ └───┴───┴───┘ ├───┼───┼───┤ + │
 * │ Caps │ A │ S │ D │ F │ G │ H │ J │ K │ L │: ;│" '│ Enter  │               │ 4 │ 5 │ 6 │   │
 * ├──────┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴────────┤     ┌───┐     ├───┼───┼───┼───┤
 * │ Shift  │ Z │ X │ C │ V │ B │ N │ M │< ,│> .│? /│  Shift   │     │ ↑ │     │ 1 │ 2 │ 3 │   │
 * ├─────┬──┴─┬─┴──┬┴───┴───┴───┴───┴───┴──┬┴───┼───┴┬────┬────┤ ┌───┼───┼───┐ ├───┴───┼───┤ E││
 * │ Ctrl│    │Alt │         Space         │ Alt│    │    │Ctrl│ │ ← │ ↓ │ → │ │   0   │ . │←─┘│
 * └─────┴────┴────┴───────────────────────┴────┴────┴────┴────┘ └───┴───┴───┘ └───────┴───┴───┘
 */

递归程序防止爆栈

在 Windows 上，通常的方法是在 编译选项 中加入 -Wl,--stack=1000000000
命令行中可以使用：g++ -static -Wl,--stack=268435456 -O2 -std=c++17 -o main main.cpp
powerShell 中可以使用：g++ -static (Tab 上的那个键)-Wl,--stack=268435456 -O2 -std=c++17 -o main main.cpp
在 Linux 上，通常的方法是在运行程序前 在终端内 执行 ulimit -s unlimited （WSL 下无法设置可惜）

Python 输入样例（以备不时之需，用 PyPy3 提交）

用 Python 过的一次大数题：490C

# 多 case 输入
for _ in range(int(input())):
# 单行输入
    n = int(input())
# 两个元素一行输入
	a, b = map(int, input().split())

# 提前结束	
exit()

产生 log 的几个原因

二分，三分
$1 + \frac{1}{2} + \cdots \frac{1}{n} \sim \log n$
树状数组，线段树
重链剖分
倍增（太好用了，特别是 nxt 数组进行加速）
莫队（离线算法）

产生根号的几个原因

朴素判断素数
整除分块：$\lfloor \frac{n}{i} \rfloor$ 的值域是 $O(\sqrt{n})$ 的
网络流中 HLPP（没读过这篇复杂度分析的论文，不懂）
分块处理，分块打表

类中静态成员定义时初始化

在前面加 inline 即可，最后都加，否则在配合 STL 的时候就会链接出错，从而无法编译通过例如：

#include <bits/stdc++.h>
struct Flow {
	static const int INF = 1e9;
	int n;
	std::vector<int> d;
	Flow(int _n) : n(_n) {}
	int test() {
		d.assign(n, INF);
		return INF;
	}
};

int main() {
	Flow g(3);
	std::cout << g.test() << std::endl;
	return 0;
}

bitset 高端压位卡常

典型应用，求传递闭包。

C++ 黑魔法：$n^2$ 过百万，编译器优化+指令集优化

博文和 ppt

但是实际上不需要这么麻烦，仅需在头部添加下面代码（codeforces 上支持）

1 2	#pragma GCC optimize("Ofast,no-stack-protector,unroll-loops") #pragma GCC target("sse,sse2,sse3,ssse3,sse4.1,sse4.2,abm,mmx,avx,avx2,popcnt,tune=native")

示例：

911G 的提交
1151G 的提交

概率问题永不 TLE 的技巧

// 概率问题永不 TLE 的技巧
auto begin = std::chrono::steady_clock::now();
while ((std::chrono::steady_clock::now() - begin).count() < 5e8) { // 数值取时限的一半
	// do something
}

int VS long long

以前当 int 类型出现带模乘法的时候，我们一般就直接使用 LL 了，不用强制类型转化了。但是后来转而使用 int 了。

注意到 LL(a) * b % n 和 1LL * a * b 在汇编意义下没有区别。两种情况哪个方便写那个。放弃 LL 的主要理由：

本身就是 int 型的变量，为什么用 LL 存储呢？
内存减少一半
在 32 位机器上更快。

const VS constexpr

从 2021-4-27 开始使用 constexpr。对于变量而言，它们的本质区别，一个编译器初始化，一个运行时初始化。

莫队加速

已经学会普通莫队，增加回滚莫队，删除回滚莫队。带修改莫队没兴趣，树上莫队就是 dfs + 莫队。

位运算的关系

异或 1 改变，异或 0 不变
某位异或位 0，表示此位相等，反之不等。
$a \oplus b = (a \mid b) \oplus (a \And b)$
$a \oplus b = (a \mid b) - (a \And b)$
$a + b = (a \mid b) + (a \And b)$
$a + b = (a \oplus b) + 2 (a \And b)$
(a & b) | c = (a | b) & (a | c)
(a | b) & c = (a & b) | (a & c)
(a | b) ^ 1 = (a ^ 1) & (b ^ 1)
(a & b) ^ 1 = (a ^ 1) | (b ^ 1)
(a | b) ^ c 和 (a & b) ^ c 可以逐位转化，因此任何一个数 x 经过任意多次的&, |, ^ 运算最终都可以写成 ((x ^ a) & b) | c。

最大最小值分配律

$\begin{aligned} \min(\max(x, a), b) = \max(\min(x, b), \min(a, b)) \\ \max(\min(x, a), b) = \min(\max(x, b), \max(a, b)) \end{aligned}$

例题：Atcoder abc196E

最高位后面位取反 $O(1)$

1
2
3

int reverseBit(int n) {
	return n ^ ((1 << 32 - __builtin_clz(n)) - 1);
}

最低位 1 置 0: `x = x & (x - 1)`

树状数组中使用的 `lowbit(x) = x & (-x)` 得到 x 的最大 2 的幂次因子。

mask 位上暴力枚举

1
2
3

for (int i = n; i; i = (i - 1) & n) {
	// do something
}

取平均值（防溢出）

1 2	(x & y) + ((x ^ y) >> 1) // 向下取整 (x \| y) - ((x ^ y) >> 1) // 向上取整

Gosper’s Hack：n 个集合中选 k 个

思路：想想怎么把 1011100 变成 110011

void GospersHack(int n, int k) {
	int cur = (1 << k) - 1;
	int limit = (1 << n);
	std::vector<int> r;
	while (cur < limit) {
		// do something
		int lb = cur & -cur;
		int r = cur + lb;
		cur = (((r ^ cur) >> 2) / lb) | r;
	}
}

动态规划

水涨船高技巧

把一个集合中所有元素加一个常数，可以不操作，加在水位线上即可。

Meet in Middle(拆半搜索法)

类似于动态规划，是一种思想。特别适合处理指数复杂度。

例题：AtCoder abc184F，当然针对此题可以深搜剪枝法。

Small to large（把小的合并到大的里面去）

例子：并查集（dus），map 的合并，树上启发式合并（dus on tree)，重链剖分。

例题：600E 的题解

倍增思想（太强了）

例子：RMQ，LCA。nxt 数组进行加速

int128 的使用

__int128 read(){
	__int128 x = 0, f = 1;
	char ch = getchar();
	while (ch < '0' || ch > '9'){
		if (ch == '-')
			f = -1;
		ch = getchar();
	}
	while (ch >= '0' && ch <= '9') {
		x = x * 10 + ch - '0';
		ch = getchar();
	}
	return x * f;
}

void printS(__int128 x){
	if (x > 9) printS(x / 10);
	putchar(x % 10 + '0');
}
void print(__int128 x){ 
	if (x < 0) {
		putchar('-');
		x = -x;
	}
	printS(x);
}

int main(){
	__int128 a = read();
	print(a * a);
	printf("\n");
	return 0;
}

cout, cerr, clog 的区别

cerr 可看作自动刷新缓冲区的 clog。 cout 和 clog 的区别就是 clog 直接打印在屏幕上，而我们文件输入输出的时候使用 freopen("out", "w", stdout) 后 cout 会输出到 out 文件中，而 clog 依然打印在屏幕中，这就很有用了。

交互式题目模板

gym101021: Guess the Number
需要 fflush(stdout);（对于 scanf/printf）或 std:::cout << std::flush （对于 std::cin/std::cout）来刷新缓冲区，不过 std::endl 会自动刷新一次缓冲区，所以此时可以省略。

注意 std::endl 和 \n 的区别是前一个刷新缓冲区，后一个不刷新。仅在交互问题或者 debug 的时候使用 std::endl;

#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
#define print(x) std::cout << (x) << std::endl
using LL = long long;

int main() {
	//freopen("in","r",stdin);
	std::ios::sync_with_stdio(false);
	std::cin.tie(nullptr);
	int l = 1, r = 1e6;
	while (l <= r) {
		int m = (l + r) / 2;
		std::cout << m << std::endl;
		std::string s;
		std::cin >> s;
		if (s[0] == '<') r = m - 1;
		else l = m + 1;
	}
	std::cout << "! " << r << std::endl;
	return 0;
}

交互题测试例子：先假设一组答案，按照自己的程序计算出每次 query 的答案，最好用文件输入输出测试。

负数下标技巧

1
2
3

const int N = 1e5 + 2;
int aa[N];
int *a = (aa + N / 2);

可用于 $O(1)$ 首尾插入或删除元素，访问第 $i$ 个元素。
当然也可以用 std::deque 加一个标号，实现上述操作

优雅的输出技巧

int n = 10;
for(int i = 0; i < n; ++i) {
	std::cout << i << " \n"[i == n - 1];
}

模正数向下取整：$\lfloor \frac{a}{n} \rfloor$ 和模正数向下取整：$\lceil \frac{a}{n} \rceil$

LL floor(LL a, LL n) { // n > 0
	return a < 0 ? (a - n + 1) / n : a / n;
}
LL ceil(LL a, LL n) { // n > 0
	return a < 0 ? a / n : (a + n - 1) / n;
}

注意到 C/C++ 中，整数除法 / 是向 0 取整(int(x)也是向 0 取整)，但是 Python（Sagemath）整数除法 // 是向下取整。在 C++ 中一定不要用 (x - 1) / n + 1 的姿势向上取整！！！

Floor_sum : $\displaystyle \sum_{i = 0}^{n - 1} \lfloor \frac{a \cdot i + b}{m} \rfloor$

注意到这个求和式表示由直线 $x = 0, x = n, y = 0, y = \frac{a}{m}x + \frac{b}{m}$ 构成的梯形内部（仅包含上边界）的整点个数。

显然我们可以先预处理，使得 $a < m$ 且 $b < m$，此时我们考虑 $y_{\max} = \lfloor \frac{a \cdot n + b}{m} \rfloor$，定义 $x_{\max} = y_{\max} \cdot m - b$，那么 $[\lceil \frac{x_{\max}}{a} \rceil, n)$ 的纵坐标取整值都相等为 $y_{\max}$，然后我们横纵坐标互换（在一张相对透明的纸上画图，然后在另一面看就是）由直线 $x = 0, x = y_{\max}, y = 0, y = \frac{m}{a}x + \frac{b’}{a}$ 构成的梯形内部（不包含边界）的整点个数。其中 $b’ = \lceil \frac{x_{\max}}{a} \rceil \cdot a - x_{\max}$（表示直线 $y = y_{\max}$ 与 $y = \frac{a}{m}x + \frac{b}{m}$ 以及 $x = \lceil \frac{x_{\max}}{a} \rceil$ 相交的长度）

LL floorSum(int n, int m, int a, int b) {
	LL r = 0;
	if (a >= m) {
		r += LL(a / m) * (n - 1) * n / 2;
		a %= m;
	}
	if (b >= m) {
		r += LL(b / m) * n;
		b %= m;
	}
	int yMax = (LL(a) * n + b) / m;
	if (yMax == 0) return r;
	LL xMax = LL(yMax) * m - b;
	r += (n - (xMax + a - 1) / a) * yMax;
	r += floorSum(yMax, a, m, (a - (xMax % a)) % a);
	return r;
}

当然我们可以考虑用长方形整点数减去上部分的整点数（要往上平移一个单位）这样搞更快。

LL floorSum(int n, int m, int a, int b) {
	LL r = 0;
	if (a >= m) {
		r += LL(a / m) * (n - 1) * n / 2;
		a %= m;
	}
	if (b >= m) {
		r += LL(b / m) * n;
		b %= m;
	}
	int yMax = (LL(a) * n + b) / m;
	if (yMax == 0) return r;
	r += LL(n - 1) * yMax;
	r -= floorSum(yMax, a, m, m - b - 1);
	return r;
}

模板例题和应用例题

带取整的函数取最值的技巧

先考虑不取整的情况，然后一般这个值是可能的最小值或者最大值
然后通过循环看是否满足取整的情况

Barrent reduction 快速模，弃用，因为并不会变快… ：wiki 镜像解释

对于给定常数 $M$ 求 a % M，并要求 $ 0 \leq a < M^2$，并且 $a < 2^k$。因此下面 $k$ 的取值还是需要注意的。

constexpr LL M  = 1e9 + 7; // too big, M should satisfy M * M < int
constexpr int  k = std::__lg(M) + 2;
constexpr LL m = (1LL << k) / M;

auto mod = [&](int a) {
	LL r = a - ((a * m) >> k) * M;
	return r >= M ? r - M : r;
};

输出全排列

void permutation(int n){
	std::vector<int> x(n);
	std::iota(x.begin(), x.end(), 1);
	do {
		std::for_each(x.begin(), x.end(),[](int i){
			std::cout << i;
		});
		std::cout << std::endl;
	} while (std::next_permutation(x.begin(), x.end()));
}

输出全排列的原理

首先初始状态从小到大排列，然后对每一个状态考虑它的后缀，如果后缀是从大到小排列，再考虑向前一位的后缀，直到不是从大到小排列，然后找比第一个位置大的最小值放在开头，其它位置排序。

import math
def permutation(n):
	ans = []
	cnt = math.factorial(n);
	r = list(range(1, n + 1))
	ans.append(r.copy())
	cnt -= 1
	while cnt > 0:
		i = n - 1
		while r[i - 1] > r[i]: i -= 1
		r[i:] = r[i:][::-1]
		j = i
		while r[j] < r[i - 1]: j += 1
		r[i - 1], r[j] = r[j], r[i - 1]
		ans.append(r.copy())
		cnt -= 1
	return ans	

for i in range(1,5):
	print(permutation(i))

初等数论

Greatest Common divisor

// 简洁写法，不推荐，推荐使用内建 __gcd
LL gcd(LL a, LL b){
	return b ? gcd(b, a % b) : a;
}
// lambda 表达式写法，开头不能是auto因为递归
std::function<int(int, int)> gcd = [&](int a, int b)->int{
	return b == 0 ? a : gcd(b, a % b);
};
// 快速版本 https://cp-algorithms.com/algebra/euclid-algorithm.html
LL gcd(LL a, LL b) {
	if (!a || !b) return a | b;
	unsigned shift = __builtin_ctzll(a | b);
	a >>= __builtin_ctzll(a);
	do {
		b >>= __builtin_ctzll(b);
		if (a > b) std::swap(a, b);
		b -= a;
	} while (b);
	return a << shift;
}
// 普通版拓展GCD
LL exGcd(LL a, LL b, LL& x, LL& y){  // ax + by = gcd(a,b)
	if (b == 0){
		x = 1; y = 0; return a;
	}
	LL d = exgcd(b, a % b, y, x);
	y -= a / b * x;
	return d;
}
// C++17版拓展GCD，优雅了不少!
std::tuple<LL, LL, LL> exGcd(LL a, LL b) {  // ax + by = gcd(a,b)
	if (b == 0) return {a, 1, 0};
	auto [d, y, x] = exGcd(b, a % b);
	return {d, x, y - a / b * x};
}

Sum of least common multiple $s_n = \sum_{i=1} ^n lcm(i,n)$

$\begin{aligned} \sum_{i=1} ^n lcm(i,n) &= \sum_{i=1} ^n \frac{i \cdot n}{\gcd(i,n)} \\ &= n \sum_{d|n} \sum_{\gcd(i,\frac{n}{d})= 1} i \\ &= \frac{n}{2} \sum_{d|n} \frac{n}{d} \psi(\frac{n}{d}) \\ &= \frac{n}{2} \sum_{d|n} d \psi(d) \end{aligned}$

所以，我们可以在 $O(n \log n)$ 处理好 $s_n$ 的前 $n$ 项。

Double sum of least common multiple $ds_n = \sum_{1 \leq i \leq j \leq n} lcm(i,j)$

本来这个也挺麻烦，但是可以借助 $s_n$ 计算：$ds_n = \sum_{j=1} ^n s_j$，所以复杂度就一致了。当然也可以直接化简成：

$\sum_{1 \leq i \leq j \leq n} lcm(i,j) = \sum_{i = 1} ^ n \frac{i}{4} \psi(i) \lfloor \frac{n}{i} \rfloor (\lfloor \frac{n}{i} \rfloor + 1)$

不借助 $s_n$ 其实也能暴力搞出来的。

模乘法逆元

1
2
3

LL inv(LL a, LL p){ // 0 < a < p and gcd(a,p) = 1
	return a == 1 ? 1 : (p - p / a) * inv(p % a, p) % p;
}

上述代码主要用于线性时间预处理所有$p$以内的逆元，对于较小的常数$a$, 可以直接试除 b p mod a == 1
用下面快速幂也可以求逆，inv 的步数平均下面显著的比快速幂小，但是由于用到了递归，因此最终它们的平均效率是一致的。
可以通过预处理小部分值达到快速的效果。

快速模乘法

// 0 <= x < p < INT_MAX
LL powMod(LL x, LL n, LL p){
	LL r = 1;
	while (n) {
		if (n&1) r = r * x % p;
		n >>= 1; x = x * x % p;
	}
	return r;
}

可以用 for 循环写的更短一点，但没必要，若 p < INT64_MAX 就需要使用 __int128 了。现代计算机都是 64 位的，因此处理 int 跟 long long 基本没有时间差异，直接用 LL 避免了强制类型转化，因此效率反而更高。但是 __int128 就不同了。

// 0 <= x < p < INT64_MAX
LL powMod(LL x, LL n, LL p){
	LL r = 1;
	while (n) {
		if (n&1) r = __int128(r) * x % p;
		n >>= 1; x = __int128(x) * x % p;
	}
	return r;
}

阶乘、组合数、Lucas 定理

const int M = 1e9 + 7;
const int N = 3e5 + 2;
LL fac[N], ifac[N];
LL inv(LL a){
	return a == 1 ? 1 : (M - M / a) * inv(M % a) % M;
}
void init() {
	fac[0] = 1;
	for (int i = 1; i < N; ++i) fac[i] = fac[i - 1] * i % M;
	ifac[N - 1] = inv(fac[N - 1]);
	for (int i = N - 1; i; --i) ifac[i - 1] = ifac[i] * i % M;
}
LL binom(int n, int k) {
	//if (n < k || n < 0) return 0;
	return fac[n] * ifac[k] % M * ifac[n - k] % M;
}
LL lucas(LL n, LL m, LL p) { // C(n,m)%p, 仅在p较少时发挥作用
	LL r = 1;
	while (n && m) {
		LL np = n % p, mp = m % p;
		if (np < mp) return 0;
		r = binom(np, mp);
		n /= p, m /= p;
	}
	return r;
}

常用组合数公式

对任意实数，定义：$\binom{\alpha}{k} = \frac{\alpha(\alpha - 1) \cdots (\alpha - k + 1)}{k !}$ 所以我们有：

$\binom{-n}{k} = (-1)^{k} \binom{n + k - 1}{k}$ $\binom{n + 1}{k + 1} = \binom{n}{k} + \binom{n}{k + 1}$

最后一个数，先还是不选，这是一个问题。

${n \choose k}{k \choose i} = {n \choose i} {n - i \choose k - i}$

组合意义理解：$n$ 个人中选出 $i$ 个一流人才， $k - i$ 个二流人才。

${n + m \choose k} = \sum_{i + j = k} {n \choose i} {m \choose j}$

组合意义理解：$n, m$ 两个堆选出 $k$ 个人。

乞丐版素数判断

bool isPrime(int n) {
	if (n == 2) return true;
	if (n == 1 || n % 2 == 0) return false;
	for (int i = 3; i * i <= n; i += 2){
		if (n % i == 0) return false;
	}
	return true;
}

$O(n \log n)$ 素数筛

const int N =  1e7+2;
bool isP[N];

void initPrime(){
	if(isP[2]) return;
	isP[2] = true;
	for (int i = 3; i < N; i += 2) isP[i] = true;
	for (int i = 3; i < N; i += 2) if (isP[i]) {
		for (int j = 3 * i; j < N; j += i * 2) {
			isP[j] = false;
		}
	}
}

欧拉线性素数筛（正式使用版）

const int N = 1e9 + 9;
bool isp[N];
std::vector<int> p;
int initPrimeP() {
	p.emplace_back(2);
	isp[2] = true;
	for (int i = 3; i < N; i += 2) isp[i] = true;
	int sq = int(std::sqrt(N + 0.1))|1;
	for (int i = 3; i <= sq; i += 2) if (isp[i]) {
		p.emplace_back(i);
		for (int j = i * i; j < N; j += i << 1) {
			isp[j] = false;
		}
	}
	for (int i = sq + 2; i < N; i += 2) if (isp[i]) p.emplace_back(i);
	return p.size();
}

欧拉线性素数筛（弃用）

const int N = 1e9 + 9;
bool isp[N];
std::vector<int> p;

void initPrimeP() {
	p.emplace_back(2);
	isp[2] = true;
	for (int i = 3; i < N; i += 2) isp[i] = true;
	for (int i = 3; i < N; i += 2) {
		if (isp[i]) p.emplace_back(i);
		for (int j = 1, t = (N - 1) / i + 1; j != p.size() && p[j] < t; ++j) { // 用除号是防止溢出
			isp[i * p[j]] = false;
			// 不要下面的一步的话，复杂度 O(nloglogn), 但是不用除法，常数小
			if (i % p[j] == 0) break;
		}
	}
}

$O(1)$ （执行）时间复杂度判断一个数是否为素数

奇技淫巧来源：https://codeforces.com/blog/entry/79941#comment-659202

template<int N>
struct Sieve {
	bool isP[N];
	constexpr Sieve(): isP() {
		isP[2] = true;
		for (int i = 3; i < N; i += 2) isP[i] = true;
		for (int i = 3; i < N; i += 2) if (isP[i]) {
				for (int j = 3 * i; j < N; j += i * 2) {
					isP[j] = false;
				}
			}
	}
};
// MAXN 默认最大值为1<<18=262144, 调节参数 -fconstexpr-loop-limit= 例如：
// g++ main.cpp -std=c++17 -fconstexpr-loop-limit=12345678 -fconstexpr-ops-limit=1234567890
// 使得 MAXN = 1e7+2
constexpr int MAXN = 1e5+2;
bool fast_is_prime(int n) {
	static constexpr Sieve<MAXN> s;
	return s.isP[n];
}

大素数 Miller-Rabin 概率判别法和大整数的最大/小因子分解

namespace PollardRho {
std::mt19937 rnd(std::chrono::steady_clock::now().time_since_epoch().count());
LL powMod(LL x, LL n, LL p) {
	LL r = 1;
	while (n) {
		if (n & 1) r = __int128(r) * x % p;
		n >>= 1;   x = __int128(x) * x % p;
	}
	return r;
}
// 1 < a < n，若 n 是素数，那么 a^(n - 1) = 1 mod n
// m - 1 = m * 2 ^ t，返回 false 表示判断失败
bool witness(LL a, LL n, LL m, int t) {
	LL x = powMod(a, m, n);
	if (x == 1 || x == n - 1) return false;
	while (t--) {
		x = __int128(x) * x %  n;
		if (x == n - 1) return false;
	}
	return true;
}
const int TIMES = 52;
bool rabin(LL n) {
	if (n < 2) return false;
	if (n == 2) return true;
	if (n % 2 == 0) return false;
	LL m = n - 1;
	int t = __builtin_ctzll(m);
	m >>= t;
	for (int cnt = 0; cnt < TIMES; ++cnt) {
		LL a = rnd() % (n - 1) + 1;
		if (witness(a, n, m, t)) return false;
	}
	return true;
}
LL pollardrho(LL n) {
	LL x = 0, y = 0, z = 1, i = 1, k = 2, c = rnd() % (n - 1) + 1;
	while (true) {
		x = (__int128(x) * x + c) % n;
		z = __int128(y - x + n) * z % n;
		// 累计 gcd 一次计算！太猛了啊 茶茶白
		if (++i == k) {
			LL d = std::__gcd(z, n);
			if (d > 1) return d;
			y = x;
			if (k > n) return n;
			k <<= 1;
		}
	}
}
LL spf(LL n) {
	if (rabin(n) || n == 1) return n;
	LL d = n;
	while (d == n) d = pollardrho(n);
	return std::min(spf(d), spf(n / d));
}
LL gpf(LL n, LL mxf = 1) {
	if (rabin(n)) return n;
	if (n <= mxf) return 1;
	LL d = n;
	while (d == n) d = pollardrho(n);
	LL res = gpf(d, mxf);
	return std::max(res, gpf(n / d, std::max(res, mxf)));
}
};

$O(n)$ 最小素因子预处理

const int N = 2e8; // 再大内存吃不消了 
int sp[N], p[N];
int spf() {
	int cnt = 1;
	p[1] = 2;
	for (int i = 2; i < N; i += 2) sp[i] = 2;
	for (int i = 1; i < N; i += 2) sp[i] = i;
	for (int i = 3; i < N; i += 2) {
		if (sp[i] == i) p[++cnt] = i;
		for (int j = 1; j <= cnt && p[j] <= sp[i] && i * p[j] < N; ++j) { //防止乘法溢出
			sp[i * p[j]] = p[j]; // 注意到sp只被赋值一次
		}
	}
	return cnt;
}

Möbius function

另类递推公式： $ \mu(i) = - \sum_{d \mid i, d < i} \mu(d) $。

// 乞丐版
int getMu(int n){
	if (n % 4 == 0) return 0;
	int r = 1;
	if (n % 2 == 0 ){
		n /= 2;
		r = -r;
	}
	for (int i = 3; i * i <= n; i += 2){
		if (n % i == 0) {
			n /= i;
			if(n % i  == 0) return 0;
			r = -r;
		}
	}
	return n > 1 ? -r : r;
}
// O(n log n) 预处理版本
const int N  = 1e6 +2;
int mu[N];
void initmu() {
	mu[1] = 1;
	for (int i = 1; i < N; ++i) {
		for (int j = i * 2; j < N; j += i) {
			mu[j] -= mu[i];
		}
	}
}

// O(n) 预处理版
const int N  = 1e6 + 2;
bool isp[N];
int p[N], mu[N];
int initmu() {
	if(mu[1] == 1) return;
	int cnt = 1;
	p[1] = 2;
	mu[1] = 1;
	isp[2] = true;
	for (int i = 3; i < N; i += 2) isp[i] = true;
	for (int i = 3; i < N; i += 2) {
		if (isp[i]) mu[i] = -1, p[++cnt] = i;
		for (int j = 1, t = (N - 1) / i + 1; j <= cnt && p[j] < t; ++j) {
			isp[i * p[j]] = false;
			if (i % p[j] == 0) break;
			mu[i * p[j]] = -mu[i];
		}
	}
	for (int i = 2; i < N; i += 4) mu[i] = -mu[i >> 1];
	return cnt;
}

Möbius function 前缀和

int sumMu[N];
void initSumMu() {
	if (mu[1] != 1) initmu();
	for (int i = 1; i < N; ++i) sumMu[i] = sumMu[i - 1] + mu[i];
}
std::map<int, int> mp;
int getSumMu(int n) { // M(n) = M(n-1) + mu(n)
	if (n < N) return sumMu[n];
	auto it = mp.find(n);
	if (it != mp.end()) return it->second;
	int r = 1;
	for (int i = 2, j; i <= n; i = j + 1) {
		j = n / (n / i);
		r -= (j - i + 1) * sumMu(n / i);
	}
	mp.insert(std::make_pair(n, r));
	return r;
}
//  Möbius function 绝对值前缀和
LL getSumAbsMu(LL n) { // Q(n) = Q(n-1) + |mu(n)|
	LL r = 0;
	for (LL i = 1, t; (t = i * i) < n; ++i) {
		r += mu[i] * (n / t);
	}
	return r;
}

Euler’s totient function

const int N = 1e7+2;
int phi[N];
void initPhi() {
	if (phi[2] == 1) return;
	for (int i = 1; i < N; i += 2) phi[i] = i;
	for (int i = 2; i < N; i += 2) phi[i] = i >> 1;
	for (int i = 3; i < N; i += 2) {
		if (phi[i] != i) continue;
		for (int j = i; j < N; j += i) phi[j] = phi[j] / i * (i - 1);
	}
}
int getPhi(int n) {
	int r = (n % 2 == 0 ? n/2 : n);
	while (n % 2 == 0) n/=2;
	for (int i = 3; i * i <= n; i += 2) {
		if (n % i  == 0) {
			r = r / i *(i-1);
			while (n % i == 0) n/=i;
		}
	}
	if (n > 1) r = r / n * (n - 1);
	return r;
}
LL sumPhi[N];
void initSumphi() {
	if (phi[2] != 1) initPhi();
	for (int i = 1; i < N; ++i) sumPhi[i] = sumPhi[i - 1] + phi[i];
}
std::map<int, LL> mp;
LL getSumphi(int n) {
	if (n < N) return (LL) sumPhi[n];
	auto it = mp.find(n);
	if (it != mp.end()) return it->second;
	LL r = LL(n + 1) * n / 2;
	for (int i = 2, j; i <= n; i = j + 1) {
		j = n / (n / i);
		r -= (j - i + 1) * getSumphi(n / i);
	}
	mp.insert(std::make_pair(n, r));
	return r;
}

拓展 Euler 定理

数论中欧拉定义说：若 $\gcd(a, m) = 1$ 则 $a^{\phi(m)} \equiv 1 \mod m$。

类似于拓展的 Fermat 小定理：$a^p \equiv a \mod p$，我们有拓展 Euler 定理：

$a^n \equiv a^{n \mod \phi(m) + phi(m)} \mod m$

证明对 m 素因子分解，再利用 Euler 函数是可乘函数，显然。

$\pi(x)$ 函数

const int M = 7;
const int PM = 2 * 3 * 5 * 7 * 11 * 13 * 17;
int phi[PM + 1][M + 1], sz[M + 1];
void init() {
	initprime();
	pi[2] = 1;
	for (int i = 3; i < N; ++i) {
		pi[i] = pi[i - 1];
		if (isp[i]) ++pi[i];
	}
	sz[0] = 1;
	for (int i = 0; i <= PM; ++i) phi[i][0] = i;
	for (int i = 1; i <= M; ++i) {
		sz[i] = p[i] * sz[i - 1];
		for (int j = 1; j <= PM; ++j) {
			phi[j][i] = phi[j][i - 1] - phi[j / p[i]][i - 1];
		}
	}
}
LL primepi(LL x);
LL primephi(LL x, int s) {
	if (s <= M) return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s];
	if (x / p[s] <= p[s]) return primepi(x) - s + 1;
	if (x / p[s] / p[s] <= p[s] && x < N) {
		int s2x = pi[(int)(sqrt(x + 0.2))];
		LL ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;
		for (int i = s + 1; i <= s2x; ++i) {
			ans += pi[x / p[i]];
		}
		return ans;
	}
	return primephi(x, s - 1) - primephi(x / p[s], s - 1);
}
LL primepi(LL x) {
	if (x < N) return pi[x];
	int ps2x = primepi(int(sqrt(x + 0.2)));
	int ps3x = primepi(int(cbrt(x + 0.2)));
	LL ans = primephi(x, ps3x) + LL(ps2x + ps3x - 2) * (ps2x - ps3x + 1) / 2;
	for (int i = ps3x + 1, ed = ps2x; i <= ed; ++i) {
		ans -= primepi(x / p[i]);
	}
	return ans;
}

$\pi(x)$ 函数计算的另一种做法

const int N = 1e7 + 2;
LL L[N], R[N];

LL primepi(LL n) {
	LL rn = (LL) sqrt(n + 0.2);
	for (LL i = 1; i <= rn; ++i) R[i] = n / i - 1;
	LL ln = n / (rn + 1) + 1;
	for (LL i = 1; i <= ln; ++i) L[i] = i - 1;
	for (LL p = 2; p <= rn; ++p) {
		if (L[p] == L[p - 1]) continue;
		for (LL i = 1, tn = std::min(n / (p * p), rn); i <= tn; ++i) {
			R[i] -= (i * p <= rn ? R[i * p] : L[n / (i * p)]) - L[p - 1];
		}
		for (LL i = ln; i >= p * p; --i) {
			L[i] -= L[i / p] - L[p - 1];
		}
	}
	return R[1];
}

求原根

首先，模 $m$ 有原根的充要条件是：$m=2,4,p^a,2p^a$，其中 $p$ 为奇素数。

对于求模奇素数 $p$ 的原根方法：对 $p-1$ 素因子分解：$p-1 = p_1^{a_1} \cdots p_s^{a_s}$ 若恒有

$g^{\frac{p-1}{p_i}} \neq 1(\mod \; p)$

则 $g$ 是模 $p$ 的原根。对于 $p^a$ 的原根为 $g$ 或 $g + p$，若 $p^a$ 的原根为 $g_a$ 而 $2p^a$ 的原根为 $g_a$ 和 $g_a + p^a$ 的中奇数者。所有原根：可以先求出一个原根，然后所有的数便是 $g^0, g^1, \cdots, g^{\phi(m) - 1}$, 所有原根就是那些 $\gcd(i, \phi(m)) = 1$ 的 $g^i$ (证明见 P150《数论基础》潘承洞)。

基于上述方法的代码实现

另外更好的做法：

我们首先求出 $m = \phi(n)$，然后一个个的搜，搜索的时候附带把很多点都给剔除了，所以很快就能找到！具体代码

数论函数的 Dirichlet 乘积

以前的代码不想贴了，不优雅，下次有题做的时候补上。

中国剩余定理

pair<LL,LL> crt2(LL a1,LL m1,LL a2,LL m2){ // x = ai mod mi, m_i >0
	LL t1,t2,ans = a2-a1;
	LL d = exgcd(m1,m2,t1,t2);
	assert(ans%d == 0);
	LL m = m1/d*m2;
	ans = (a1+ans/d*t1%m2*m1)%m; // %m2 是避免溢出
	return make_pair(ans>0?ans:ans+m,m);
}
const int N = 22;
LL a[N],m[N];
pair<LL,LL> crt(int n){ // x = a[i] mod m[i], m[i] >0
	pair<LL,LL> ans = make_pair(a[0]%m[0],m[0]);
	for(int i=1;i<n;++i){
		ans = crt2(ans.first,ans.second,a[i],m[i]);
	}
	return ans;
}

离散对数

LL baby_step_giant_step(LL a, LL b, LL p) { // a^x = b mod p
	a %= p, b %= p;
	if (a == 0) return b % p ? -1 : 1;
	if (b == 1) return 0;
	LL cnt = 0, t = 1;
	for (LL g = std::__gcd(a, p); g != 1; g = std::__gcd(a, p)) {
		if (b % g) return -1;
		p /= g, b /= g, t = t * (a / g) % p;
		++cnt;
		if (b == t) return cnt;
	}
	std::map<LL, LL> mp;
	LL m = LL(std::sqrt(p + 0.1) + 1);
	LL base = b;
	for (LL i = 0; i != m; ++i) {
		mp[base] = i;
		base = base * a % p;
	}
	base = powMod(a, m, p);
	for (LL i = 1; i <= m + 1; ++i) {
		t = t * base % p;
		if (mp.count(t)) return i * m - mp[t] + cnt;
	}
	return -1;
}

模素数开根号

LL modsqrt(LL a, LL p) { // find x s.t x*x=a mod p;
	a = (p + a % p) % p;
	if (a == 0) return 0;
	if (p == 2) return (a & 1) ? 1 : 0;
	LL q = (p - 1) >> 1;
	if (powMod(a, q, p) != 1) return -1;
	if (q & 1) return powMod(a, (q + 1) >> 1, p);
	LL b, cnt = 1;
	while (powMod(b = rand() % p, q, p) == 1); //find a non quadratic residue
	while (!(q & 1)) ++cnt, q >>= 1;
	b = powMod(b, q, p);
	LL x = powMod(a, (q + 1) >> 1, p);
	for (LL s = 1, t = powMod(a, q, p); t != 1; s = 1) { //keep x*x=t*a;t^{2^{cnt-1}}=1
		for (LL tt = t * t % p; s < cnt && tt != 1; ++s) tt = tt * tt % p;
		LL d = powMod(b, 1 << (cnt - s - 1), p);
		x = (x * d) % p;
		b = d * d % p;
		t = t * b % p;
		cnt = s;
	}
	return x;
}
LL mod_sqrt(LL a, LL p, int k = 1) { //find smallest x>=0 s.t x*x=a mod p^k
	if (a == 0) return 0;
	int ka = 0;
	while (a % p == 0) a /= p, ++ka, --k;
	if (k <= 0) return 0;
	if (ka & 1) return -1;
	auto pow = [](LL x, int n) {
		LL r=1;
		for (; n; n >>= 1, x *= x)  if (n&1) r *= x;
		return r;
	};
	LL n = pow(p, k), x;
	if (p == 2) {
		if (k == 1 || k == 2) x = a == 1 ? 1 : -1;
		else if (a % 8 != 1) x = -1;
		else {
			x = 1;
			for (int i = 4; i <= k; ++i) {
				if ((x * x) % (1 << i) == a % (1 << i)) continue;
				x += 1 << (i - 2);
			}
		}
	} else x = mod_sqrt_p(a, n, p, k);
	return x == -1 ? -1 : pow(p, ka >> 1) * (x < n - x ? x : n - x);
}

对于模一般的 $n$，先素因子分解分别求出答案，然后用中国剩余定理求最终解。

自然数方幂和 $O(k)$ 算法

const int N = 1e6 + 6;
int sp[N], p[N];
LL inv[N], AP[N], AS[N], f[N];
LL getPowSum(LL n, int k, LL mod) { // mod > k
	if (k == 0) return n % mod;
	if (p[0] != 2) spf();
	int nk = k + 1;
	LL tmp = 1;
	for (int i = 2; i <= nk; ++i) tmp = tmp * i % mod;
	inv[nk] = powMod(tmp, mod - 2, mod);
	for (int i = nk - 1; i >= 0; --i) inv[i] = inv[i + 1] * (i + 1) % mod;
	f[0] = 0;
	f[1] = 1;
	for (int i = 2; i <= nk; ++i) {
		if (sp[i] == i) f[i] = powMod(i, k, mod);
		else f[i] = f[sp[i]] * f[i / sp[i]] % mod;
	}
	for (int i = 1; i <= nk; ++i) {
		f[i] += f[i - 1];
		if (f[i] >= mod) f[i] -= mod;
	}
	if (n <= nk) return f[n];
	AP[0] = AS[nk] = 1;
	for (int i = 1; i <= nk; ++i) AP[i] = AP[i - 1] * (n + 1 - i) % mod;
	for (int i = nk - 1; i >= 0; --i) AS[i] = AS[i + 1] * (n - i - 1) % mod;
	LL res = 0;
	for (int i = 1; i <= nk; ++i) { // because f(i)=0
		LL x = f[i] * AP[i] % mod * AS[i] % mod * inv[i] % mod * inv[nk - i] % mod;
		if ((nk - i) & 1) res -= x; // be careful
		else res += x;
	}
	return (res % mod + mod) % mod;
}

自然数方幂和精确版

#include <boost/multiprecision/cpp_int.hpp>
using BINT = boost::multiprecision::cpp_int;

BINT f[N];
BINT getPowSum(LL n, int k) { // k<1000
	if (k == 0) return BINT(n);
	if (p[1] != 2) spf();
	int nk = 2 * k + 1;
	f[0] = 0;
	f[1] = 1;
	auto bPow = [](BINT x, int n) -> BINT {
		BINT r(1);
		for (; n; x *= x, n >>= 1) if (n&1) r *= x;
		return r;
	};
	for (int i = 2; i <= nk + 1; ++i) {
		if (sp[i] == i) f[i] = bPow(BINT(i), k);
		else f[i] = f[sp[i]] * f[i / sp[i]];
	}
	for (int i = 1; i <= nk; ++i) f[i] += f[i - 1];
	if (n <= nk) return f[n];
	BINT res = 0, tl = 1, tr = 1;
	for (int i = nk - 1; i >= 0; --i) tr = tr * (n - i - 1) / (nk - i);
	for (int i = 0; i <= nk; ++i) {
		if ((nk - i) & 1) res -= f[i] * tl * tr;
		else res += f[i] * tl * tr;
		tl = tl * (n - i) / (i + 1);
		tr = tr * (nk - i) / (n - i - 1);
	}
	return res;
}

需要下载boost 包类似的包还有 NTL，GMP

生成函数

有位伟人说过：生成函数就像袋子，把一个个零碎的小部件放在一个袋子里，就可以优雅高效的只关注袋子了。

在 codeforces 上 zscoder 大佬给了一个入门教程和进阶教程还有 MiFaFaOvO 的终极教程

生成函数分两种：Original generating function，Expentional generating function，选择哪一种是看问题中是否牵扯组合数。无论哪一种都能保存原数列的全部信息，并且由于级数可以使用微积分和常微分方程的技术，所以会变得更好处理。然后大概率可以优化算法复杂度 $O(n^2) \to O(n \log n)$

关于生成函数多项式的处理：https://cp-algorithms.com/algebra/polynomial.html

多项式高效运算模板：https://github.com/e-maxx-eng/e-maxx-eng-aux/blob/master/src/polynomial.cpp

生成函数一般的处理思路：计算生成函数，分解成有分母不超过二次的分式之和，然后每一个二次的分母部分找一个递推数列来搞定。

OI-wiki 多项式运算

多项式计数杂谈很值得学习一下。

多项式取对数和指数

$B(z) = e^{A(z)}$，即 $A(z) = \ln B(z)$ (不妨假设 $A(0) = 0$ 或等价地 $B(0) = 1$)

那么 $B’(z) = A’(z) \cdot B(z)$, 所以 $[z^{n - 1}] B’(z) = \sum_{k = 0}^{n - 1} [z^k] A’(z) \cdot B(z) [z^{n - 1 - k}] = \sum_{k = 1}^{n} [z^{k - 1}] A’(z) \cdot B(z) [z^{n-k}]$，从而

$n [z^n] B(z) = \sum_{k = 1}^n k [z^k] A(z) \cdot B(z) [z^{n - k}]$

上式等价于

$n [z^n] A(z) = n [z^n] B(z) - \sum_{k = 1}^{n - 1} k [z^k] A(z) \cdot B(z) [z^{n - k}]$

参考：Soulist

NTT 正式可用版(last updated: 2020/7/9)

const LL M = 998244353, ROOT = 3;
LL powMod(LL x, LL n) {
	LL r(1);
	while (n) {
		if (n & 1) r = r * x % M;
		n >>= 1;
		x = x * x % M;
	}
	return r;
}
void bitreverse(std::vector<LL> &a) {
	for (int i = 0, j = 0; i != a.size(); ++i) {
		if (i > j) std::swap(a[i], a[j]);
		for (int l = a.size() >> 1;
			(j ^= l) < l; l >>= 1);
	}
}
void ntt(std::vector<LL> &a, bool isInverse = false) {
	LL g = powMod(ROOT, (M - 1) / a.size());
	if (isInverse) {
		g = powMod(g, M - 2);
		LL invLen = powMod(LL(a.size()), M - 2);
		for (auto & x: a) x = x * invLen % M;
	}
	bitreverse(a);
	std::vector<LL> w(a.size(), 1);
	for (int i = 1; i != w.size(); ++i) w[i] = w[i - 1] * g % M;
	auto addMod = [](LL x, LL y) {
		return (x += y) >= M ? x -= M : x;
	};
	for (int step = 2, half = 1; half != a.size(); step <<= 1, half <<= 1) {
		for (int i = 0, wstep = a.size() / step; i != a.size(); i += step) {
			for (int j = i; j != i + half; ++j) {
				LL t = (a[j + half] * w[wstep * (j - i)]) % M;
				a[j + half] = addMod(a[j], M - t);
				a[j] = addMod(a[j], t);
			}
		}
	}
}
void mul(std::vector<LL>& a, std::vector<LL> b) {
	int sz = 1, tot = a.size() + b.size() - 1;
	while (sz < tot) sz *= 2;
	a.resize(sz);
	b.resize(sz);
	ntt(a);
	ntt(b);
	for (int i = 0; i != sz; ++i) a[i] = a[i] * b[i] % M;
	ntt(a, 1);
	a.resize(tot);
}

siyuan 的 FWT 模板（弃用）

#include <cstdio>
#include <algorithm>
#include <vector>

const int P = 998244353;

void add(int &x, int y) {
    (x += y) >= P && (x -= P);
}
void sub(int &x, int y) {
    (x -= y) < 0 && (x += P);
}
struct FWT {
    int extend(int n) {
        int N = 1;
        for (; N < n; N <<= 1);
        return N;
    }
    void FWTor(std::vector<int> &a, bool rev) {
        int n = a.size();
        for (int l = 2, m = 1; l <= n; l <<= 1, m <<= 1) {
            for (int j = 0; j < n; j += l) for (int i = 0; i < m; i++) {
                if (!rev) add(a[i + j + m], a[i + j]);
                else sub(a[i + j + m], a[i + j]);
            }
        }
    }
    void FWTand(std::vector<int> &a, bool rev) {
        int n = a.size();
        for (int l = 2, m = 1; l <= n; l <<= 1, m <<= 1) {
            for (int j = 0; j < n; j += l) for (int i = 0; i < m; i++) {
                if (!rev) add(a[i + j], a[i + j + m]);
                else sub(a[i + j], a[i + j + m]);
            }
        }
    }
    void FWTxor(std::vector<int> &a, bool rev) {
        int n = a.size(), inv2 = (P + 1) >> 1;
        for (int l = 2, m = 1; l <= n; l <<= 1, m <<= 1) {
            for (int j = 0; j < n; j += l) for (int i = 0; i < m; i++) {
                int x = a[i + j], y = a[i + j + m];
                if (!rev) {
                    a[i + j] = (x + y) % P;
                    a[i + j + m] = (x - y + P) % P;
                } else {
                    a[i + j] = 1LL * (x + y) * inv2 % P;
                    a[i + j + m] = 1LL * (x - y + P) * inv2 % P;
                }
            }
        }
    }
    std::vector<int> Or(std::vector<int> a1, std::vector<int> a2) {
        int n = std::max(a1.size(), a2.size()), N = extend(n);
        a1.resize(N), FWTor(a1, false);
        a2.resize(N), FWTor(a2, false);
        std::vector<int> A(N);
        for (int i = 0; i < N; i++) A[i] = 1LL * a1[i] * a2[i] % P;
        FWTor(A, true);
        return A;
    }
    std::vector<int> And(std::vector<int> a1, std::vector<int> a2) {
        int n = std::max(a1.size(), a2.size()), N = extend(n);
        a1.resize(N), FWTand(a1, false);
        a2.resize(N), FWTand(a2, false);
        std::vector<int> A(N);
        for (int i = 0; i < N; i++) A[i] = 1LL * a1[i] * a2[i] % P;
        FWTand(A, true);
        return A;
    }
    std::vector<int> Xor(std::vector<int> a1, std::vector<int> a2) {
        int n = std::max(a1.size(), a2.size()), N = extend(n);
        a1.resize(N), FWTxor(a1, false);
        a2.resize(N), FWTxor(a2, false);
        std::vector<int> A(N);
        for (int i = 0; i < N; i++) A[i] = 1LL * a1[i] * a2[i] % P;
        FWTxor(A, true);
        return A;
    }
} fwt;

int main() {
    int n;
    scanf("%d", &n);
    std::vector<int> a1(n), a2(n);
    for (int i = 0; i < n; i++) scanf("%d", &a1[i]);
    for (int i = 0; i < n; i++) scanf("%d", &a2[i]);
    std::vector<int> A;
    A = fwt.Or(a1, a2);
    for (int i = 0; i < n; i++) {
        printf("%d%c", A[i], " \n"[i == n - 1]);
    }
    A = fwt.And(a1, a2);
    for (int i = 0; i < n; i++) {
        printf("%d%c", A[i], " \n"[i == n - 1]);
    }
    A = fwt.Xor(a1, a2);
    for (int i = 0; i < n; i++) {
        printf("%d%c", A[i], " \n"[i == n - 1]);
    }
    return 0;
}

FWT 和 FMT 本质上是一致的，只是写法姿势不同。FWT 是用递归的思想写的，FMT 是用逐位计算写的。别人的总归用起来不舒服，还是用下面自己写的 FMT 模板吧！

FMT 模板

考虑数论函数的 Dirichlet 积中的 Möbius 变换（本质就是带条件的求和）这里的变换也就是的各个位上的带条件求和。
例题：洛谷 P6097，参考：yijan

namespace FMT {
const int M = 998244353;
const int inv2 = (M + 1) / 2;
auto add = [](int &x, int y) {
	(x += y) >= M && (x -= M);
};
auto sub = [](int &x, int y) {
	(x -= y) < 0 && (x += M);
};
auto extend = [](int n) {
	int r = std::log(n);
	while ((1 << r) < n) ++r;
	return r; 
};
auto FMTor = [](std::vector<int> &a, bool isRev) {
	int n = extend(a.size());
	a.resize(1 << n);
	for (int i = 0; i < n; ++i) {
		for (int j = 0; j < (1 << n); ++j) if ((j >> i) & 1) {
			if (isRev) sub(a[j], a[j ^ (1 << i)]);
			else add(a[j], a[j ^ (1 << i)]);
		}
	}
};
auto FMTand = [](std::vector<int> &a, bool isRev) {
	int n = extend(a.size());
	a.resize(1 << n);
	for (int i = 0; i < n; ++i) {
		for (int j = (1 << n) - 1; j >= 0; --j) if ((j >> i) & 1) {
			if (isRev) sub(a[j ^ (1 << i)], a[j]);
			else add(a[j ^ (1 << i)], a[j]);
		}
	}
};
auto FMTxor = [](std::vector<int> &a, bool isRev) {
	int n = extend(a.size());
	a.resize(1 << n);
	for (int i = 0; i < n; ++i) {
		for (int j = (1 << n) - 1; j >= 0; --j) if ((j >> i) & 1) {
			int u = a[j], v = a[j ^ (1 << i)];
			a[j] = (v - u + M) % M;
			a[j ^ (1 << i)] = (u + v) % M;
		}
		if (isRev) for (auto &x : a) x = LL(inv2) * x % M;
	}
};
auto fun = [](std::function<void(std::vector<int> &, bool)> f, std::vector<int> a, std::vector<int> b) {
	int n = extend(std::max(a.size(), b.size()));
	a.resize(1 << n); b.resize(1 << n);
	f(a, 0); f(b, 0);
	std::vector<int> c(1 << n);
	for (int i = 0; i < (1 << n); ++i) c[i] = LL(a[i]) * b[i] % M;
	f(c, 1);
	return c;
};
auto Or = [](std::vector<int> a, std::vector<int> b) {
	return fun(FMTor, a, b);
};
auto And = [](std::vector<int> a, std::vector<int> b) {
	return fun(FMTand, a, b);
};
auto Xor = [](std::vector<int> a, std::vector<int> b) {
	return fun(FMTxor, a, b);
};
// i = j | k and j & k = 0
auto OrAnd = [](std::vector<int> a, std::vector<int> b) {
	int n = extend(std::max(a.size(), b.size()));
	a.resize(1 << n); b.resize(1 << n);
	std::vector<std::vector<int>> sa(n + 1, std::vector<int>(1 << n));
	auto sb = sa, sc = sa;
	for (int i = 0; i < (1 << n); ++i) sa[__builtin_popcount(i)][i] = a[i];
	for (int i = 0; i < (1 << n); ++i) sb[__builtin_popcount(i)][i] = b[i];
	for (int i = 0; i <= n; ++i) {
		FMTor(sa[i], 0);FMTor(sb[i], 0);
	}
	for (int i = 0; i <= n; ++i) {
		for (int j = 0; j <= i; ++j) {
			for (int k = 0; k < (1 << n); ++k) {
				add(sc[i][k], LL(sa[j][k]) * sb[i - j][k] % M);
			}
		}
		FMTor(sc[i], 1);
	}
	std::vector<int> c(1 << n);
	for (int i = 0; i < (1 << n); ++i) c[i] = sc[__builtin_popcount(i)][i];
	return c;
};
}

namespace 真香！

MEX

给定集合 S，求集合 MEX(S) 为 S 中最小未出现的自然数。更进一步的，我们可以求 MEX(S, x) 即 S 中所有数与 x 异或后的 MEX 值。

模板例题：842D

class MEX {
	// 具体数值要对应修改。
	inline static const int B = 20;
	std::array<std::map<int, int>, B> mp;
	std::set<int> S;
public:
	void insert(int x) {
		if (S.count(x)) return;
		S.insert(x);
		int mask = 0;
		for (int i = B - 1; i >= 0; --i) {
			mask |= 1 << i;
			++mp[i][x & mask];
		}
	}
	void erase(int x) {
		if (!S.count(x)) return;
		S.erase(x);
		int mask = 0;
		for (int i = B - 1; i >= 0; --i) {
			mask |= 1 << i;
			--mp[i][x & mask];
		}
	}
	// find mex(a_i ^ x)
	int solve(int x = 0) {
		int mask = 0, r = 0;
		for (int i = B - 1; i >= 0; --i) {
			mask |= x & (1 << i);
			if (mp[i][mask] == (1 << i)) {
				mask ^= 1 << i;
				r |= 1 << i;
			}
		}
		return r;
	}
};

代码解释：mp[i][j]，首先 j 的值域为 $2^i[0, 2^{B-i}]$，即只在乎从 i 到 B 这些位上的值。如果 $mp[i][j] = (1 << i)$ 就说明 j 这一位下面的值都被填满了。在理解 solve 的原理的时候，我们可以先让 x = 0，此时可以看出一丝端倪。即如果 j 这一位被填满了，那么，我们异或 x，也会被填满。

矩阵类

LL mod;
class Matrix {
public:
	const static int N = 102;
	LL a[N][N];
	int n;
	Matrix() {}
	Matrix(int _n, int x = 0): n(_n) { // xIn
		all(0);
		for (int i = 0; i < n; ++i) {
			a[i][i] = x;
		}
	}
	void all(int x) {
		for (int i = 0; i < n; ++i) {
			for (int j = 0; j < n; ++j) {
				a[i][j] = x;
			}
		}
	}
	Matrix operator + (const Matrix & A) const {
		Matrix R(n);
		for (int i = 0; i < n; ++i) {
			for (int j = 0; j < n; ++j) {
				R.a[i][j] = a[i][j] + A.a[i][j];
				if (R.a[i][j] >= mod) R.a[i][j] -= mod;
			}
		}
		return R;
	}
	Matrix operator * (const Matrix & A) const {
		Matrix R(n);
		for (int i = 0; i < n; ++i) {
			for (int k = 0; k < n; ++k) {
				for (int j = 0; j < n; ++j) {
					R.a[i][j] = (R.a[i][j] + a[i][k] * A.a[k][j]) % mod;
				}
			}
		}
		return R;
	}
	void print() {
		for (int i = 0; i < n; ++i) {
			for (int j = 0; j < n; ++j) {
				std::cout << a[i][j] << " ";
			}
			std::cout << std::endl;
		}
	}
};
Matrix pow(Matrix A, int n) {
	Matrix R(A.n, 1);
	while (n) {
		if (n&1) R = R * A;
		n >>= 1; A = A * A;
	}
	return R;
}

注意到，矩阵乘法一定要写成上面的循环形式，这样利用高速缓存执行时间是原有的 $\frac{1}{4}$

另外有序数组的累和要比无序的快很多，也是因为高速缓存（这个不太懂原理）

Gauss 消元法浮点版

#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
using LL = long long;
using pii = std::pair<int, int>;
using pll = std::pair<LL, LL>;

std::vector<double> Gauss(std::vector<std::vector<double>> A, std::vector<double> b) {
	int n = A.size(), m = A[0].size();
	std::vector<double> x(m);
	std::vector<int> p(m);
	std::iota(p.begin(), p.end(), 0);
	const double eps = 1e-12;
	auto findNonZero = [&](int i) { // 实际上找最大的比较好
		for (int row = i; row < n; ++row) if (fabs(A[row][i]) > eps) return row;
		return n;
	};
	auto triangleGauss = [&](int sz) { // A[i][i] = 1
		std::vector<double> x(sz);
		for (int i = sz - 1; i >=0; --i) {
			x[i] = b[i];
			for (int row = 0; row < i; ++row) b[row] -= A[row][i] * x[i];
		}
		x.resize(A[0].size());
		return x;
	};
	int sz = n;
	for (int i = 0, row; i < n; ++i) {
		while (i < m) {
			row = findNonZero(i);
			if (row != n) break;
			for (int j = 0; j < n; ++j) A[j][i] = A[j][m - 1];
			std::swap(p[i], p[--m]);
		}
		if (i == m) {
			for (int row = m; row < n; ++row) if (fabs(b[row])) {
				// std::cout << "\nNo answer\n";
				return std::vector<double>();
			}
			sz = i;
			break;
		}
		if (row != i) {
			std::swap(A[row], A[i]);
			std::swap(b[row], b[i]);
		}
		b[i] /= A[i][i];
		for (int j = m - 1; j >= i; --j) A[i][j] /= A[i][i];
		for (int row = i + 1; row < n; ++row) {
			b[row] -= A[row][i] * b[i];
			for (int j = m - 1; j >= i; --j) {
				A[row][j] -= A[row][i] * A[i][j];
			}
		}
	}
	// if (sz != A[0].size()) std::cout << "\nInfinite answer\n";
	auto xt = triangleGauss(sz);
	for (int t = 0; t < A[0].size(); ++t) x[p[t]] = xt[t];
	return x;
}

Gauss 消元法有限域版

// make sure 0 <= A[i][j], b[i] < M
#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
using LL = long long;
using pii = std::pair<int, int>;
using pll = std::pair<LL, LL>;

std::vector<LL> Gauss(std::vector<std::vector<LL>> A, std::vector<LL> b) {
	int n = A.size(), m = A[0].size();
	std::vector<LL> x(m);
	std::vector<int> p(m);
	std::iota(p.begin(), p.end(), 0);
	const LL M = 998244353;
	std::function<LL(LL)> inv = [&](LL a) -> LL {
		return a == 1 ? 1 : (M - M / a) * inv(M % a) % M; 
	};
	auto sub = [](LL &x, LL y) {
		(x -= y) < 0 && (x += M);
	};
	auto triangleGauss = [&](int sz) { // A[i][i] = 1
		std::vector<LL> x(sz);
		for (int i = sz - 1; i >=0; --i) {
			x[i] = (b[i] + M) % M;
			for (int row = 0; row < i; ++row) sub(b[row], A[row][i] * x[i] % M);
		}
		x.resize(A[0].size());
		return x;
	};
	auto findNonZero = [&](int i) {
		for (int row = i; row < n; ++row) if (A[row][i]) return row;
		return n;
	};
	int sz = n;
	for (int i = 0, row; i < n; ++i) {
		while (i < m) {
			row = findNonZero(i);
			if (row != n) break;
			for (int j = 0; j < n; ++j) A[j][i] = A[j][m - 1];
			std::swap(p[i], p[--m]);
		}
		if (i == m) {
			for (int row = m; row < n; ++row) if (b[row]) {
				// std::cout << "\nNo answer\n";
				return std::vector<LL>();
			}
			sz = i;
			break;
		}
		if (row != i) {
			std::swap(A[i], A[row]);
			std::swap(b[i], b[row]);
		}
		LL inva = inv(A[i][i]);
		(b[i] *= inva) %= M;
		for (int j = m - 1; j >= i; --j) (A[i][j] *= inva) %= M;
		for (int row = i + 1; row < n; ++row) {
			sub(b[row], A[row][i] * b[i] % M);
			for (int j = m - 1; j >= i; --j) {
				sub(A[row][j], A[row][i] * A[i][j] % M);
			}
		}
	}
	// if (sz != A[0].size()) std::cout << "\nInfinite answer\n";
	auto xt = triangleGauss(sz);
	for (int t = 0; t < A[0].size(); ++t) x[p[t]] = xt[t];
	return x;
}

线性规划单纯形法：例题

#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
using LL = long long;
using pii = std::pair<int, int>;
using pll = std::pair<LL, LL>;

using VD = std::vector<double>;
const double eps = 1e-10;
const double inf = 1e10;
// make sure that A = (I, A') and b >= 0, compute max cx
VD simplexCore(VD c, std::vector<VD> A, VD b) {
	int n = A.size(), m = c.size();
	std::vector<int> p(m);
	std::iota(p.begin(), p.end(), 0);
	for (int i = 0; i < n; ++i) A[i].emplace_back(b[i]);
	c.emplace_back(0);
	A.emplace_back(c);
	for (int j = n; j <= m; ++j) {
		for (int i = 0; i < n; ++i) {
			A[n][j] -= A[n][i] * A[i][j];
		}
	}
	auto check = [&]() -> bool {
		for (int j = n; j < m; ++j) if (A[n][j] > eps) {
			bool flag = false;
			for (int i = 0; i < n; ++i) if (A[i][j] > eps) {
				flag = true;
				break;
			}
			if (!flag) return false;
		}
		return true;
	};
	while (1) {
		int ch = std::max_element(A[n].begin() + n, A[n].begin() + m) - A[n].begin(), hc;
		if (A[n][ch] < eps) break;	
		assert(check()); // otherwise unbounded, no max solution
		double theta = DBL_MAX;
		for (int i = 0; i < n; ++i) if (A[i][ch] > eps && A[i].back() / A[i][ch] < theta) {
			theta = A[i].back() / A[i][ch];
			hc = i;
		}
		std::swap(p[ch], p[hc]);
		double tmp = 1 / A[hc][ch];
		for (int j = n; j <= m; ++j) A[hc][j] *= tmp;
		for (int i = 0; i <= n; ++i) if (i != hc) {
			for (int j = n; j <= m; ++j) if (j != ch) {
				A[i][j] -= A[i][ch] * A[hc][j];
			}
		}
		for (int i = 0; i <= n; ++i) A[i][ch] *= -tmp;
		A[hc][ch] = tmp;
	}
	VD x(m);
	for (int i = 0; i < n; ++i) x[p[i]] = A[i].back();
	// watch(-A.back().back()); // max_val
	return x; // point Corresponds to max_val
}
// compute max cx, with Aqx = bq and Alq x <= blq, end of 0 can be ommit in A and Aq
VD simplex(VD c, std::vector<VD> Aq, VD bq, std::vector<VD> Alq, VD blq) {
	assert(Aq.size() == bq.size());
	assert(Alq.size() == blq.size());
	int n = Aq.size() + Alq.size();
	int m = c.size();
	for (int i = 0; i < bq.size(); ++i) if (bq[i] < -eps) {
		for (auto &x : Aq[i]) x = -x;
		bq[i] = -bq[i];
	}
	for (int i = 0; i < blq.size(); ++i) if (blq[i] < -eps) {
		for (auto &x : Alq[i]) x = -x;
		++m;
	}
	std::vector<VD> A(n, VD(n + m));
	VD f(n + m), b(n);
	int now = n + c.size();
	for (int i = 0; i < n; ++i) A[i][i] = 1;
	for (int i = 0; i < Aq.size(); ++i) {
		for (int j = 0; j < Aq[i].size(); ++j) A[i][n + j] = Aq[i][j];
		b[i] = bq[i];
		f[i] = -inf;
	}
	for (int i = 0; i < Alq.size(); ++i) {
		for (int j = 0; j < Alq[i].size(); ++j) A[i + Aq.size()][n + j] = Alq[i][j];
		if (blq[i] < -eps) {
			A[i + Aq.size()][now++] = -1;
			f[i + Aq.size()] = -inf;
		}
		b[i + Aq.size()] = fabs(blq[i]);
	}
	for (int i = 0; i < c.size(); ++i) f[n + i] = c[i];
	auto x = simplexCore(f, A, b);
	return VD(x.begin() + n, x.begin() + n + c.size());
}

差分约束

$n$ 个变量，$m$ 个约束条件，每个约束条件都形如 $x_i - x_j \leq c_k$，此时我们从节点 j 向 i 连一条长度为 $c_k$ 的有向边，（如果有等于号，我们就连两条），设 dist[0] = 0，然后 0 节点向所有节点连一条长度为 0 的有向边。跑单源最短路，如果环中有负环，那么无解，否则 $x_i = dist[i]$ 为一组解。

可用图论中 Bellman-Ford 算法，或 spfa（随笔图跑的快），例题：LOJ P1993，spfa 做法，Bellman-Ford 做法

变式：$\frac{x_i}{x_j} \leq c_k$（取 log 即可）

求$m$阶线性递推关系式第$n$项：时间复杂度$O(m^2 \log n)$，本质是做多项式的带模乘法

#include<bits/stdc++.h>
using LL=long long;
const LL mod = 1e9 + 7;
const int M = 1003;
LL c[M], ans[2 * M];

class RecSeq {
public:
	LL a[2 * M];
	int m;
	RecSeq(int _m, LL x = 0): m(_m) {
		std::memset(a, 0, sizeof (a));
		a[0] = x;
	}
	RecSeq operator * (const RecSeq & A) const {
		RecSeq R(m);
		for (int i = 0; i < m; ++i) {
			for (int j = 0; j < m; ++j) {
				R.a[i + j] = (R.a[i + j] + a[i] * A.a[j]) % mod;
			}
		}
		for (int i = 2 * m - 2; i >= m; --i) {
			for (int j = 0; j < m; ++j) {
				R.a[i - m + j] += (R.a[i] * c[j]) % mod;
			}
			R.a[i] = 0;
		}
		for (int i = 0; i < m; ++i) {
			R.a[i] %= mod;
		}
		return R;
	}
};
template < typename T >
	T tPow(T & A, LL n) {
		T R(A.m, 1);
		while (n) {
			if (n & 1) R = R * A;
			n >>= 1;
			A = A * A;
		}
		return R;
	}
void initC(int m) {
	std::memset(c, 0, sizeof (c));
	c[0] = c[m - 1] = 1;
	for (int i = 1; i < m; ++i) {
		ans[i] = 1;
	}
	ans[m] = 2;
	for (int i = m + 1; i < 2 * m; ++i) {
		ans[i] = ans[i - 1] + ans[i - m];
	}
}
LL getans(LL n, int m) {
	initC(m);
	RecSeq A(m);
	A.a[1] = 1;
	RecSeq R = tPow(A, n - m);
	LL r = 0;
	for (int i = 0; i < m; ++i) {
		r += (R.a[i] * ans[i + m]) % mod;
	}
	return r % mod;
}

利用多项式带模除法：Division with remainder的 $O(m \log m)$ 算法，可优化到 $O(m \log m \log n)$，

但是如果递推关系中仅有常数个不为 0，比如通常是两个，也可以不用多项式带模除法来搞，只需 NTT 就可以优化到 $O(m \log m \log n + m^2)$ （暂时不知道如何去掉 $m^2$）

博弈

数据结构

并查集（Disjoint Set Union）

// 初始情形: p[i] = i
int find(int x) {
	if (p[x] != x) p[x] = find(p[x]);
	return p[x];
}
void merge(int i, int j) {
	// p[find(j)] = p[find(i)];
	// In general we should write below, and merge small to big
	// int fi = find(i), fj = find(j);
	// if (fi != fj) p[fi] = fj;
}

离散化

// 返回的是离散化之后的数组值对应的原始值
template<typename T>
std::vector<T> discrete(std::vector<T>& a) {
	auto b = a;
	std::sort(b.begin(), b.end());
	b.erase(std::unique(b.begin(), b.end()), b.end());
	std::vector<T> r(b.size());
	for (auto & x : a) {
		int id = std::lower_bound(b.begin(), b.end(), x) - b.begin();
		r[id] = x;
		x = id;
	}
	return r;
}

逆序数

直接求 $O(n^2)$ 没啥好写的。
把原数组每个位置进行编号，排序，然后每次把最大的数的编号丢进树状数组中，丢进去先看这个编号前面有多少个数，累加一下就可以了，$O(n^2)$，结合下面树状数组的知识还是很简单的。
带离散化的树状数组（就是如果元素的数值特别大，树状数组内存就不够了，所以需要离散化一下）
归并的求（不会也不想搞 0.0）
逐位处理（代码如下）

// all number in a are nonegetive
LL inverseNumber(std::vector<int>& a) {
	LL r = 0;
	int n = a.size();
	int dn = std::log(*std::max_element(a.begin(), a.end())) + 2;
	for (int d = 1; d <= dn; ++d) {
		std::vector<int> p(n);
		std::iota(p.begin(), p.end(), 0);
		std::sort(p.begin(), p.end(),[&](int i, int j) {
			return (a[i] >> d) < (a[j] >> d) || ((a[i] >> d) == (a[j] >> d) && i < j);
		});
		for (int i = 0, j = 0; i < n; i = j) {
			int x1 = 0;
			while (j < n && (a[p[i]] >> d) == (a[p[j]] >> d)) {
				if ((a[p[j]] >> (d - 1)) & 1) {
					++x1;
				} else {
					r += x1;
				}
				++j;
			}			
		}
	}
	return r;
}

树状数组（单点更新，从 1 到某个位置的最小值）

struct TreeArrayMin {
	std::vector<int> s;
	TreeArrayMin() {}
	TreeArrayMin(int n) : s(n + 1, INT_MAX) {}
	int lowbit(int n) { 
		return n & (-n);
	}
	void modify(int id, int p) {
		while (id < s.size()) {
			s[id] = std::min(s[id], p);
			id += lowbit(id);
		}
	}
	// 计算区间 [1, id] 的最小值
	int min(int id) {
		int r = INT_MAX;
		while (id) {
			r = std::min(r, s[id]);
			id -= lowbit(id);
		}
		return r;
	}
};

树状数组（单点更新，区间求和，编号从 1 开始）

struct TreeArray {
	std::vector<LL> s;
	TreeArray() {}
	TreeArray(int n) : s(n + 1) {}
	int lowbit(int n) { 
		return n & (-n);
	}
	void add(int id, int p) {
		while (id < s.size()) {
			s[id] += p;
			id += lowbit(id);
		}
	}
	LL sum(int id) {
		LL r = 0;
		while (id) {
			r += s[id];
			id -= lowbit(id);
		}
		return r;
	}
	// find minimal index s.t. sum(id) >= x, sum must be increased
	int search(LL val) {
		LL sum = 0;
		int id = 0;
		for (int i = std::__lg(s.size()); ~i; --i) {
			if (id + (1 << i) < s.size() && sum + s[id + (1 << i)] < val) {
				id += (1 << i);
				sum += s[id];
			}
		}
		return ++id;
	}
};

理解了树状数组存储方式就能掌握上述 $O(n)$ 的搜索技巧，并且由此我们可以用树状数组 O(n \log n) 求动态第 k 大的值，并且如果元素很大，那么我们可以离线后离散化再建树状数组。

树状数组（区间更新，区间求和，编号从 1 开始）

有了单点更新的树状数组，只需简单利用差分就可以变成区间的更新了。
设原始数组为 a[1 ~ n]，定义 c[i] = a[i] - a[i - 1], (a[0] = 0) 显然

$\sum_{i = 1}^m a_i = \sum_{i = 1}^m (m - i + 1) c_i = m \sum_{i = 1}^m c_i - \sum_{i = 1}^m (i - 1) c_i$

比如对区间 [l, r] 做更新，那么就只需更新两点：r + 1, l ，套用之前的类就行了。

注意在树状数组中搜索本来应该是 $O(\log ^2 n)$，但是因为在 $2^i$ 的位置搜索时，一步到位。所以复杂度会降到 $O(\log n)$：理论依据

class TreeArrayPlus {
	int n;
	// c[i] = a[i] - a[i - 1], b_i = (i - 1) * c_i
	TreeArray B, C;
	void add(int id, int p) {
		C.add(id, p);
		B.add(id, (id - 1) * p);
	}
public:
	TreeArrayPlus() {}
	TreeArrayPlus(int _n) : n(_n), B(n), C(n) {}
	void add(int l, int r, int p) {
		add(l, p);
		add(r + 1, -p);
	}
	LL sum(int id) {
		return id * C.sum(id) - B.sum(id);
	}
	LL sum(int l, int r) {
		return sum(r) - sum(l - 1);
	}
	// find minimal index s.t. sum(id) >= x, sum must be increased
	int search(LL val) {
		LL sumB = 0, sumC = 0;
		int id = 0;
		for (int i = std::__lg(n); ~i; --i) if (int idi = id + (1 << i); idi <= n) {
			if (idi * (sumC +  C.s[idi]) - B.s[idi] - sumB < val) {
				id = idi;
				sumB += B.s[id];
				sumC += C.s[id];
			}
		}
		return ++id;
	}
};

线段树（正式版）

首先显然总节点 $m$ 上界为 $4n$，并且可以证明 $\frac{m}{n}$ 的上确界为 $4$，下确界为 $2$ 注意到如果 $n = 2^k + 2^{j + 1}$ 时，则 $m = 2 ^{k + 1} + 2^k + \cdots 2^{k - j} + 1$，所以 $\frac{m}{n} = \frac{4 - 2^{-j} + 2^{-k}}{1 + 2^{j + 1 - k}}$，对任意 $\epsilon > 0$ 存在 $j$ 使得 $4 - 2 ^{-j} > 4 - \epsilon$，然后让 $k$ 趋于无穷，那么显然 $\frac{m}{n}$ 上极限为 $4$.（$n = 40$ 时， $\frac{m}{n} > 3$，$n = 2^20 + 2^10 = 1049600$ 时，$\frac{m}{n} > 3.99$）

和与最大值的线段树模板（如果单纯求和，可以用树状数组），现在使用左闭右开线段树，且使用吉老师 pushTag 版本

struct SegmentTree {
	int n;
	std::vector<int> mn, tag;
	std::vector<LL> sm;
	#define lson l, m, 2 * p
	#define rson m + 1, r, 2 * p + 1
	void resize() {
		mn.resize(4 * n);
		tag.resize(4 * n);
		sm.resize(4 * n);
	}
	SegmentTree(int _n) : n(_n) {
		resize();
	}
	SegmentTree(const std::vector<int> &a) {
		n = a.size();
		resize();
		std::function<void(int, int, int)> build = [&](int l, int r, int p) {
			if (l == r) {
				mn[p] = sm[p] = a[l - 1];
				return;
			}
			int m = (l + r) / 2;
			build(lson);
			build(rson);
			pull(p);
		};
		build(1, n, 1);
	}
	void pull(int p) {
		mn[p] = std::min(mn[2 * p], mn[2 * p + 1]);
		sm[p] = sm[2 * p] + sm[2 * p + 1];
	}
	void push(int l, int r, int p) {
		if (tag[p]) {
			int m = (l + r) / 2;
			set(lson, tag[p]);
			set(rson, tag[p]);
			tag[p] = 0;
		}
	}
	void set(int l, int r, int p, int v) {
		tag[p] = mn[p] = v;
		sm[p] = LL(r - l + 1) * v;
	}
	void rangeSet(int L, int R, int v, int l, int r, int p) {
		if (L <= l && R >= r) {
			set(l, r, p, v);
			return;
		}
		int m = (l + r) / 2;
		push(l, r, p);
		if (L <= m) rangeSet(L, R, v, lson);
		if (R > m) rangeSet(L, R, v, rson);
		pull(p);
	}
	// 以下内容根据需要修改
	int query(int x, int& y, int l, int r, int p) {
		if (mn[p] > y) return 0;
		if (x <= l && sm[p] <= y) {
			y -= sm[p];
			return r - l + 1;
		}
		int m = (l + r) / 2;
		push(l, r, p);
		int ans = 0;
		if (x <= m) ans += query(x, y, lson); 
		ans += query(x, y, rson);
		return ans;
	}
	int bounded(int v, int l, int r, int p) {
		if (mn[p] >= v) return r + 1;
		if (l == r) return l;
		int m = (l + r) / 2;
		if (mn[2 * p] >= v) return bounded(v, rson);
		return bounded(v, lson);
	}
	void modify(int x, int y) {
		int l = bounded(y, 1, n, 1);
		if (l <= x) rangeSet(l, x, y, 1, n, 1);
	}
	int query(int x, int y) {
		return query(x, y, 1, n, 1);
	}
};

代码参考了 Jiangly 的模板，他的是左闭右开的。

线段树（弃用版）

#define lrt rt << 1
#define rrt rt << 1 | 1
#define lson l, m, lrt
#define rson m + 1, r, rrt
const int N = 1e5 + 5;
LL sum[N * 3], col[N * 3];
void pushUp(int rt) {
	sum[rt] = sum[lrt] + sum[rrt];
}
void pushDown(int rt, int m) {
	if (col[rt]) {
		col[lrt] += col[rt];
		col[rrt] += col[rt];
		sum[lrt] += (m - (m >> 1)) * col[rt];
		sum[rrt] += (m >> 1) * col[rt];
		col[rt] = 0;
	}
}
void build(int l, int r, int rt) {
	if (l == r) {
		scanf("%I64d", & sum[rt]);
		return;
	}
	col[rt] = 0;
	int m = (l + r) >> 1;
	build(lson);
	build(rson);
	pushUp(rt);
}
void update(int L, int R, int p, int l, int r, int rt) {
	if (L <= l && R >= r) {
		sum[rt] += p * (r - l + 1);
		col[rt] += p;
		return;
	}
	pushDown(rt, r - l + 1);
	int m = (l + r) >> 1;
	if (L <= m) update(L, R, p, lson);
	if (R > m) update(L, R, p, rson);
	pushUp(rt);
}
LL query(int L, int R, int l, int r, int rt) {
	if (L <= l && R >= r) return sum[rt];
	pushDown(rt, r - l + 1);
	int m = (l + r) >> 1;
	LL ans = 0;
	if (L <= m) ans += query(L, R, lson);
	if (R > m) ans += query(L, R, rson);
	return ans;
}

吊打线段树的珂朵莉树（Chtholly Tree）

RMQ 求区间最大值（弃用，RMQ 和 spfa 已经死了）

class RMQ {
	static
	const int N = 100005;
	public:
		int n;
	int a[N][20]; //a[i][j]=max(s[i]---s[i+2^j-1])
	RMQ(int * s, int _n): n(_n) {
		for (int i = 0; i != n; ++i) a[i][0] = s[i];
		int len = (int)(log(double(n)) / log(2.0)) + 2;
		for (int i = 1; i != len; ++i) {
			for (int j = 0; j <= n - (1 << i); ++j) {
				a[j][i] = max(a[j][i - 1], a[j + (1 << (i - 1))][i - 1]);
			}
		}
	}
	int Max(int l, int r) { // 0 <= l <= r < n
		int len = (int)(log(double(r - l + 1)) / log(2.0));
		return max(a[l][len], a[r - (1 << len) + 1][len]);
	}
};

我们经常用二分答案的思想，但是其实二分答案是仅仅知道其单调的情况下的策略，实际上，对于具体的问题，我们完全可以对 $m$ 的值进行不同的处理，而非单纯的 $m = (l + r) >> 1 $。

三分法简单版

// for a given covex function f: (f(a) + f(b)) / 2 >= f((a + b) / 2)
int l = 1, r = 1e9;
while (r > l) {
	int m = (r - l) / 3;
	int lm = l + m, rm = r - m;
	if (f(lm) < f(rm)) r = rm - 1;
	else l = lm + 1;
}
return l;

标准三分法用黄金分割的原因

我们不妨设原始区间为 [0, 1]，我们在其中选两个点 0 < a < b < 1，然后比较 f(a) 和 f(b)，然后再相应改变区间。然后重复上述过程。如果我们能充分利用计算过的值，也就是说假设更新后的区间为 [0, b] 那么我们自然想让 a 的计算值充分被利用，所以我们想再选的两个点的其中一个是 a，如果更新后区间为 [a, 1] 同理。也就是说我们有策略

$\frac{a}{b} = b, \frac{b - a}{1 - a} = a$

化简可得 $b(1 + b) = 1$，即 $b = \frac{\sqrt{5} - 1}{2}, a = b ^ 2 = \frac{3 - \sqrt{5}}{2} = 1 - b$。

注意到上述 $b$ 的值正好是黄金分割 0.618…

三分法黄金分割版（可用版）)

// for a given covex function f: (f(a) + f(b)) / 2 >= f((a + b) / 2)
auto tupleCut = [&](int l, int r) {
	const double phiL = (3 - std::sqrt(5)) / 2;
	const double phiR = (std::sqrt(5) - 1) / 2;
	auto getLeft = [&](int l, int r) -> int {
		return l + (r - l) * phiL;
	};
	auto getRight = [&](int l, int r) -> int {
		return l + (r - l) * phiR;
	};
	int lm = getLeft(l, r), rm = getRight(l, r);
	LL fl = f(lm), fr = f(rm);
	while (r - l > 2) {
		if (fl < fr) {
			r = rm;
			rm = lm;
			fr = fl;
			fl = f(lm = getRight(l, rm));
		} else {
			l = lm;
			lm = rm;
			fl = fr;
			fr = f(rm = getLeft(lm, r));
		}
	}
	fl = f(l), fr = f(r);
	while (l < r) {
		if (fl < fr) fr = f(--r);
		else fl = f(++l);
	}
	return fr;
};

注意返回的是最值而不是最值点

注意到一定要用上述写法，由于取整带来的误差，所以必须充分利用“左分点”是“右分点”的“右分点”，“右分点”是“左分点”的“左分点”（用来保证 $l \leq lm \leq rm \leq r$），然后如果单次求 $f$ 的复杂度特别高，两段各自最后一步也是可以优化一下，少算一次 $f$，但是会很不优雅，所以还是算了。

subarray（连续部分） VS subsequence（不要求连续部分）

最长（严格）递增子序列

下面数组 b 的意义：b[i] 表示当前长为 i 的子列的末位的最小值。
这个子序列还能给附加限制（比如都要大于某个值）

int LIS(std::vector<int>& a) { // length of longest increasing subsquence
	std::vector<int> b;
	for (auto x : a) {
		if (auto it = std::lower_bound(b.begin(), b.end(), x); it == b.end()) {
			b.emplace_back(x);
		} else *it = x;
	}
	return b.size();
}
int LNDS(std::vector<int>& a) { // length of longest increasing subsquence
	std::vector<int> b;
	for (auto x : a) {
		if (auto it = std::upper_bound(b.begin(), b.end(), x); it == b.end()) {
			b.emplace_back(x);
		} else *it = x;
	}
	return b.size();
}
auto LISP(std::vector<int>& a) { // longest increasing subsquence
	std::vector<int> b, pb, pa(a.size());
	std::iota(pa.begin(), pa.end(), 0);
	for (int i = 0; i != a.size(); ++i) {
		if (auto it = std::upper_bound(b.begin(), b.end(), a[i]); it == b.end()) {
			if (!pb.empty()) pa[i] = pb.back();
			b.emplace_back(a[i]);
			pb.emplace_back(i);
		} else {
			*it = a[i];
			int t = it - b.begin();
			pb[t] = i;
			if (t > 0) pa[i] = pb[t - 1];
		}
	}
	std::stack<int> c;
	int now = pb.back();
	c.push(a[now]);
	while (now != pa[now]) {
		now = pa[now];
		c.push(a[now]);
	}
	return c;
}
// lower_bound(first,end,val) 表示在单增 [frist,end) 中首次大于等于 val 的位置
// upper_bound(first,end,val) 表示在单增 [frist,end) 中首次大于 val 的位置

背包

const int MAX = 1e5 + 5;
int r[MAX];
void pack(int cash, int num, int v, int w) {
	if (num == 0 || w == 0 || v == 0) return;
	if (num == 1) { //  0-1背包
		for (int i = cash; i >= v; --i)
			r[i] = max(r[i], r[i - v] + w);
		return;
	}
	if (num * v >= cash - v + 1) { //完全背包
		for (int i = v; i <= cash; ++i)
			r[i] = max(r[i], r[i - v] + w);
		return;
	}
	int q[MAX], s[MAX], head, tail;
	for (int j = 0; j < v; ++j) { //多重背包
		q[0] = r[j];
		s[0] = head = tail = 0;
		for (int i = 1, k = j + v; k <= cash; ++i, k += v) {
			q[i] = r[k] - i * w;
			while (s[head] < i - num) ++head;
			while (head <= tail && q[tail] < q[i]) --tail;
			s[++tail] = i;
			q[tail] = q[i];
			r[k] = q[head] + i * w;
		}
	}
}

堆与 STL 优先队列

可以使用 C++STL 的 priority_queue，查找可用 lower_bound 和 upper_bound。C++ STL 中优先队列是用堆来实现的。用途十分广泛，例如加速最小生成树，拓扑排序，等等。
堆的实现一般是用数组。我们可以用 1 作为树的根，对每一个节点 $x$，它的两个节点分别就是 $2x$ 和 $2x + 1$ 平时都用 x << 1, x << 1 | 1 表示。堆只支持三个操作:

插入一个节点(我们实现时是插入最尾部，这样保证了是一个完全二叉树) $O(\log n)$
删除最大键值节点（删除根元素的值） $O(\log n)$
输出最大键值节点（查看根元素的值） $O(1)$

我掌握的也很喜欢的几个数据结构：单调队列，单调栈，树状数组，堆，线段树，红黑树

单调队列：解决滑动窗口问题（固定长度内的最值问题）

知乎 Pecco 讲的很好（建议直接去看它的讲解）：

如果一个选手比你小还比你强，你就可以退役了。——单调队列的原理

// 求每个长度为 m 的区间最大值的编号
std::vector<int> monicDequeMax(std::vector<int> &a, int m) {
	std::vector<int> r;
	std::deque<int> Q;
	for (int i = 0; i < a.size(); ++i) {
		if (!Q.empty() && i - Q.front() >= m) Q.pop_front();
		// 如果求最小值，大于号改成小于号即可
		while (!Q.empty() && a[i] > a[Q.back()]) Q.pop_back();
		Q.push_back(i);
		// 如果需要返回值，就在下面加入 a[Q.front()]
		if (i >= m - 1) r.emplace_back(Q.front());
	}
	return r;
}

例题：LOJ P2216：这个是二维的，我们可以一维一维的处理

单调队列优化 DP

例题：LOJ P2034：取数字使得和最大，但是不能取连续 k 个数

肯定是 dp 问题，如果把 dp[i] 定义成取第 i 个，前 i 个结果的最值，会发现很搞。
因此我们反过来考虑。考虑删除若干个数，且删除的间隔不超过 k，求删除的最小和。最终答案就是总和减去最小和。设 dp[i] 表示删除 i，且满足性质的前 i 个数的答案。那么显然 $dp[i] = a[i] i \leq k$，$dp[i] = a[i] + \min_{i - k \leq j \leq i - 1} dp[j]$。注意最终答案不是总和减去 dp 的最小值，而是 $dp[n - k - 2, \cdots, n - 1]$ 的最小值。

#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
using LL = long long;
int main() {
	// freopen("in", "r", stdin);
	std::ios::sync_with_stdio(false);
	std::cin.tie(nullptr);
	int n, k;
	std::cin >> n >> k;
	std::vector<int> a(n), dp(n);
	LL s = 0;
	for (auto &x : a) {
		std::cin >> x;
		s += x;
	}
	std::deque<int> Q;
	for (int i = 0; i < n; ++i) {
		dp[i] = a[i];
		if (i >= k + 1) dp[i] += dp[Q.front()];
		if (!Q.empty() && i - Q.front() >= k + 1) Q.pop_front();
		while (!Q.empty() && dp[i] <= dp[Q.back()]) Q.pop_back();
		Q.push_back(i);
	}
	std::cout << s - *std::min_element(dp.end() - k - 1, dp.end());
	return 0;
}

单调栈：形式更简单应用更广

知乎 Pecco 的精彩讲解：维护一个栈，当我们碰上一个新元素，我们知道，越靠近栈顶的元素离新元素位置越近。所以不断比较新元素与栈顶，如果新元素比栈顶大，则可断定新元素就是栈顶的下一个更大元素，于是弹出栈顶并继续比较。直到新元素不比栈顶大，再将新元素压入栈。显然，这样形成的栈是单调递减的。

应用一：求下一个比自身大的元素位置（下可以改成上，大可以改成小）

洛谷模板题：LOJ P5788

应用二：两个元素间所有元素均（不）大/小于这二者。

洛谷进阶题：LOJ P1823，问有多少对元素，它们之间没有比它们都大的元素。
代码放在这里了

单调栈优化 DP

应用一优化 DP 例题：1313C2，首先最优答案肯定时先递增后递减的。相当于有一个制高点，枚举制高点，自然有 $O(n^2)$ 的算法。但是可以优化到 $O(n)$
代码放在这里了

应用二优化 DP 例题：1407D，每次跳跃，它们之间的元素都严格大于它们或者严格小于它们。首先设 dp[i] 为到达 i 最小跳跃数，那么显然 $\displaystyle dp[i] = \min_{j \to i} dp[j] + 1$。我们可以用两个单调栈来看那些 j 能跳到 i。这里了

红黑树 red-black tree

class RBT {
	typedef int elemType;
	struct Node {
		const static bool RED = 0;
		const static bool BLACK = 1;
		Node * ch[2], * fa; // x->fa->ch[x->rs] = x
		int sz;
		elemType key;
		bool color, rs; // is rightson
	};
	Node * root; // root has no father
	void faSon(Node * x, Node * y, bool rs) {
		y->fa = x;
		y->rs = rs;
		x->ch[rs] = y;
	}
	void newNode(Node * x, elemType val, bool rs) { // x->ch[rs]=null
		Node * p = new Node;
		p->ch[0] = p->ch[1] = null;
		p->sz = 1;
		p->key = val;
		p->color = Node::RED;
		faSon(x, p, rs);
		++x->sz;
	}
	void rotate(Node * x, bool rs) { // x must not null
		Node * y = x->ch[!rs];
		if (y == null) return;
		if (x == root) root = y;
		else faSon(x->fa, y, x->rs);
		faSon(x, y->ch[rs], !rs);
		faSon(y, x, rs);
		y->sz = x->sz;
		x->sz = x->ch[0]->sz + x->ch[1]->sz + 1;
	}
	void insMaintain(Node * x) { // x->color is RED
		if (x == root || x->fa->color == Node::BLACK) return;
		if (x->fa->fa->ch[!x->fa->rs]->color == Node::BLACK) {
			if (x->rs ^ x->fa->rs) rotate(x->fa, x->fa->rs);
			else x = x->fa;
			x->color = Node::BLACK;
			x->fa->color = Node::RED;
			rotate(x->fa, !x->rs);
		} else {
			x = x->fa->fa;
			x->color = Node::RED;
			x->ch[0]->color = x->ch[1]->color = Node::BLACK;
			insMaintain(x);
		}
	}
	void delCase1(Node * x, Node * y) {
		y->color = Node::BLACK;
		y->fa->color = Node::RED;
		y = y->ch[!y->rs];
		rotate(x->fa, x->rs);
		delCase2(x, y);
	}
	void delCase2(Node * x, Node * y) {
		if (y->ch[y->rs]->color == Node::BLACK) {
			if (y->ch[!y->rs]->color == Node::BLACK) {
				y->color = Node::RED;
				delMaintain(y->fa);
			} else {
				y->color = Node::RED;
				y->ch[!y->rs]->color = Node::BLACK;
				rotate(y, y->rs);
				delCase3(y->fa);
			}
		} else delCase3(y);
	}
	void delCase3(Node * y) {
		y->color = y->fa->color;
		y->ch[y->rs]->color = y->fa->color = Node::BLACK;
		rotate(y->fa, !y->rs);
	}
	void delMaintain(Node * x) {
		if (x == root || x == null) return;
		if (x->color == Node::RED) {
			x->color = Node::BLACK;
			return;
		}
		Node * y = x->fa->ch[!x->rs];
		if (y->color == Node::RED) delCase1(x, y);
		else delCase2(x, y);
	}
	Node * pred(Node * x, elemType val) { // max elem <= val
		if (x == null) return null;
		if (x->key > val) return pred(x->ch[0], val);
		else if (x->ch[1] == null) return x;
		return pred(x->ch[1], val);
	}
	Node * succ(Node * x, elemType val) { // min elem >= val
		if (x == null) return null;
		if (x->key < val) return succ(x->ch[1], val);
		else if (x->ch[0] == null) return x;
		return succ(x->ch[0], val);
	}
	int rank(Node * x, elemType val) { // count elem <= val
		if (x == null) return 0;
		if (x->key > val) return rank(x->ch[0], val);
		return x->ch[0]->sz + 1 + rank(x->ch[1], val);
	}
	Node * select(Node * x, int k) { // k-th smallest elem
		if (x == null || x->sz < k) return null;
		int sz = x->ch[0]->sz + 1;
		if (sz == k) return x;
		if (sz < k) return select(x->ch[1], k - sz);
		return select(x->ch[0], k);
	}
	void clear(Node * x) {
		if (x->ch[0] != null) clear(x->ch[0]);
		if (x->ch[1] != null) clear(x->ch[1]);
		if (x != null) delete x;
	}
	void print(Node * x) {
		printf("key = %d, sz = %d ", x->key, x->sz);
		puts(x->color == Node::RED ? "RED" : "BLACK");
		if (x->ch[0] != null) print(x->ch[0]);
		if (x->ch[1] != null) print(x->ch[1]);
	}
	public:
		const static int INF = 0x3f3f3f3f;
	Node * null;
	RBT() {
		null = new Node; // no key, rs, father
		null->ch[0] = null->ch[1] = null;
		null->sz = 0;
		null->color = Node::BLACK;
		root = null;
		null->key = INF; // for convenient
	}
	Node * search(elemType val) {
		Node * x = root;
		while (x != null) {
			if (val == x->key) return x;
			x = x->ch[val >= x->key];
		}
		return null;
	}
	void insert(elemType val) {
		if (root == null) {
			root = new Node; // no father, rs
			root->ch[0] = root->ch[1] = null;
			root->sz = 1;
			root->color = Node::BLACK;
			root->key = val;
			return;
		}
		Node * x = root;
		while (x->ch[val >= x->key] != null) {
			++x->sz;
			x = x->ch[val >= x->key];
		}
		newNode(x, val, val >= x->key);
		insMaintain(x->ch[val >= x->key]);
		root->color = Node::BLACK;
	}
	void del(elemType val) {
		Node * x = search(val), * y;
		if (x == null) return;
		while (x->ch[0] != null || x->ch[1] != null) {
			bool rs = 0;
			if (x->ch[rs] == null) rs = !rs;
			y = x->ch[rs];
			while (y->ch[!rs] != null) y = y->ch[!rs];
			std::swap(x->key, y->key);
			x = y;
			if (x->color == Node::RED) break;
		}
		delMaintain(x);
		root->color = Node::BLACK;
		y = x;
		while (y != root) {
			y = y->fa;
			--y->sz;
		}
		if (x == root) root = null;
		else if (x->ch[0] != null) faSon(x->fa, x->ch[0], x->rs);
		else if (x->ch[1] != null) faSon(x->fa, x->ch[1], x->rs);
		else x->fa->ch[x->rs] = null;
		delete x;
	}
	elemType pred(elemType val) {
		return pred(root, val)->key;
	}
	elemType succ(elemType val) {
		return succ(root, val)->key;
	}
	int rank(elemType val) {
		return rank(root, val);
	}
	elemType select(int k) {
		return select(root, k)->key;
	}
	void clear() {
		clear(root);
		root = null;
	}
	int size() {
		return root->sz;
	}
	void print() {
		if (root != null) print(root);
	}
	int getans(int a, int k) { // for particular use
		return select(root, rank(root, a) + k)->key;
	}
};

几何

二维凸包正式版

using Point = std::pair<double, double>;
bool crossLeft(const Point &op, const Point &sp, const Point &ep) {
	return (sp.first - op.first) * (ep.second - op.second) 
	<= (sp.second - op.second) * (ep.first - op.first);
}
std::vector<Point> convexHull(std::vector<Point> p) {
	std::sort(p.begin(), p.end());
	p.erase(std::unique(p.begin(), p.end()), p.end());
	int n = p.size();
	std::vector<Point> q(n + 1);
	int top = 0;
	for (int i = 0; i < n; ++i) {
		while (top > 1 && crossLeft(q[top - 1], p[i], q[top - 2])) --top;
		q[top++] = p[i];
	}
	int len = top;
	for (int i = n - 2; i >= 0; --i) {
		while (top > len && crossLeft(q[top - 1], p[i], q[top - 2])) --top;
		q[top++] = p[i];
	}
	top -= n > 1;
	q.resize(top);
	return q;
}

三维凸包模板可以参考上交模板

旋转卡壳求凸包直径（即平面上点的最远欧式距离）

整数版模板例题：LOJ P1452
实数版模板例题：LOJ P6247

using Point = std::pair<double, double>;
double cross(const Point &op, const Point &sp, const Point &ep) {
	return (sp.first - op.first) * (ep.second - op.second) 
	- (sp.second - op.second) * (ep.first - op.first);	
}
double dist2(const Point &p, const Point &q) {
	double x = q.first - p.first, y = q.second - p.second;
	return x * x + y * y;
};
double diameter(std::vector<Point> p) {
	auto q = convexHull(p);
	if (q.size() < 2) return 0;
	if (q.size() == 2) return dist2(q[0], q[1]);
	int n = q.size();
	q.emplace_back(q[0]);
	LL ans = 0;
	for (int i = 0, j = 2; i < n; ++i) {
		while (cross(q[i], q[i + 1], q[j]) < cross(q[i], q[i + 1], q[j + 1])) j = (j + 1) % n;
		ans = std::max({ans, dist2(q[i], q[j]), dist2(q[i + 1], q[j])});
	}
	return std::sqrt(ans);
}

Manhattan 距离的话仅需考虑 $x_i + y_i$ 的最大最小值之差与 $x_i - y_i$ 的最大最小值之差的最大值即可。

分治法求平面最短距离（任何距离都适用）

首先根据横坐标排序，然后取中位数假设处理好了左右两边的值，然后合并中间的值，首先距离中心点的横坐标不能超过已知的最小值，然后把筛出来的点按照纵坐标排序，然后 $O(n)$ 更新答案。总题复杂度 $O(n \log^2 n)$，如果使用归并排序理论复杂度为 $O(n \log n)$，但是实际效果并不如直接排序。

例题：[https://www.luogu.com.cn/problem/P1429] 和 LOJ P6247

using Point = std::pair<double, double>;
// 这里不要用 dist2，否则很多比较的地方都要平方，反而不优雅了。
double dist (const Point& p, const Point &q) {
	double x = q.first - p.first, y = q.second - p.second;
	return std::sqrt(x * x + y * y);
};
double minDist(std::vector<Point> a) {
	double d = DBL_MAX;
	int n = a.size();
	if (n <= 1) return d;
	std::sort(a.begin(), a.end());
	std::function<void(int, int)> merge = [&](int l, int r) {
		if (r - l <= 1) return;
		if (r - l == 2) {
			d = std::min(d, dist(a[l], a[l + 1]));
			return;
		}
		int m = (l + r) / 2;
		merge(l, m);
		merge(m + 1, r);
		std::vector<Point> p;
		for (int i = m - 1; i >= l && a[m].first - a[i].first < d; --i) {
			p.emplace_back(a[i].second, a[i].first);
		}
		for (int i = m; i < r && a[i].first - a[m].first < d; ++i) {
			p.emplace_back(a[i].second, a[i].first);
		}
		std::sort(p.begin(), p.end());
		for (int i = 0; i < p.size(); ++i) {
			for (int j = i + 1; j < p.size() && p[j].first - p[i].first < d; ++j) {
				d = std::min(d, dist(p[i], p[j]));
			}
		}
	};
	merge(0, n);
	return d;
}

分治法还能求两个点作为对交的矩阵的最大面积（详见动态规划之分治算法）

三维偏序之陈丹琪分治（仅支持离线查询，还是直接暴力好用~）

如果带更新怎么处理呢？先预处理求出，之后更新一个计算一个更新？（那也不太行呀）

一般地，我们考虑可以考虑 $k$ 维偏序，设有 $n$ 个 $k$ 维向量，$a_j \leq a_i$ 当且仅当所有的下标都满足小于等于关系，想知道对任意 $i$ 有多少个 $j \neq i$ 使得 $a_j \leq a_i$。

有复杂度 $O(n \log^k n)$ 的算法，因此在 $k > 3$ 时，我们会选择直接 $O(n^2)$ 暴力解决问题（见下小节）。

$k = 1$ 时，我们直接排序，假设没有相同元素，那么它们排完序之后的位置就是答案，有相同的数字的话可以先合并，也可以用 upper_bound 查找出结果。复杂度 $O(n \log n)$
$k = 2$ 时，我们先对第一个坐标偏序，再来一个树状数组，一个个的加入元素，加入之前可以查询结果。这也是求逆序数的操作(如果数据值域范围很大，可以离散化处理一下，仅需对要加入树状数组的那一维离散化，排序可以使用下标排序，就可以避免使用 tuple）。

因此三维偏序是一个空缺的问题，就有大名鼎鼎的 cdq 分治。

模板例题：LOJ P3810，这个题的题解中，有人讲的很好，echo6342：

cdq 分治每次计算前一半对后一半的影响。具体地，假设三维分别是 $x, y, z$，先按 $x$ 排序。分治时每次将前半边、后半边分别按 $y$ 排序。虽然现在 $x$ 的顺序被打乱了，但是前半边还是都小于后半边的，所以要是只计算前半边对后半边的偏序关系，是不会受到 $x$ 的影响的。维护后一半的指针 i，前一半的指针 j，每次将 i 后移一位时，若 $y[j] \leq y[i]$ 则不断后移 j，并不断将 z[j] 加入树状数组。然后再查询树状数组中有多少数小于等于 z[i]。最后要清空树状数组（注意清空的时候不能直接清空，而是根据更新的命令，反向一次命令来清空，否则一直开树状数组耗时的），还有就是要去重贼麻烦，还是弃用吧。

struct Node {
	int x, y, z, id, w;
	bool operator<(const Node &A) const {
		if (x == A.x) return y == A.y ? z < A.z : y < A.y;
		return x < A.x;
	}
};
// ans[i] 表示 小于或等于 a[i] 的元素个数
std::vector<int> cdq(std::vector<Node> &a, int k) {
	// 先按照 y 排序，免得后面代码写的太麻烦
	std::vector<int> ans(a.size());
	std::sort(a.begin(), a.end());
	// 去重操作
	int last = 0;
	for (int i = 1; i < a.size(); ++i) {
		if (a[i].x != a[i - 1].x || a[i].y != a[i - 1].y || a[i].z != a[i - 1].z) {
			int t = i - last - 1;
			for (int j = last; j < i; ++j) {
				ans[a[j].id] = t;
				a[j].w = 0;
			}
			a[i - 1].w = i - last;
			last = i;
		}
	}
	int t = a.size() - last - 1;
	for (int i = last; i < a.size(); ++i) {
		ans[a[i].id] = t;
		a[i].w = 0;
	}
	a.back().w = a.size() - last;
	TreeArray A(k);
	auto cmpy = [](const Node &lhs, const Node &rhs) {
		return lhs.y < rhs.y;
	};
	std::function<void(int, int)> divide = [&](int l, int r) {
		if (r - l <= 1) return;
		int m = (l + r) / 2;
		divide(l, m);
		divide(m, r);
		std::sort(a.begin() + l, a.begin() + m, cmpy);
		std::sort(a.begin() + m, a.begin() + r, cmpy);
		int t = l;
		for (int i = m; i < r; ++i) {
			while (t < m && a[t].y <= a[i].y) {
				A.add(a[t].z, a[t].w);
				++t;
			}
			ans[a[i].id] += A.sum(a[i].z);
		}
		for (int i = l; i < t; ++i) A.add(a[i].z, -a[i].w);
	};
	divide(0, a.size());
	return ans;
}

k 维偏序（暴力 bitset 优化，分块时间换空间） $O(\frac{k n^2}{w})$

const int N = 4e4 + 2;
// a 是 k * n 矩阵表示 n 个 k 维向量，输出每个小于自身的向量个数
std::vector<int> partialOrder(std::vector<std::vector<int>> &a) {
// 直接暴力不太行，所以需要时间换空间，具体说就是分块。
	int k = a.size(), n = a[0].size();
	using Node = std::vector<std::pair<int, int>>;
	std::vector<Node> f(k, Node(n));
	for (int i = 0; i < k; ++i) {
		for (int j = 0; j < n; ++j) f[i][j] = {a[i][j], j};
		std::sort(f[i].begin(), f[i].end());
	}
	int sn = std::sqrt(n);
	using Data = std::vector<std::bitset<N>>;
	std::vector<Data> bs(k, Data(n / sn + 1));;
	for (int i = 0; i < k; ++i) {
		std::bitset<N> now;
		for (int j = 0; j < n; ++j) {
			if (j % sn == 0) bs[i][j / sn] = now;
			now.set(f[i][j].second);
		}
		if (n % sn == 0) bs[i][n / sn] = now;
	}
	auto getbst = [&](int i, int val) -> std::bitset<N> {
		// 如果求小于或等于的个数，这里要改成 upper_bound 并且要用 INT_MAX，还有最终答案减 1（去掉自身）
		int j = std::lower_bound(f[i].begin(), f[i].end(), std::make_pair(val, INT_MIN)) - f[i].begin();
		std::bitset<N> r = bs[i][j / sn];
		for (int t = j / sn * sn; t < j; ++t) r.set(f[i][t].second);
		return r;
	};
	std::vector<int> r(n);
	for (int j = 0; j < n; ++j) {
		std::bitset<N> now; now.set();
		for (int i = 0; i < k; ++i) {
			now &= getbst(i, a[i][j]);
		}
		r[j] = now.count();
	}
	return r;
}

模板例题：LOJ U66865，可参考小蒟蒻 yyb 的博客中的 ppt 实现。其它例题：偏序++

虽然三维偏序问题用 cdq 分治更好，但是用 bitset 暴力过题还是没啥问题的，例如 LOJ P3810

几类根号算法（原来这个叫整除分块）

首先注意到 $\max \{x \mid \lfloor \frac{n}{x} \rfloor = \lfloor \frac{n}{i} \rfloor \} = \lfloor \frac{ n }{ \lfloor \frac{n}{i} \rfloor } \rfloor$

Proof：由于 $x \cdot \lfloor \frac{n}{x} \rfloor \leq n$ 所以 $x \leq \frac{n}{\lfloor \frac{n}{x} \rfloor}$，所以 $x \leq \lfloor \frac{n}{\lfloor \frac{n}{x} \rfloor} \rfloor$，所以取 $x = i$ 和 $x = \lfloor \frac{n}{i} \rfloor$，则 $\lfloor \frac{n}{\lfloor \frac{ n }{ \lfloor \frac{n}{i} \rfloor } \rfloor} \rfloor = \lfloor \frac{n}{i} \rfloor$，另一方面，若 $\lfloor \frac{n}{x} \rfloor = \lfloor \frac{n}{i} \rfloor$, 则 $x \leq \lfloor \frac{n}{\lfloor \frac{n}{i} \rfloor} \rfloor$，证毕。

$ s(n) = \sum_{i=1}^{n} \lfloor \frac{n}{i} \rfloor $

LL getsum(int n) {
	LL sum = 0;
	for (int i = 1, j; i <= n; i = j + 1) {
		j = n / (n / i);
		sum += LL(j - i + 1) * (n / i);
	}
	return sum;
}

$\sum_{i=1}^n \lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor f(i)$

LL getsum(int n, int m) { // g[i]=f[i]+g[i-1]
	LL sum = 0;
	for (int i = 1, j; i <= min(n, m); i = j + 1) {
		j = min(n / (n / i), m / (m / i);
		sum += LL(n / i) * (m / i) * (g[j] - g[i - 1]);
	}
	return sum;
}

例题：牛客，大致题意是：$n=p \times k + m, 0 \leq m < k$，求 $\sum_{p = 1}^n km$，注意到 $m = n - pk$，所以用上述公式计算即可。

再来一题：武汉理工

gcd 经典推导 1, 2

$h(n) = \frac{n(n-1)(n-2)}{3} - \sum_{i=2}^n h(\lfloor \frac{n}{i} \rfloor)$

std::map<int, int> mp;
int getans(int n) {
	auto it = mp.find(n);
	if (it != mp.end()) return it->second;
	int r = LL(n) * (n - 1) % M * (n - 2) % M * inv3 % M;
	for (int i = 2, j; i <= n; i = j + 1) {
		j = n / (n / i);
		r -= LL(j - i + 1) * getans(n / i) % M;
		if (r < 0) r += M;
	}
	mp.insert(std::make_pair(n, r));
	return r;
}

递归算法复杂度分析

下图取自算法导论