卡常 Tips

此页面的 Telegram 频道。

科普一些 O2 优化背景下的常数优化方法，以及一些底层优化的基础知识。

ST 表 (Sparse Table)

2023-11-02

• ST 表的数组，哪维放前面更快？
• $\lfloor\log_2(r-l+1)\rfloor$ 怎么算最快？

对于一个序列 $A$，定义 $t_{k,i}=\min A_{i-2^k+1\sim i}$，注意区间是 $[i-2^k+1,i]$ 而非 $[i,i+2^k-1]$，这样定义的好处是代码短一点，通过下面的式子预处理：

$$t_{k+1,i}=\min(t_{k,i-2^k},t_{k,i})$$

这里有两个问题：

• $t$ 数组应该把哪一维放前面？
• 预处理的两个 for 应该把哪个放前面？

最优的组合

for (int k = 0; k < log_n; k++)
    for (int i = 2 << k; i <= n; i++)
        t[k + 1][i] = min(t[k][i - (1 << k)], t[k][i]);

解释 (访问连续性, 指令级并行)

这里有两个知识，访问连续性 和 指令级并行。

• 访问连续性
• 指令级并行

当数组下标连续递增或递减时，访问的平均时间会很小，远小于随机访问。注意，在多个数组上连续访问也是有效的，数组的数量不超过 $256$ 应该都是可以接受的（有个例外是 4K Aliasing，它不重要）。

int sum = 0;
for (int i = 1; i <= n; i++) sum += a[i] + b[i] + c[i];

当简单循环中的每条语句的 互相独立 时，可以获得一种加速叫 指令级并行，多条机器指令可以同时执行。

互相独立是指执行顺序不重要，并且要修改的变量没有一样的。

// 顺序很重要
for (int i = 1; i <= n; i++) a[i] += a[i - 1];

// 都要修改 sum
int sum = 0;
for (int i = 1; i <= n; i++) sum += a[i];

// 互相独立
for (int i = 1; i <= n; i++) a[i] += val;

指令级并行 其实无处不在，只是循环中的语句互相独立的时候，更容易同时执行。

不难看出，最优组合同时具备 访问连续性 和 指令级并行。值得一提的是，下面这种组合也具备 访问连续性。

// 为避免下标为负数，加上一个常量 OFFSET
for (int i = 1; i <= n; i++)
    for (int k = 0; k < log_n; k++)
        t[k + 1][i + OFFSET] = min(t[k][i - (1 << k) + OFFSET], t[k][i + OFFSET]);

想一想为什么？

对于每一个 $i$，访问的值有

t[0][i] t[1][i] ... t[log_n][i]
t[0][i-1] t[1][i-2] ... t[log_n][i-(1<<log_n)]

每一个都是随着 $i$ 连续访问的，常数个指针同时连续移动也是一种 访问连续性。

询问的时候有两个问题：

• 数组哪一维放前面会比较快。
• $\lfloor\log_2(r-l+1)\rfloor$ 怎么算最快。

仍然是 $\log n$ 的一维放前面。

解释 (访问分布, 缓存)

不管哪一维放前面，数组都是随机访问的。

但是 $\log(r-l+1)$ 通常不是均匀分布的，不管是随机数据还是大部分构造数据，$\log(r-l+1)$ 倾向于集中于某个值，缓存 让访问得以加速。

缓存

缓存是一个缩小版的内存，比起内存有更高的访问速度，但可用空间也更小，访问频率最高的数据会自动留在缓存中。

因此当某些数据的访问频率较高时，缓存就会起作用，让平均访问时间变短。

当然，如果特意让 $\log_2(r-l+1)$ 在 $[0,\log_2 n]$ 均匀生成，哪一维放前面速度就一样了。

最后是计算 $\lfloor\log_2(r-l+1)\rfloor$ 的问题，最常见的方法是预处理一个 log_2 数组，更好的方法是使用 std::__lg(r - l + 1) 或 31 - __builtin_clz(r - l + 1)。

最快的方法

31 ^ __builtin_clz(r - l + 1)

当 $x$ 在 $[1,10^6]$ 均匀随机时，速度比较如下：

