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

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

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

此篇代码不再更新,已经发现不少 bug。下面是代码存储位置(并不包含全部此篇内容):

通用代码块

竞赛入门书籍推荐:GuideToCompetitiveProgramming

佛祖保佑,永无 bug

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
// #define print(x) std::cout << (x) << std::endl
// #define println std::cout << std::endl
using LL = long long;

int main(){
//freopen("in","r",stdin);
std::ios::sync_with_stdio(false);std::cin.tie(nullptr);
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 │ . │←─┘│
* └─────┴────┴────┴───────────────────────┴────┴────┴────┴────┘ └───┴───┴───┘ └───────┴───┴───┘
*/

VScode 下 debeg

需要在默认生成的 task.json 中在 "args": 下添加一行 "-std=c++17", 不然要报 warning 很烦。另外要使用绝对路径。

1
2
// 不能使用相对路径 debug,也可能是设置问题
freopen("C:/Users/dna049/cf/in", "r", stdin);

递归程序防止爆栈

在 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

1
2
3
4
5
6
7
8
9
# 多 case 输入
for _ in range(int(input())):
# 单行输入
n = int(input())
# 两个元素一行输入
a, b = map(int, input().split())

# 提前结束
exit()

产生 log 的几个原因

  1. 二分,三分
  2. $1 + \frac{1}{2} + \cdots \frac{1}{n} \sim \log n$
  3. 树状数组,线段树
  4. 重链剖分
  5. 倍增(太好用了,特别是 nxt 数组进行加速)
  6. 莫队(离线算法)

产生根号的几个原因

  1. 朴素判断素数
  2. $\lfloor \frac{n}{i} \rfloor$ 的值域是 $O(\sqrt{n})$ 的
  3. 网络流中 HLPP(没读过这篇复杂度分析的论文,不懂)
  4. 分块处理

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#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")

示例:

莫队加速

位运算的关系

  • 异或 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

最大最小值分配律

例题: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
}

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

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

1
2
3
4
5
6
7
8
9
10
11
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 的使用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
__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;
}

交互式题目模板

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#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 加一个标号,实现上述操作

优雅的输出技巧

1
2
3
4
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$

1
2
3
4
5
6
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$ 相交的长度)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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;
}

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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$ 的取值还是需要注意的。

1
2
3
4
5
6
7
8
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;
};

输出全排列

1
2
3
4
5
6
7
8
9
10
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()));
}

输出全排列的原理

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
// 简洁写法,不推荐,推荐使用内建 __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)$

所以,我们可以在 $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$,所以复杂度就一致了。当然也可以直接化简成:

不借助 $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 的步数平均下面显著的比快速幂小,但是由于用到了递归,因此最终它们的平均效率是一致的。
可以通过预处理小部分值达到快速的效果。

快速模乘法

1
2
3
4
5
6
7
8
9
// 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 就不同了。

1
2
3
4
5
6
7
8
9
// 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 定理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
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 !}$ 所以我们有:

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

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

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

乞丐版素数判断

1
2
3
4
5
6
7
8
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)$ 素数筛

1
2
3
4
5
6
7
8
9
10
11
12
13
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;
}
}
}

欧拉线性素数筛(正式使用版)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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();
}

欧拉线性素数筛(弃用)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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 概率判别法 和 大整数的最 大/小 因子分解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
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)$ 最小素因子预处理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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;
}

Mobius function

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
// 乞丐版
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;
}

Mobius function 前缀和

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
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;
}
// Mobius 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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
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 定理:

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

$\pi(x)$ 函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
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)$ 函数计算的另一种做法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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$ 是 模 $p$ 的原根。对于 $p^a$ 可以用 $p$ 的原根简单构造,而 $2p^a$ 的原根为 $p^a$ 的原根与 $p^a$ 的原根和 $p^a$的和中奇数者。(证明见 P150《数论基础》潘承洞)

求所有原根见

数论函数的 Dirichlet 乘积

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

中国剩余定理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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;
}

离散对数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
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;
}

模素数开根号

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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)$ 算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
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;
}

自然数方幂和精确版

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#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}]$,从而

上式等价于

参考:Soulist

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
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 nft(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);
nft(a);
nft(b);
for (int i = 0; i != sz; ++i) a[i] = a[i] * b[i] % M;
nft(a, 1);
a.resize(tot);
}

siyuan 的 FWT 模板 (弃用)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
#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 积 中的 Mobius 变换(本质就是带条件的求和)这里的变换也就是的各个位上的带条件求和。
例题:洛谷 P6097, 参考:yijan

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
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,也会被填满。

矩阵类

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
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 消元法浮点版

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#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 消元法有限域版

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
// 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;
}

线性规划单纯形法例题

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#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 P1993spfa 做法Bellman-Ford 做法

变式:$\frac{x_i}{x_j} \leq c_k$(取 log 即可)

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#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,比如通常是两个,也可以不用多项式带模除法来搞,只需 NFT 就可以优化到 $O(m \log m \log n + m^2)$ (暂时不知道如何去掉 $m^2$)

博弈

数据结构

并查集(Disjoint Set Union)

1
2
3
4
5
6
7
8
9
10
11
// 初始情形: 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;
}

离散化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
// 返回的是离散化之后的数组值对应的原始值
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;
}

逆序数

  1. 直接求 $O(n^2)$ 没啥好写的。
  2. 把原数组每个位置进行编号,排序,然后每次把最大的数的编号丢进树状数组中,丢进去先看这个编号前面有多少个数,累加一下就可以了,$O(n^2)$,结合下面树状数组的知识还是很简单的。
  3. 带离散化的树状数组(就是如果元素的数值特别大,树状数组内存就不够了,所以需要离散化一下)
  4. 归并的求(不会也不想搞 0.0)
  5. 逐位处理(代码如下)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// 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 到某个位置的最小值)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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 开始)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
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) 显然

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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 版本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
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 的模板,他的是左闭右开的。

线段树 (弃用版)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#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 已经死了)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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 $。

三分法简单版

1
2
3
4
5
6
7
8
9
// 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] 同理。也就是说我们有策略

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

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

三分法黄金分割版(可用版))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
// 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 的子列的末位的最小值。
这个子序列还能给附加限制(比如都要大于某个值)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
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 的位置

背包

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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 表示。 堆只支持三个操作:

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

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

单调队列:解决滑动窗口问题(固定长度内的最值问题)

知乎 Pecco 讲的很好(建议直接去看它的讲解):

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
// 求每个长度为 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]$ 的最小值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
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;
}
};

几何

二维凸包正式版

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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;
}

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
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]。 最后要清空树状数组(注意清空的时候不能直接清空,而是根据更新的命令,反向一次命令来清空,否则一直开树状数组耗时的),还有就是要去重贼麻烦,还是弃用吧。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
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})$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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 $

1
2
3
4
5
6
7
8
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)$

1
2
3
4
5
6
7
8
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$,所以用上述公式计算即可。

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

1
2
3
4
5
6
7
8
9
10
11
12
13
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;
}

递归算法复杂度分析

下图取自算法导论

complexAnalysis

To Be Continue

日常表白 zly