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

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

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

通用代码块

佛祖保佑,永无 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 下无法设置可惜)

产生 log 的几个原因

  1. 二分,三分
  2. $1 + \frac{1}{2} + \cdots \frac{1}{n} ~ \log n$
  3. 树状数组,线段树
  4. 重链剖分
  5. 倍增

产生根号的几个原因

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

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

在前面加 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 高端压位卡常

位运算的关系

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

Meet in Middle(拆半搜索法)

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

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

Small to large(把小的合并到大的里面去)

例子:并查集(dus),map 的合并,树上启发式合并(dus on tree),重链剖分。

例题:600E题解

倍增思想

例子:RMQ,LCA。

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 也没啥用,仅仅记录
int gcd(int a, int b) {
if (!a || !b) return a | b;
unsigned shift = __builtin_ctz(a | b);
a >>= __builtin_ctz(a);
do {
b >>= __builtin_ctz(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
10
template<typename T,typename U>
T powMod(T x,U n,T p){
T r(1);
while (n) {
if (n&1) r = r * x % p;
n >>= 1; x = x * x % p;
}
return r;
}
// 可以用for循环写的更短一点,但没必要

快速模加法下的快速模乘法(弃用,利用 __int128 即可,类似的还有快速模加法也弃用了)

1
2
3
4
5
6
7
8
LL powMod(LL x, LL n, LL p){
LL r(1);
while (n){
if (n&1) r = mulmod(r, x, p);
n >>= 1; x = mulmod(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;
}

常用组合数公式

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

组合意义理解:$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;

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;
}
}
}

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

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();
}

$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
std::mt19937 rnd(std::chrono::steady_clock::now().time_since_epoch().count());
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 = mulMod(x, x, n);
if (x == n - 1) return false;
}
return true;
}
bool rabin(LL n) {
if (n < 2) return false;
if (n == 2) return true;
if (!(n & 1)) return false;
LL m = n - 1;
int t = 0, cnt = 33;
while (!(m & 1)) {
++t;
m >>= 1;
}
while (cnt--) {
LL a = rnd() % (n - 1) + 1;
if (Witness(a, n, m, t)) return false;
}
return true;
}

大整数的最小素因子

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
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;
}
bool rabin(LL n) {
if (n < 2) return false;
if (n == 2) return true;
if (!(n & 1)) return false;
LL m = n - 1;
int t = 0, cnt = 33;
while (!(m & 1)) {
++t;
m >>= 1;
}
while (cnt--) {
LL a = std::rand() % (n - 1) + 1;
if (witness(a, n, m, t)) return false;
}
return true;
}

大整数的最小素因子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
LL pollardrho(LL n) {
LL x = std::rand() % (n - 1) + 1;
LL y = x, i = 1, k = 2, c = std::rand() % (n - 1) + 1;
while (true) {
x = (__int128(x) * x + c) % n;
LL d = std::__gcd(y - x + n, n);
if (d > 1) return d;
if (++i == k) {
y = x;
if (k > n) return n;
k <<= 1;
}
}
}
LL findp(LL n) {
if (rabin(n)) return n;
LL p = n;
while (p == n) p = pollardrho(n);
return std::min(findp(p), findp(n/p));
}

$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;
}

$\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 真香!

矩阵类

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());
}

求$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
12
13
14
15
16
17
// 初始情形: p[i] = i
int find(int x) {
int ans = x;
while (ans != p[ans]) ans = p[ans];
while (x != ans) {
int t = p[x];
p[x] = ans;
x = t;
}
return ans;
}
void merge(int i, int j) {
// p[find(j)] = p[find(i)];
// In general we should write below
// 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;
}
};

树状数组(区间更新,区间求和,编号从 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$)

和与最大值的线段树模板(如果单纯求和,可以用树状数组)

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
10
11
// for a given covex function f: (f(a) + f(b)) / 2 >= f((a + b) / 2)
int l = 1, r = 1e9;
while (l + 2 < r) {
int lm = (2ll * l + r) / 3, rm = (l + 2ll * r + 2) / 3;
if (f(lm) < f(rm)) r = rm;
else l = lm;
}
while (l < r) {
if (f(l) < f(r)) --r;
else ++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$,但是会很不优雅,所以还是算了。

最长(严格)递增子序列

下面数组 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 pdd = std::pair<double, double>;
bool crossLeft(const pii &op, const pii &sp, const pii &ep) {
return (sp.first - op.first) * (ep.second - op.second)
< (sp.second - op.second) * (ep.first - op.first);
}
std::vector<pdd> convexHull(std::vector<pdd> p) {
std::sort(p.begin(), p.end());
p.erase(std::unique(p.begin(), p.end()), p.end());
int n = p.size();
std::vector<pdd> 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
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
const int N = 501;
int vc[N];
struct Node {
int x, y, id;
bool operator != (const Node & A) const {
return x != A.x || y != A.y;
}
bool operator < (const Node & A) const {
if (y == A.y) return x < A.x;
return y < A.y;
}
}p[N], q[N];
bool crossLeft(const Node & op,
const Node & sp,
const Node & ep) {
return (sp.x - op.x) * (ep.y - op.y) > (sp.y - op.y) * (ep.x - op.x);
}
int graham(int n) {
sort(p, p + n);
if (n == 0) return 0;
q[0] = p[0];
if (n == 1) return 1;
q[1] = p[1];
if (n == 2) return 2;
q[2] = p[2];
int top = 1;
for (int i = 2; i != n; ++i) {
while (top && crossLeft(q[top], p[i], q[top - 1])) --top;
q[++top] = p[i];
}
int len = top;
q[++top] = p[n - 2];
for (int i = n - 3; i != -1; --i) {
while (top != len && crossLeft(q[top], p[i], q[top - 1])) --top;
q[++top] = p[i];
}
return top;
}

几类根号算法

首先注意到 $\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