方法	时间
31 ^ __builtin_clz(x)	$1.00$
__lg(x)	$1.40$
log_2[x]	$8\sim9$

解释

• log_2[x] 需要访问较大数组，所以速度慢。
• __lg(x) 等价于 31 - __builtin_clz(x)，而 __builtin_clz(x) 会编译为 31 ^ __builtin_ia32_bsrsi(x)。
• 31 ^ __builtin_clz(x) 编译后两个异或会抵消，结果为 __builtin_ia32_bsrsi(x)，比前一种方法少了一次减法和异或。

set 和 map

2023-11-03

std::set 最常用的函数有：

• 用于修改的 insert(x), erase(x), erase(it)（指定迭代器）。
• 用于查询的 find(x), lower_bound(x), upper_bound(x)。
• 迭代器 it++, it--。

想减少 set 的时间开销，一般有两种方法：

利用 insert(x) 和 erase(it) 的返回值，减少函数调用。

• insert(x)
• erase(it)

insert(x) 的返回值

insert 函数的返回值类型为 pair<iterator, bool>，bool 代表元素是否插入成功，iterator 则是一个迭代器，指向所插入元素，或者指向原本就在容器中的元素。

实现一个函数，如果 $x$ 不在 set 中就插入 $x$，否则删除 $x$。

void foobar(set<int>& s, int x) {
    auto [it, t] = s.insert(x);
    if (!t) s.erase(it);
}

erase(it) 的返回值

erase(it) 函数的返回值为迭代器 it 的后继，即 ++it。

常见用法为用 it = s.erase(it) 替代 s.erase(it++)。

实现一个函数，将 set 中属于 $[l,r]$ 的元素全部删除，并统计删除的个数。

int foobar(set<int>& s, int l, int r) {
    auto it = s.lower_bound(l);
    int res = 0;
    while (it != s.end() && *it <= r)
        it = s.erase(it), res++;
    return res;
}

多用 insert(it, x), erase(it), it++, it-- 等高效函数。

insert(it, x) 的用法

当我们想插入 x，并且知道插入的位置恰好在迭代器 pos 之前，就可以使用这个函数，复杂度 $O(1)$，返回值为插入元素的迭代器。

如果 x 没有插入在迭代器 pos 之前，它就会重新调用 insert(x)，所以复杂度会退化到 $O(\log n)$。

给定一个数组 A[0..n-1]，把它插入到 set 中。

set<int> s;
sort(A, A + n);
int m = unique(A, A + n) - A;
for (int i = 0; i < m; i++) s.insert(s.end(), A[i]);

同时插入 $x$ 和 $x-1$。

void foobar(set<int>& s, int x) {
    auto it = s.insert(x).first;
    s.insert(it, x - 1);
}

erase(it) 的妙用

erase(it) 函数复杂度是 $O(1)$。如果我们记录了每次插入产生的迭代器，就可以 $O(1)$ 删除之前插入的元素。

实现插入、删除和 lower_bound 操作，保证操作合法，元素值域为 $[0,n)$。

set<int> s;
set<int>::iterator its[n];
void insert(int x) {
    its[x] = s.insert(x);
}
void erase(int x) {
    s.erase(its[x]);
}
void lower_bound(int x) {
    return *s.lower_bound(x);
}

it++, it-- 的复杂度

理论上它们在随机数据下的复杂度是 $O(1)$ 的，但是大量测试的拟合复杂度为 $O(\log n)$，我猜测可能和缓存有关。尽管如此，它的速度仍然是 find, lower_bound 等函数的 $10$ 倍左右。

一些便利的函数

• prev(it) 返回迭代器 it 的前驱，返回值等于 --it，但不会修改 it。next(it) 返回后继。
• map<int, int> 和 set<pair<int, int>> 可以使用 emplace(x, y) 来插入元素对。知道位置恰好在 it 之前的时候，使用 emplace_hint(it, x, y) 可以 $O(1)$ 插入。

map<int, int> m;
set<pair<int, int>> s;
m.emplace(1, 2); // 等价于 m.insert(make_pair(1, 2))
// map::insert 不会覆盖已有的值，不同于 [] 赋值
s.emplace(1, 2); // 等价于 s.insert(make_pair(1, 2))
m.emplace_hint(m.end(), 2, 3);
s.emplace_hint(m.end(), 2, 3);

初探循环展开

11-06

本节在 O2 优化下讨论，因为 O3 会自动循环展开，编译器版本为 GCC 11.4。

问题引入

下面是一段简单数组求和代码，总共求了 $10^6$ 次和。

unsigned a[5000], sum = 0;
iota(a, a + 5000, 0);
for (int k = 0; k < 1'000'000; k++)
    for (int i = 0; i < 5000; i++)
        sum += a[i];
assert(sum); // 防止编译器把这段代码删了

在我电脑上的运行时间大概为 1.55s。

接下来进行一个简单地循环展开。

// 前面省略
    for (int i = 0; i < 5000; i += 4) {
        sum += a[i + 0];
        sum += a[i + 1];
        sum += a[i + 2];
        sum += a[i + 3];
    }
// 后面省略

运行时间大概为 1.16s，似乎并没有获得显著的加速。

换一个写法：

// 前面省略
    for (int i = 0; i < 5000; i += 4)
        sum += a[i] + a[i + 1] + a[i + 2] + a[i + 3];
// 后面省略

运行时间大概为 0.65s，直接比第一种循环展开快了 $1$ 倍。

循环展开的原理

从上面的例子可以看出，运行时间并不是完全由运算量决定的。

一个极端的例子

$O(n^3)$ 矩阵乘法的运算量是固定的，但是 python 的矩阵乘法真的很快。

python

import time
import numpy as np
a = np.random.rand(5000, 5000)
b = np.random.rand(5000, 5000)

start = time.time()
c = np.matmul(a, b)
end = time.time()
print(end - start)

$5000\times 5000$ 的 double 矩阵乘法用时 1.1s。

决定 CPU 实际效率的因素主要有三个：

• 指令流水线 (instruction pipelining)。
• 缓存 (cache)。
• SIMD (Single Instruction Multiple Data)

指令流水线 是实现 指令级并行 的技术，也是循环展开可以带来加速的重要原因。

执行任何一条机器指令，需要众多步骤，不过可以简单地分为 CPU 前端的流程和后端的流程。CPU 前端在内存中一次性读取多条机器代码，然后加入 指令流水线 的末尾，后端同时处理着 指令流水线 中的每个流程，尽可能使众多执行单元保持忙碌。

事实上流水线不可能完全通畅，流水线会因为三种原因阻塞：

• 同一类型指令太多，相应的执行单元不够用。
• 当前指令需要得到上一条指令的结果才能开始执行。
• 分支预测失败。

我们来解释一下这个求和代码在 CPU 中的执行细节。

for (int i = 0; i < n; i++) sum += a[i];

它编译后的汇编等效于这段 C++ 代码：

void *rax = a, *rcx = a + n;
unsigned edx = 0;
loop:
    edx += *(unsigned*)rax;
    rax += 4;
    if (rax != rcx)
        goto loop;

实际的汇编是这样的：

loop:
    add edx, DWORD PTR [rax]
    add rax, 4
    cmp rax, rcx
    jne loop

另外，机器代码和汇编是一一对应的，所以机器代码也由这些语句构成。CPU 前端负责解析机器代码和完成循环的跳转 (jne)，CPU 后端则执行两条 add 和一条 cmp，由于 CPU 的执行单元足够多，所以三条指令会占用不同的执行单元。

第一次的循环展开变快很好解释，因为指令变少了，加速不显著是因为这些减少的指令原本就不和 add edx, DWORD PTR [rax] 共用执行单元。

第二次的循环展开变快的原因就很有意思了，指令数完全没有变少，但是加速很显著。第一次的循环展开只是把 sum += a[i + ?] 重复了四遍，它们构成了一条依赖链，因为 sum += a[i + 1] 必须要等上一条语句把 sum 算出来后才能执行，所以流水线会阻塞。而 a[i] + a[i + 1] + a[i + 2] + a[i + 3] 不存在这样严重的依赖链，对于每个 $i$，a[i] + a[i + 1] + a[i + 2] + a[i + 3] 都可以并发地计算，虽然 sum += 的部分确实构成了依赖链，但是长度短得多。

最后给出一种推荐的写法：

unsigned s0 = 0, s1 = 0, s2 = 0, s3 = 0;
for (int i = 0; i < 5000; i += 4) {
    s0 += a[i + 0];
    s1 += a[i + 1];
    s2 += a[i + 2];
    s3 += a[i + 3];
}
sum += s0 + s1 + s2 + s3;

运行时间大概为 0.62s，比第二次的循环展开有微小的进步。

拓展

除了求和之外，还有三个典型的例子：

整体加

unsigned a[5000] = {};
for (int k = 0; k < 1'000'000; k++)
    for (int i = 0; i < 5000; i++)
        a[i]++;
assert(*max_element(a, a + 5000));

循环展开后有显著加速，因为 a[i]++ 可以并发计算。

前缀和

下面两段代码速度差异很大。

• 第一种
• 第二种

unsigned a[5000];
std::iota(a, a + 5000, 0);
for (int k = 0; k < 1'000'000; k++)
    for (int i = 1; i < 5000; i++)
        a[i] += a[i - 1];
assert(*max_element(a, a + 5000));

unsigned a[5000];
std::iota(a, a + 5000, 0);
for (int k = 0; k < 1'000'000; k++)
    for (int i = 0; i + 1 < 5000; i++)
        a[i + 1] += a[i];
assert(*max_element(a, a + 5000));

哪个快建议试一下，我觉得比较慢的那种是编译器锅了（GCC 12 修了），没优化到。

然后比较快的那种编译器会优化成这个样子：

int s = 0;
for (int i = 0; i < 5000; i++)
    s += a[i], a[i] = s;

循环展开后有较小加速。对于循环体很小的循环，展开后总是会有加速，因为指令数会减半。s += a[i] 构成了依赖链，没法并行，但 a[i] = s 能在一定程度上并行。

小型矩阵乘法

void mul(unsigned a[8][8], unsigned b[8][8], unsigned c[8][8]) {
    for (int i = 0; i < 8; i++)
        for (int j = 0; j < 8; j++) {
            uint64_t sum = 0;
            for (int k = 0; k < 8; k++)
                sum += (uint64_t)a[i][k] * b[k][j];
            c[i][j] = sum % 998244353;
        }
}

把最内层循环打开后效果很好，这样乘法就可以并行。

使用 pragma 自动展开

全局自动展开就不说了，我觉得那个比较鸡肋。

推荐的做法是，在你想要展开的 for 前面加一句 #pragma GCC unroll n，其中 n 是一个数字，表示打开循环的步长。

以矩阵乘法为例：

void mul(unsigned a[8][8], unsigned b[8][8], unsigned c[8][8]) {
    for (int i = 0; i < 8; i++)
        for (int j = 0; j < 8; j++) {
            uint64_t sum = 0;
            // 完全展开
            #pragma GCC unroll 8
            for (int k = 0; k < 8; k++)
                sum += (uint64_t)a[i][k] * b[k][j];
            c[i][j] = sum % 998244353;
        }
}

使用 C++ template 自动展开

由于 OI 不能用 #pragma，所以再给一种方法，大概定义一个这样的函数：

template <int L, int R, int M = (L + R) / 2>
__attribute__((always_inline)) inline int foobar() {
    if (R - L == 1) return get_val(L);
    return foobar<L, M>() + foobar<M, R>();
}

__attribute__((always_inline)) inline 表示强制内联，再配合 template 的递归就可以达到展开的效果。

还是以矩阵乘法为例：

template <int L, int R, int M = (L + R) / 2>
__attribute__((always_inline)) inline uint64_t foobar(unsigned a[8][8], unsigned b[8][8], int i, int j) {
    if (R - L == 1) return (uint64_t)a[i][L] * b[L][j];
    return foobar<L, M>(a, b, i, j) + foobar<M, R>(a, b, i, j);
}
void mul(unsigned a[8][8], unsigned b[8][8], unsigned c[8][8]) {
    for (int i = 0; i < 8; i++)
        for (int j = 0; j < 8; j++)
            c[i][j] = foobar<0, 8>(a, b, i, j) % 998244353;
}

当然这样不太实用，每个循环都需要单独写个函数，可以添加一个匿名函数作为参数。

template <int L, int R, int M = (L + R) / 2, class F>
__attribute__((always_inline)) inline void foobar(int i, F lambda) {
    if (R - L == 1) lambda(i + L);
    else foobar<L, M>(i, lambda), foobar<M, R>(i, lambda);
}

不太需要担心匿名函数的调用开销，因为匿名函数比普通函数更倾向于 inline，使用方法如下：

for (int i = 0; i < 5000; i += 4)
    foobar<0, 4>(i, [](int i) { a[i]++; });

int sum = 0;
for (int i = 0; i < 5000; i += 4)
    foobar<0, 4>(i, [&](int i) { sum += a[i]; });

循环的合并

11-10

如果需要对数组 a 做两遍前缀和，最直接的写法是这样：

int sum0 = 0;
for (int i = 0; i < 5000; i++)
    sum0 += a[i], a[i] = sum0;
int sum1 = 0;
for (int i = 0; i < 5000; i++)
    sum1 += a[i], a[i] = sum1;

完整测试代码

#include <bits/stdc++.h>

using namespace std;

int main() {
    unsigned a[5000];
    iota(a, a + 5000, 0);
    for (int k = 0; k < 1'000'000; k++) {
        int sum0 = 0;
        for (int i = 0; i < 5000; i++)
            sum0 += a[i], a[i] = sum0;
        int sum1 = 0;
        for (int i = 0; i < 5000; i++)
            sum1 += a[i], a[i] = sum1;
    }
    assert(*max_element(a, a + 5000));
}

运行时间大概 3.5s。

如果把两个循环合并起来：

int sum0 = 0, sum1 = 0;
for (int i = 0; i < 5000; i++)
    sum0 += a[i], sum1 += sum0, a[i] = sum1;

那么运行时间会减少到 2s。

合并循环最容易看到的好处是循环变量 i 的开销减半了，跳转的数量也减半了，不过这不是最主要的好处，真要减少循环开销还得用循环展开。以下两点才是合并循环的主要好处：

• 内存读写减少了，当我们使用 a[i] 的值时，就需要读取内存，给 a[i] 赋值时，就需要写入内存，即使访问连续，内存读写也是有固定开销的。注意，这并不是说 a[i] + a[i] 会读取内存两次，O2 优化会智能地减少内存读写的次数。
• 缓解依赖关系。由于前缀和每次操作都依赖于上一次操作，所以指令级并行程度较低，合并循环后，仍然存在依赖关系，但是没那么严重了，更容易指令级并行。

我们把循环合并前和合并后的代码都循环展开，对比一下速度：

	合并前	合并后
不循环展开	3.5s	2s
循环展开	2.45s	1.36s

不合并循环 + 循环展开的代码

#include <bits/stdc++.h>

using namespace std;

int main() {
    unsigned a[5000];
    iota(a, a + 5000, 0);
    for (int k = 0; k < 1'000'000; k++) {
        int sum0 = 0;
        #pragma GCC unroll 8
        for (int i = 0; i < 5000; i++)
            sum0 += a[i], a[i] = sum0;
        int sum1 = 0;
        #pragma GCC unroll 8
        for (int i = 0; i < 5000; i++)
            sum1 += a[i], a[i] = sum1;
    }
    assert(*max_element(a, a + 5000));
}

合并循环 + 循环展开的代码

#include <bits/stdc++.h>

using namespace std;

int main() {
    unsigned a[5000];
    iota(a, a + 5000, 0);
    for (int k = 0; k < 1'000'000; k++) {
        int sum0 = 0, sum1 = 0;
        #pragma GCC unroll 8
        for (int i = 0; i < 5000; i++)
            sum0 += a[i], sum1 += sum0, a[i] = sum1;
    }
    assert(*max_element(a, a + 5000));
}

当数组很大的时候，内存读写很慢，即使连续访问也会变慢，这时合并循环更重要了：

	合并前	合并后
$5000$ 数组，迭代 $10^6$ 次（循环展开）	2.45s	1.36s
$5\times 10^6$ 数组，迭代 $1000$ 次（循环展开）	4.3s	2.3s

合并循环在 OI 中的应用

主要应用在 DP 和组合计数上，举两个例子：

$dp_{i,j}=3dp_{i-1,j}+2dp_{i-1,j-1}+dp_{i-1,j-2}$

慢的写法：

for (int j = 0; j <= n; j++) dp[i][j] = 3LL * dp[i - 1][j] % mod;
for (int j = 1; j <= n; j++) dp[i][j] = (dp[i][j] + 2LL * dp[i - 1][j - 1]) % mod;
for (int j = 2; j <= n; j++) dp[i][j] = (dp[i][j] + 1LL * dp[i - 1][j - 2]) % mod;

快的写法是把转移写在一个循环里，单独处理边界：

for (int j = 2; j <= n; j++) dp[i][j] = (3LL * dp[i - 1][j] + 2 * dp[i - 1][j - 1] + dp[i - 1][j - 2]) % mod;
dp[i][0] = 3LL * dp[i - 1][0] % mod;
dp[i][1] = (3LL * dp[i - 1][1] + 2 * dp[i - 1][0]) % mod;

在取模的题中，合并循环还有个好处是减少取模次数。

高精度 A + B - C。

虽然代码实现本身并不复杂，但高精度通常会被封装为 struct 进行使用。为了能够合并循环，我们不能使用 $\pm$ 运算符，这确实有些遗憾，于是就需要自动化合并循环了。

合并循环自动化

这东西代码很长，OI 不实用，但它可以让高精度类的 $\pm$ 表达式自动合并循环。

如果想实现这样一个高精度类，请参考 C++ expression templates。

循环的分裂

11-11

这节内容很简单，就是当循环里有 if 或 min,max 而且它们开销比较大的时候，想办法把它们去掉。

预处理阶乘和阶乘逆元。

慢的写法：

for (int i = 0; i <= n; i++)
    fac[i] = i == 0 ? 1 : 1LL * fac[i - 1] * i % mod;
for (int i = n + 1; i; i--)
    ifac[i - 1] = i > n ? Pow(fac[n], mod - 2) : 1LL * ifac[i] * i % mod;

快的写法是单独初始化 fac[0] 和 ifac[n]，循环中就不用三目了。

cpp

fac[0] = 1;
for (int i = 0; i <= n; i++) fac[i] = 1LL * fac[i - 1] * i % mod;
ifac[n] = Pow(fac[n], mod - 2);
for (int i = n + 1; i; i--) ifac[i - 1] = 1LL * ifac[i] * i % mod;

 

计算 $S=\sum_{i=0}^nA_i\times\min(i,n-i)$

把 $\min(i,n-i)$ 拆开：

int S = 0;
for (int i = 1; i <= m / 2; i++) S += A[i] * i;
for (int i = m / 2 + 1; i < m; i++) S += A[i] * (m - i);

这样还不够好，可以把乘法去掉，注意到 $nA_1+(n-1)A_2+\cdots+2A_{n-1}+A_n$ 等于 $A$ 序列的前缀和之和，利用这个公式可以去掉乘法，最后来个循环展开：

int S = 0;
int sum0 = 0;
#pragma GCC unroll 8
for (int i = m / 2; i; i--) S += sum0 += A[i];
int sum1 = 0;
#pragma GCC unroll 8
for (int i = m / 2 + 1; i < m; i++) S += sum1 += A[i];

常见无效卡常

11-11

会持续更新。在 O2 优化下，很多事情编译器自己会做，有些“优化”其实没有效果。

快读中 x = x * 10 + c - 48 写成 x = (x << 3) + (x << 1) + (c ^ 48)。

x * 10 写成 (x << 3) + (x << 1) 真没什么意义，编译出来没区别，另外 x * 10 的最优实现并不是两个左移一次加法，而是一次 lea 一次加法，至于 c - 48 改成 c ^ 48 有的电脑上变快了有的变慢了，目前还没法解释。

模算术优化 I

11-1?