数学与计算

数论

数论主要研究整数的性质。

筛法

筛法是一种在数论中通过系统排除不满足条件的数来识别符合条件的数集合的方法，使用筛法可以快速筛选出指定区间的素数，也可以拓展到求积性函数。

埃拉托斯特尼筛

如果我们从小到大考虑每个数，将每个遍历到的素数的倍数（大于等于自己的平方，即该素数作为最小质因数）记为合数，那么运行结束的时候没有被标记的数就是素数，注意 0 和 1 要特殊处理。时间复杂度 O(nlog log n)。

注意到筛法求素数的同时也得到了每个数的最小质因子，因此筛法可以在 O(nlog n) 的时间复杂度内对 [1, n] 区间所有数分解质因数，或在 O(nk) 的时间复杂度内将 [1, n] 区间的所有数分解出至多 k 个质因数，同时判定每个数是否分解完全。

基于埃拉托斯特尼筛法可以还可以实现区间筛，由于每个合数只会被其最小质因数标记，且最小质因数不大于 $\sqrt{n}$，因此只需筛至 $\sqrt{n}$ 即可，在筛的过程中标记所求区间素数即可，复杂度 $O(\sqrt{n}\log{\log{n}})$ 。

以下代码实现了埃氏筛，其中 pr 为素数列表， f[i] 代表 i 的最小质因数，判断质数的条件为 f[i] =  = i 。

1
2
3
4
5
6
7
8
9
10
vector<int> f(mxv + 1), pr;  
void init() {  
    for (i64 i = 2; i <= mxv; i++) {  
        if (f[i]) continue;  
        f[i] = i, pr.push_back(i);  
        for (i64 j = i * i; j <= mxv; j += i) {  
            if (!f[j]) f[j] = i;  
        }  
    }  
}

欧拉筛/线性筛

欧拉筛也称线性筛，每个合数只会被其最小质因数标记一次，时间复杂度 O(n) 。

由于合数 ni = i × pr[j] ，其中 pr[j] 为 i × x 的最小质因数，只有在枚举到数 i 时被更新，而 i 可能大于 $\sqrt{n}$ ，因此不可以筛至平方根。

以下代码实现了欧拉筛，其中 pr 为素数列表， f[i] 代表 i 的最小质因数，判断质数的条件为 f[i] =  = i 。

1
2
3
4
5
6
7
8
9
10
vector<int> f(mxv + 1), pr;  
void init() {  
    for (int i = 2; i <= mxv; i++) {  
        if (!f[i]) f[i] = i, pr.push_back(i);  
        for (int j = 0, ni; j < pr.size() && (ni = i * pr[j]) <= mxv; j++) {  
            f[ni] = pr[j];  
            if (i % pr[j] == 0) break;  
        }  
    }  
}

线性筛欧拉函数

对正整数 n，欧拉函数 φ(n) 是小于等于 n 的正整数中与 n 互质的数的数目。

  1. 如果 n = 1，则 φ(1) = 1 。因为 1 与任何数（包括自身）都构成互质关系。
  2. 如果 n 是质数，如果 n 是质数，则 φ(n) = n − 1。因为质数与小于它的每一个数，都构成互质关系。
  3. 如果 n 是质数的 k 次方，即 n = pk ，则 $\varphi(p^k) = p^k - p^{k-1} = p^k (1 - \frac{1}{p})$ 。
  • 这是因为一个数只要不包含质因数 p，就与 n 互质，而满足 k × p ≤ pk 的数一共有 pk − 1 个 。
  1. 如果 n 可以分解成两个互质的整数之积 p1 × p2 ，则 φ(n) = φ(p1p2) = φ(p1)φ(p2)
  • 这是因为当 a 与 p1 互质，b 与 p2 互质，c 与 p1p2 互质，则 c 与数对 (a, b) 是一一对应关系。
  1. 任何一个大于 1 的正整数，都可以写一系列质数的积。

由以上推导，可以得出欧拉函数的计算公式：

$$ \begin{equation} \begin{split} \varphi(n) &= \varphi(p_1^{k_1} p_2^{k_2} \ldots p_n^{k_n}) \\ &= \varphi(p_1^{k_1}) \varphi(p_2^{k_2}) \ldots \varphi(p_n^{k_n}) \\ &= p_1^{k_1} p_2^{k_2} \ldots p_n^{k_n} (1 - \frac{1}{p_1})(1 - \frac{1}{p_2}) \ldots (1 - \frac{1}{p_n}) \\ &= n(1 - \frac{1}{p_1})(1 - \frac{1}{p_2}) \ldots (1 - \frac{1}{p_n}) \end{split} \end{equation} $$

根据欧拉函数的定义，设 n 是一个合数，p1 是 n 的最小质因子，$n' = \frac{n}{p_1}$ ：

如果 n′ mod  p1 = 0，那么 n′ 包含了 n 的所有质因子。

$$ \varphi(n) = n \times \prod_{i=1}^s \frac{p_i - 1}{p_i} $$

$$ = p_1 \times n' \times \prod_{i=1}^s \frac{p_i - 1}{p_i} $$

 = p1 × φ(n′)

那如果 n′ mod  p1 ≠ 0 呢，这时 n′ 和 p1 是互质的，根据欧拉函数性质，我们有：

φ(n) = φ(p1) × φ(n′)

 = (p1 − 1) × φ(n′)

若 n 是质数，有 μ(n) = n − 1 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
vector<i64> phi(mxv + 1), f(mxv + 1), pr;  
void init() {  
    phi[1] = 1;  
    for (int i = 2; i <= mxv; i++) {  
        if (!f[i]) f[i] = i, phi[i] = i - 1, pr.push_back(i);  
        for (int j = 0, ni; j < pr.size() && (ni = i * pr[j]) <= mxv; j++) {  
            f[ni] = pr[j];  
            if (i % pr[j] == 0) {  
                phi[ni] = phi[i] * pr[j];  
                break;  
            }  
            phi[ni] = phi[i] * phi[pr[j]];  
        }  
    }  
}

线性筛莫比乌斯函数

莫比乌斯函数定义如下：

$$ \mu(n) = \begin{cases} 1 & n = 1 \\ 0 & \exists d > 1, d^{2} \mid n \\ (-1)^{\omega(n)} & \text{otherwise} \end{cases} $$

其中 ω(n) 表示 n 的本质不同质因子个数。

根据莫比乌斯函数的定义，设 n 是一个合数，p1 是 n 的最小质因子，$n' = \frac{n}{p_1}$，有：

$$ \mu(n) = \begin{cases} 0 & n' \mod p_1 = 0 \\ -\mu(n') & \text{otherwise} \end{cases} $$

若 n 是质数，有 μ(n) = −1 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
vector<i64> mu(mxv + 1), f(mxv + 1), pr;  
void init() {  
    mu[1] = 1;  
    for (int i = 2; i <= mxv; i++) {  
        if (!f[i]) f[i] = i, mu[i] = -1, pr.push_back(i);  
        for (int j = 0, ni; j < pr.size() && (ni = i * pr[j]) <= mxv; j++) {  
            f[ni] = pr[j];  
            if (i % pr[j] == 0) {  
                mu[ni] = 0;  
                break;  
            }  
            mu[ni] = -mu[i];  
        }  
    }  
}

数论分块

使用数论分块求 $\sum_{i=1}^{n}\lfloor\frac{n}{i}\rfloor$ ，将 $\lfloor\frac{n}{i}\rfloor$ 相同的数打包处理，复杂度 $O(\sqrt{n})$ 。

1
2
3
4
5
i64 ans = 0;  
for (i64 l = 1, r; l <= n; l = r + 1) {  
    r = n / (n / l);  
    ans += n / l * (r - l + 1);  
}

素数与质因数

Miller-Rabin 素性测试

进阶的素数判定方法，快速判断指定数字是否为质数，在 long long 范围内测试时间复杂度为 O(klog n) 。

可在底数列表中添加随机数以进行更多轮试验，常用于对高精度数进行测试，此时时间复杂度是 O(klog3n) ，其中 k 为试验轮数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// vector<int> base{2, 7, 61}; int范围内有效
vector<i64> base{2, 325, 9375, 28178, 450775, 9780504, 1795265022};

bool miller_rabin(i64 n) {
    if (n < 3 || n % 2 == 0) return n == 2;
    if (n % 3 == 0) return n == 3;
    i64 u = n - 1, t = 0;
    while (u % 2 == 0) u /= 2, ++t;
    for (auto a: base) {
        if (a % n >= -1 && a % n <= 1) continue;
        i64 v = qpow(a % n, u, n);
        if (v == 1) continue;
        i64 s;
        for (s = 0; s < t; ++s) {
            if (v == n - 1) break;
            v = qmul(v, v, n);
        }
        if (s == t) return false;
    }
    return true;
}

Pollard-Rho 分解因数

Pollard-Rho 算法是一种随机化算法，可以在 $O(\sqrt p)=O(n^{1/4})$ 的期望复杂度获得一个非平凡因子（非平凡因子不一定是素因子），其中 p 是 n 的最小素因子。

有概率分解失败，退出并返回 n 本身。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
i64 pollard_rho(i64 x) {
    i64 s = 0, t = 0;
    i64 c = (i64) rnd() % (x - 1) + 1;
    i64 val = 1, step = 0;
    for (int goal = 1;; goal *= 2, s = t, val = 1) {
        for (step = 1; step <= goal; ++step) {
            t = (qmul(t, t, x) + c) % x;
            val = qmul(val, abs(t - s), x);
            if (step % 127 == 0) {
                i64 d = gcd(val, x);
                if (d > 1) return d;
            }
        }
        i64 d = gcd(val, x);
        if (d > 1) return d;
    }
}

分解质因数

利用 Miller-Rabin 素性测试与 Pollard-Rho 分解质因数，期望复杂度 O(n1/4logn) 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
i64 maxf = 0;
auto fac = [&](auto self, i64 x) {
    if (x <= maxf || x < 2) return;
    if (miller_rabin(x)) {
        maxf = max(maxf, x);
        return;
    }
    i64 p = x;
    while (p >= x) p = pollard_rho(x);
    while (x % p == 0) x /= p;
    if (x < p) swap(x, p);
    self(self, x), self(self, p);
};
fac(fac, n);

最大公约数

Stein 算法

由普通的辗转相减法，有 gcd (a, b) = gcd (a − b, b) ，直接这样操作的最坏时间复杂度为 O(n)。

先计算质因数 2 的贡献，通过右移除 2 保证每次辗转相减前 a 和 b 均为奇数，则每次 a − b 一定为偶数，至少能除掉一个对答案没有贡献的 2 。

这个过程可以通过 buildin 位运算优化，最坏时间复杂度是 O(log n) ，但实际运行效率很高。

1
2
3
4
5
6
7
8
int qgcd(int a, int b) {  
    int az = __builtin_ctz(a), bz = __builtin_ctz(b), z = min(az, bz);  
    for (b >>= bz; a >>= az; az = __builtin_ctz(a)) {  
        if (a < b) swap(a, b);  
        a = a - b;  
    }  
    return b << z;  
}

扩展欧几里得算法

扩展欧几里得算法可以求得 ax + by = gcd (x, y) 的一组解，时间复杂度 O(log n) 。

1
2
3
4
5
6
7
8
9
int exgcd(int a, int b, int& x, int& y) {
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    int d = exgcd(b, a % b, x, y), tmp = x;
    x = y, y = tmp - a / b * y;
    return d;
}

求解线性同余方程

线性同余方程 ax ≡ c (mod  b) 可以改写为线性不定方程 ax + by = c ，求解步骤如下：

  1. 令 d = gcd (x, y)，如果 d 不能整除 c，则方程无解
  2. 化简方程：$a'=\frac{a}{d}$, $b'=\frac{b}{d}$, $c'=\frac{c}{d}$
  3. 使用扩展欧几里得算法求得 a′x + b′y = c′ 的一组特解 (x′, y′)
  4. 求得原方程的一组特解 x0 = c′x′, y0 = c′y′

原方程的所有解可以表示为 x = x0 + tb, y = y0 − ta ，其中 t 为任意整数。

1
2
3
4
5
6
7
bool liEu(int a, int b, int c, int& x, int& y) {
    int d = exgcd(a, b, x, y);
    if (c % d != 0) return false;
    int k = c / d;
    x *= k, y *= k;
    return true;
}

扩展中国剩余定理

扩展中国剩余定理通过扩展欧几里得算法合并同余方程，时间复杂度 O(log n) 。

1
2
3
4
5
6
7
8
9
10
bool excrt(i64 a1, i64 m1, i64 a2, i64 m2, i64 &a, i64 &m) {  
    i64 x, y, c = a2 - a1;  
    i64 d = exgcd(m1, m2, x, y);  
    if (c % d != 0) return false;  
    i64 t = m2 / d;  
    x = (i128) (x % t + t) % t * (c / d) % t;  
    a = a1 + x * m1, m = m1 / d * m2;  
    a = (a % m + m) % m;  
    return true;  
}

线性求逆元

以下方法可以在 O(n) 线性时间内求解出 1 − n 或任意 n 个数的逆元。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
vector<int> inv(n + 1);

void get_invs(int n) {
    inv[1] = 1;
    for (int i = 2; i <= n; ++i) {
        inv[i] = (p - p / i) * inv[p % i] % p;
    }
}

vector<i64> a(n + 1); // 所求任意n个数的列表
void get_invs(int n) {
    vector<i64> s(n + 1), sv(n + 1);
    s[0] = 1;
    for (int i = 1; i <= n; ++i) s[i] = s[i - 1] * a[i] % p;
    sv[n] = qpow(s[n], p - 2);
    for (int i = n; i >= 1; --i) sv[i - 1] = sv[i] * a[i] % p;
    for (int i = 1; i <= n; ++i) inv[i] = sv[i] * s[i - 1] % p;
}

min_25 筛

min_25 筛是一种用于快速计算积性函数前缀和的亚线性筛法，其在 1013 数据范围内的复杂度为常数很小的 $O(\frac{n^{3/4}}{\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
auto F = [](i64 x) { return x %= p, x * (x - 1 + p) % p; };  
auto S1 = [](i64 x) { return x %= p, x * (x + 1) / 2 % p; };  
auto S2 = [](i64 x) { return x %= p, x * (x + 1) / 2 % p * (2 * x + 1) % p * inv3 % p; };  
  
auto cal = [&](i64 n) {  
    i64 m = ceil(sqrt(n));  
  
    vector<i64> vis(m + 1), pr(1);  
    for (i64 i = 2; i <= m; i++) {  
        if (vis[i]) continue;  
        pr.push_back(i);  
        for (i64 j = i * i; j <= m; j += i) vis[j] = 1;  
    }  
  
    vector<i64> v(1), g1(1), g2(1), id1(m + 1), id2(m + 1);  
    auto id = [&](i64 x) -> i64 &{ return x < m ? id1[x] : id2[n / x]; };  
    for (i64 l = 1; l <= n; l = n / (n / l) + 1) {  
        i64 x = v.size();  
        v.push_back(n / l);  
        g1.push_back((S1(v[x]) - 1 + p) % p);  
        g2.push_back((S2(v[x]) - 1 + p) % p);  
        id(v[x]) = x;  
    }  
  
    for (i64 i = 1; i < pr.size(); i++) {  
        for (i64 j = 1; j < v.size() && pr[i] * pr[i] <= v[j]; j++) {  
            g1[j] = (g1[j] - pr[i] * (g1[id(v[j] / pr[i])] - g1[id(pr[i - 1])]) % p + p) % p;  
            g2[j] = (g2[j] - pr[i] * pr[i] % p * (g2[id(v[j] / pr[i])] - g2[id(pr[i - 1])]) % p + p) % p;  
        }  
    }  
  
    auto dfs = [&](auto &&self, i64 x, i64 y) -> i64 {  
        if (pr[y] >= x) return 0;  
        i64 res = (g2[id(x)] - g2[id(pr[y])] - g1[id(x)] + g1[id(pr[y])] + 2 * p) % p;  
        for (i64 i = y + 1; i < pr.size() && pr[i] * pr[i] <= x; i++) {  
            for (i64 j = pr[i]; j * pr[i] <= x; j *= pr[i]) {  
                res = (res + F(j) * self(self, x / j, i) + F(j * pr[i])) % p;  
            }  
        }  
        return res % p;  
    };  
  
    return (dfs(dfs, n, 0) + 1) % p;  
};

朴素 min_25 筛求素数计数函数

每个合数的最小质因数不大于 $\sqrt{n}$，因此只需枚举 $[1,\sqrt{n}]$ 的每个素数 i ，减去最小质因子为 i 的合数个数即可。

定义 v[i][j] 为区间 [1, j] 划去最小质因子小于等于 i 的合数后的素数个数，可以使用埃氏筛类似的方法转移，复杂度 $O(\sqrt{n}\log{\log{n}})$ 。

定义 p[i][j] 为区间 $[1,\frac{n}{j}]$ 划去最小质因子小于等于 i 的合数后数字个数，有 $p[i][\frac{n}{j}]$ 等于 v[i][j] ，$p[\sqrt{n}][1]$ 即为所求。

枚举素数 i 和所需更新的 p[i][j] ，则划去的数字个数可以表示为满足 $t*i≤\frac{n}{j}$（即 $t<\frac{n}{i*j}$）且 t 的最小质因数不小于 i 的 t 的个数。注意到 p[i][i * j] − v[i][i − 1] 即为所求。

然后分析复杂度，对所有需要枚举的 $\pi(\sqrt{n})$ 个素数 pk：

  • 有 $\pi(n^\frac{1}{4})$ 个素数 pk 满足 $i≤n^\frac{1}{4}$ ，需要转移 $\sqrt{n}$ 次，复杂度 $\pi(n^\frac{1}{4})\cdot\sqrt{n}=O(\frac{n^\frac{3}{4}}{\log{n}})$ ；
  • 有 $\pi(\sqrt{n})-\pi(n^\frac{1}{4})$ 个素数 pk 满足 $i>n^\frac{1}{4}$，需要转移 $\frac{n}{p_k^2}$ 次，复杂度 $\sum_{k=\pi(n^\frac{1}{4})}^{\pi(\sqrt{n})}{\frac{n}{p_k^2}}=O(\frac{n^\frac{3}{4}}{\log{n}})$

故算法总复杂度为 $O(\frac{n^\frac{3}{4}}{\log{n}})$ ，可以求 1011 量级的输入。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
i64 pi(i64 n) {  
    i64 m = ceil(sqrt(n));  
    vector<i64> p(m + 1), v(m + 1);  
    for (i64 i = 1; i <= m; i++) p[i] = n / i - 1, v[i] = i - 1;  
    for (i64 i = 2; i <= m; i++) {  
        if (v[i] == v[i - 1]) continue;  
        for (i64 j = 1; j <= min(m - 1, n / i / i); j++) {  
            if (i * j < m) p[j] -= p[i * j] - v[i - 1];  
            else p[j] -= v[n / i / j] - v[i - 1];  
        }  
        for (i64 j = m; j >= i * i; j--) v[j] -= v[j / i] - v[i - 1];  
    }  
    return p[1];  
}

树状数组优化 min_25 筛求素数计数函数

首先有复杂度为 $O(\frac{n^\frac{3}{4}}{\log{n}})$ 的经典 DP：

$$ G(n,k)=G(n,k−1)−G\left(\frac{n}{p_{k}}​,k−1 \right)+G(p_{k},k−1) $$

考虑根号分治，分治线为 B1 和  B2，令  B2 以下的 G 暴力统计，B1 以下的 G 使用树状数组维护，对于 B1 以上的 G 使用 DP。

可以证明，算法总复杂度为 $O((\frac{n}{\log{n}})^\frac{2}{3})$ ，可以在数秒内求 1013 量级的输入。

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
constexpr f80 eps = 1e-8;  
  
i64 pi(i64 n) {  
    if (n <= 2) return n == 2;  
    i64 tot = 1, l2 = sqrtl(n) + 1;  
    i64 lim = max<i64>({pow(n / __lg(n), 0.66), l2, 10000});  
    lim = min(lim, n);  
    while (lim * tot < n) tot++;  
    tot--;  
  
    vector<f80> inv(l2 + 1);  
    vector<i64> w(l2 + 1), g(tot + 1);  
    for (int i = 1; i <= l2; i++) w[i] = i - 1, inv[i] = 1.L / i;  
    for (int j = 1; j <= tot; j++) g[j] = n * inv[j] + eps - 1;  
  
    vector<char> vis(lim + 1);  
    Fenwick<int> ft(lim + 1);  
    vis[1] = true, ft.add(1, 1);  
    auto query = [&](i64 x) -> i64 {  
        x = min(x, lim);  
        return x - ft.sum(x);  
    };  
  
    for (i64 i = 2; i * i <= n; i++) {  
        if (vis[i]) continue;  
        i64 x0 = w[i - 1], t = n * inv[i] + eps, r = i * i;  
        i64 tl = min(tot, t / i), tl2 = min(tl, n / (l2 * i)), tl3 = min(tl2, tot / i);  
        for (i64 j = 1; j <= tl3; j++) g[j] -= g[j * i] - x0;  
        for (i64 j = tl3 + 1; j <= tl2; j++) g[j] -= query(t * inv[j] + eps) - x0;  
        for (i64 j = tl2 + 1; j <= tl; j++) g[j] -= w[t * inv[j] + eps] - x0;  
        for (i64 j = l2; j >= r; j--) w[j] -= w[j * inv[i] + eps] - x0;  
        if (r <= lim) {  
            for (i64 j = r; j <= lim; j += i) {  
                if (vis[j]) continue;  
                vis[j] = true, ft.add(j, 1);  
            }  
        }  
    }  
    if (tot == 0) return query(n);  
    return g[1];  
}

Meissel-Lehmer 算法

一个连模板题都过不去的 Meissel-Lehmer ，也许只能求 n = 1012 量级的素数计数函数 π(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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
constexpr int N = 3.2e6; // 视题目范围调整为 sqrt(n)
constexpr int M = 7; // 二次筛法层级  
constexpr int PM = 2 * 3 * 5 * 7 * 11 * 13 * 17; // 前 M 个质数乘积  
  
bool np[N];  
int prime[N], pi[N];  
int phi[PM + 1][M + 1], sz[M + 1];  
  
int get_prime() {  
    int cnt = 0;  
    np[0] = np[1] = true;  
    pi[0] = pi[1] = 0;  
    for (int i = 2; i < N; ++i) {  
        if (!np[i]) prime[++cnt] = i;  
        pi[i] = cnt;  
        for (int j = 1; j <= cnt && i * prime[j] < N; ++j) {  
            np[i * prime[j]] = true;  
            if (i % prime[j] == 0) break;  
        }  
    }  
    return cnt;  
}  
  
void init() {  
    get_prime();  
    sz[0] = 1;  
    for (int i = 0; i <= PM; ++i) phi[i][0] = i;  
    for (int i = 1; i <= M; ++i) {  
        sz[i] = prime[i] * sz[i - 1];  
        for (int j = 1; j <= PM; ++j)  
            phi[j][i] = phi[j][i - 1] - phi[j / prime[i]][i - 1];  
    }  
}  
  
int sqrt2(i64 x) {  
    i64 r = (i64) sqrt(x - 0.1);  
    while (r * r <= x) ++r;  
    return (int) r - 1;  
}  
  
int sqrt3(i64 x) {  
    i64 r = (i64) cbrt(x - 0.1);  
    while (r * r * r <= x) ++r;  
    return (int) r - 1;  
}  
  
i64 get_phi(i64 x, int s) {  
    if (s == 0) return x;  
    if (s <= M) return phi[x % sz[s]][s] + x / sz[s] * phi[sz[s]][s];  
    if (x <= prime[s] * prime[s]) return pi[x] - s + 1;  
    if (x <= prime[s] * prime[s] * prime[s] && x < N) {  
        int s2x = pi[sqrt2(x)];  
        i64 ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;  
        for (int i = s + 1; i <= s2x; ++i) ans += pi[x / prime[i]];  
        return ans;  
    }  
    return get_phi(x, s - 1) - get_phi(x / prime[s], s - 1);  
}  
  
i64 get_pi(i64 x) {  
    if (x < N) return pi[x];  
    i64 ans = get_phi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1;  
    for (int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; ++i)  
        ans -= get_pi(x / prime[i]) - i + 1;  
    return ans;  
}  
  
i64 lehmer_pi(i64 x) {  
    if (x < N) return pi[x];  
    int a = lehmer_pi(sqrt2(sqrt2(x)));  
    int b = lehmer_pi(sqrt2(x));  
    int c = lehmer_pi(sqrt3(x));  
    i64 sum = get_phi(x, a) + (i64) (b + a - 2) * (b - a + 1) / 2;  
    for (int i = a + 1; i <= b; ++i) {  
        i64 w = x / prime[i];  
        sum -= lehmer_pi(w);  
        if (i > c) continue;  
        int lim = lehmer_pi(sqrt2(w));  
        for (int j = i; j <= lim; ++j)  
            sum -= lehmer_pi(w / prime[j]) - (j - 1);  
    }  
    return sum;  
}

组合数学

常用公式

组合数计算公式：$\binom{n}{m} = \binom{n-1}{m-1} + \binom{n-1}{m} = \frac{n!}{m!(n-m)!}$

组合数递推公式：$\binom{n}{m} = \frac{n}{n-m} \binom{n-1}{m} = \frac{n-m+1}{m} \binom{n}{m-1} = \frac{n}{m} \binom{n-1}{m-1}$

曲棍球棒恒等式：$\sum_{i=m}^n\binom{i}{m} = \binom{n+1}{m+1}$

组合数计算

递推法可以直接计算单个组合数 $\binom{n}{k}$，时间复杂度 O(min (k, n − k)) 。

1
2
3
4
5
6
7
i64 C(int n, int k) {
    if (k < 0 || k > n) return 0;
    k = min(k, n - k);
    i64 res = 1;
    for (int i = 0; i < k; ++i) res = res * (n - i) / (i + 1);
    return res;
}

线性预处理出阶乘及其逆元，可以在 O(1) 的时间复杂度内计算组合数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
constexpr int maxa = 100000; // 预处理范围  
  
vector<i64> A(maxa + 1);  
vector<i64> IA(maxa + 1);  
vector<i64> inv(maxa + 1);  
  
void init() {  
    A[0] = IA[0] = inv[1] = 1;  
    for (int i = 1; i <= maxa; i++) A[i] = A[i - 1] * i % p;  
    for (int i = 2; i <= maxa; i++) inv[i] = (p - p / i) * inv[p % i] % p;  
    for (int i = 1; i <= maxa; i++) IA[i] = IA[i - 1] * inv[i] % p;  
}
i64 C(int n, int k) {  
    if (k < 0 || k > n) return 0;  
    return A[n] * IA[k] % p * IA[n - k] % p;  
}

Lucas 定理

设 p 为素数，n 和 k 为非负整数，其 p 进制展开为：

n = n0 + n1p + n2p2 + … + nmpm

k = k0 + k1p + k2p2 + … + kmpm

其中 0 ≤ ni, ki < p ，则组合数模 p 满足：

$$ \dbinom{n}{k} \equiv \prod_{i=0}^m \dbinom{n_i}{k_i} \pmod{p} $$

此过程可以使用递归或循环实现，至多进行 logpn 次。

1
2
3
4
5
6
7
8
9
10
i64 lucas(i64 n, i64 k, i64 p) {
    i64 res = 1;
    while (n > 0 || k > 0) {
        i64 ni = n % p, ki = k % p;
        if (ki > ni) return 0;
        res = res * C(ni, ki, p) % p;
        n /= p, k /= p;
    }
    return res;
}

容斥原理

容斥原理适用于解决重叠集合的计数问题，对于全集 U 下的 集合的并 可以使用容斥原理计算：

$$ \left|\bigcup_{i=1}^{n}S_i\right|=\sum_{m=1}^n(-1)^{m-1}\sum_{a_i<a_{i+1} }\left|\bigcap_{i=1}^mS_{a_i}\right| $$

集合的交则用全集减去补集的并集求得：

$$ \left|\bigcap_{i=1}^{n}S_i\right|=|U|-\left|\bigcup_{i=1}^n\overline{S_i}\right| $$

二项式反演

二项式反演通过交替容斥修正重复计数，将“至少”与“恰好”相互转化，解决“恰好满足若干条件”这类计数问题。

若 $f(n) = \sum_{k=0}^n \binom{n}{k} g(k)$，即至多满足钦定 n 个条件的计数，则反演为：

$$ g(n) = \sum_{k=0}^n (-1)^{n-k} \binom{n}{k} f(k) $$

若 $f(n) = \sum_{k=n}^m \binom{k}{n} g(k)$，即至少满足钦定 n 个条件的计数，则反演为：

$$ g(n) = \sum_{k=n}^m (-1)^{k-n} \binom{k}{n} f(k) $$

本质：通过交替容斥修正重复计数，将“至少”与“恰好”相互转化。例如，已知“至少 k 次满足条件（可能满足更多条件）”的函数 f(k)，通过可反演求出“恰好 k 次”的 g(k) 。

Min-Max 容斥

Min-Max 容斥是一种将最小值与最大值相互转化的技巧，在期望上也成立，常用于概率期望和计数问题中。其基本形式如下：

对于全集 U 中的元素 x1, x2, …, xn，有：

maxi ∈ Sxi = ∑T ⊆ S, T ≠ ⌀(−1)|T|−1minj ∈ Txj

mini ∈ Sxi = ∑T ⊆ S, T ≠ ⌀(−1)|T|−1maxj ∈ Txj

Kth Min-Max 容斥是 Min-Max 容斥的推广形式，用于求解第 k 大或第 k 小的值。其形式如下：

对于全集 U 中的元素 x1, x2, …, xn，第 k 大值（kthmax）和第 k 小值（kthmin）满足：

$$ \text{kthmax}_{i \in S} x_i = \sum_{T \subseteq S, |T| \geq k} (-1)^{|T| - k} \binom{|T| - 1}{k - 1} \min_{j \in T} x_j $$

$$ \text{kthmin}_{i \in S} x_i = \sum_{T \subseteq S, |T| \geq k} (-1)^{|T| - k} \binom{|T| - 1}{k - 1} \max_{j \in T} x_j $$

根据 Min-Max 容斥，我们还可以得到以下关于最小公倍数（LCM）和最大公约数（GCD）的等式：

lcmi ∈ Sxi = ∏T ⊆ S(gcdj ∈ Txj)(−1)|T|−1

其中，lcm  对应集合的 max ，gcd  对应集合的 min ，指数上的乘法 a1 对应加法 +，倒数 a−1 对应减法 − 。

第二类斯特林数

1
2
3
4
5
6
7
8
9
10
11
12
13
constexpr int maxn = 100000;  
constexpr int maxk = 100000;  
  
vector S(maxn + 1, vector<i64>(maxn + 1));  
  
void init() {  
    S[0][0] = 1;  
    for (int i = 1; i <= maxn; ++i) {  
        for (int j = 1; j <= min(i, maxk); ++j) {  
            S[i][j] = (S[i - 1][j - 1] + j * S[i - 1][j]) % p;  
        }  
    }  
}

多项式与生成函数

生成函数是将组合对象（如序列或结构）的计数信息编码为幂级数的系数，是组合数学中通过解析方法研究离散问题的核心工具。

在组合数学与算法分析中，多项式可以作为生成函数的载体，通过其代数运算（如乘法、系数提取）将离散序列的组合关系转化为解析形式，从而揭示序列规律、求解递推关系或枚举计数问题。

多项式

多项式常用于计算高精度乘法、生成函数计数、以及优化动态规划。

多项式基本运算

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
using poly = vector<i64>;

poly operator+(const poly &a, const poly &b) {  
    poly res(max(a.size(), b.size()));  
    for (int i = 0; i < a.size(); i++) res[i] = a[i];  
    for (int i = 0; i < b.size(); i++) res[i] = (res[i] + b[i]) % mod;  
    return res;  
}  
  
poly operator-(const poly &a, const poly &b) {  
    poly res(max(a.size(), b.size()));  
    for (int i = 0; i < a.size(); i++) res[i] = a[i];  
    for (int i = 0; i < b.size(); i++) res[i] = (res[i] - b[i] + mod) % mod;  
    return res;  
}  
  
poly operator*(const poly &a, const poly &b) {  
    poly res(a.size() + b.size() - 1);  
    for (int i = 0; i < a.size(); i++) {  
        for (int j = 0; j < b.size(); j++) {  
            res[i + j] = (res[i + j] + (i64) a[i] * b[j]) % mod;  
        }  
    }  
    return res;  
}  
  
poly shift(const poly &p, int k) {  
    poly res(k);  
    res.insert(res.end(), p.begin(), p.end());  
    return res;  
}

Karatsuba 算法

Karatsuba 算法利用分治策略将多项式乘法的时间复杂度优化至 O(nlog23) ，约为 O(n1.585)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
poly karatsuba(const poly &a, const poly &b) {  
    if (a.empty() || b.empty()) return poly();  
    int n = max(a.size(), b.size());  
    if (n <= 32) return a * b;  
  
    int k = n >> 1;  
    poly a0(a.begin(), a.begin() + min(k, (int) a.size()));  
    poly a1(a.begin() + min(k, (int) a.size()), a.end());  
    poly b0(b.begin(), b.begin() + min(k, (int) b.size()));  
    poly b1(b.begin() + min(k, (int) b.size()), b.end());  
  
    poly a0b0 = karatsuba(a0, b0);  
    poly a1b1 = karatsuba(a1, b1);  
    poly a0a1 = a0 + a1, b0b1 = b0 + b1;  
    poly mid = karatsuba(a0a1, b0b1) - a0b0 - a1b1;  
  
    poly res = shift(a1b1, 2 * k) + shift(mid, k) + a0b0;  
    return res;  
}

朴素循环卷积

以下代码实现了两个长度为 m 的向量在模 m 意义下的循环卷积（按模 m 下标相加的乘法），时间复杂度为 O(m2)。

稀疏多项式传入形参 a 可以有效优化执行速度，此时间复杂度为 O(nm)，其中 n 为实参 a 非零项的数量。

1
2
3
4
5
6
7
8
9
10
11
12
poly mul(const poly &a, const poly &b, i64 m) {  
    poly res(m);  
    for (i64 i = 0; i < m; ++i) {  
        if (a[i] == 0) continue;  
        for (i64 j = 0; j < m; ++j) {  
            if (b[j] == 0) continue;  
            i64 s = (i + j) % m;  
            res[s] = (res[s] + a[i] * b[j]) % p;  
        }  
    }  
    return res;  
}

拉格朗日插值

通过已知的 n 个点 (xi, yi)，构造一个 n − 1 次多项式 L(x)，使得 L(xi) = yi 。

其核心是基函数 $\ell_i(x)=\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}$ 的线性组合，即 L(x) = ∑yiℓi(x) 。

下列代码直接计算拉格朗日插值公式，对每个点 i 暴力求解基函数 ℓi(k) 的值，通过模逆元处理分母，复杂度 O(n2) 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
i64 lagrange(const poly &x, const poly &y, i64 k) {  
    int n = x.size();  
    i64 res = 0;  
    for (int i = 0; i < n; i++) {  
        i64 s1 = 1, s2 = 1;  
        for (int j = 0; j < n; j++) {  
            if (i == j) continue;  
            s1 = s1 * (k - x[j]) % p;  
            s2 = s2 * (x[i] - x[j]) % p;  
        }  
        res = (res + y[i] * s1 % p * qpow(s2, p - 2)) % p;  
    }  
    return (res + p) % p;  
}

拉格朗日插值法求多项式系数

以下代码可以计算拉格朗日插值多项式的系数，分为三部分，时间复杂度均为 O(n2) ：

  1. 分子多项式 num(X) = ∏i(X − xi) ：逐步乘入每个 (X − xi)，循环递推计算多项式系数；
  2. 分母 den[i] = yi ⋅ ∏j ≠ i(xi − xj)−1 mod  p ：直接计算每个点的拉格朗日基函数分母；
  3. 提取系数： 通过霍纳法则反向递推计算 $\frac{\text{num}(X)}{X - x_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
vector<i64> lagrange_coeff(const poly &x, const poly &y) {  
    int n = x.size();  
    vector<i64> num(n + 1), den(n, 1);  
    num[0] = 1;  
    for (int i = 0; i < n; ++i) {  
        for (int j = i + 1; j > 0; j--) {  
            num[j] = (num[j - 1] + num[j] * (p - x[i])) % p;  
        }  
        num[0] = num[0] * (p - x[i]) % p;  
        for (int j = 0; j < n; j++) {  
            if (i != j) den[i] = den[i] * (x[i] - x[j] + p) % p;  
        }  
        den[i] = y[i] * qpow(den[i], p - 2) % p;  
    }  
  
    poly res(n);  
    for (int i = 0; i < n; i++) {  
        i64 q = num.back();  
        for (int j = n - 1; j >= 0; j--) {  
            res[j] = (res[j] + den[i] * q) % p;  
            q = (num[j] + x[i] * q) % p;  
        }  
    }  
    return res;  
}

快速傅里叶变换

一个 n 次多项式可以被 n 个点唯一确定，当取固定 n 个自变量 x 时，我们可以很方便地实现乘法（只需要将对应位置值相乘即可）。

FFT 可以快速计算多项式 $A(x) = \sum_{k=0}^{n-1} a_k x^k$ 在 n 次单位根 ωnk = e2πik/n 处的值。

为方便，不妨令 n 为 2 的幂，有如下性质：

  • 对称性：ωnk + n/2 = −ωnk
  • 平方后仍为 2 的幂次单位根：(ωnk)2 = ωn/2k

任何一个多项式都可以表示为偶次项和奇次项两部分：

  • 偶次项： A0(x) = a0 + a2x + ⋯ + an − 2xn/2 − 1
  • 奇次项： A1(x) = a1 + a3x + ⋯ + an − 1xn/2 − 1
  • 原多项式：A(x) = A0(x2) + x ⋅ A1(x2)

通过单位根的对称性，可以将问题变为奇偶两个规模减半的子问题。

$$ \begin{align*} A(\omega_n^k) &= A_0(\omega_{n/2}^k) + \omega_n^k A_1(\omega_{n/2}^k) \\ A(\omega_n^{k+n/2}) &= A_0(\omega_{n/2}^k) - \omega_n^k A_1(\omega_{n/2}^k) \end{align*} $$

逆变换可以相同的逻辑实现，只需将单位根变为 $\omega_n^k = \frac{1}{n} e^{-2\pi i k / n}$ 。

朴素递归实现

使用标准库的复数和递归的最朴素实现，注意先将多项式不满 2 的幂次项补 0 ，时间复杂度 O(nlog n) 。

inv 传入 false 表示正变换，传入 true 表示逆变换，逆变换后需要将所有系数手动除以 n ；

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
using cpx = complex<f64>;  
using poly = vector<cpx>;  
  
void fft(vector<cpx> &a, int inv = false) {  
    int n = a.size();  
    if (n == 1) return;  
    vector<cpx> a1(n >> 1), a2(n >> 1);  
    for (int i = 0; i < n; i += 2) a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];  
    fft(a1, inv), fft(a2, inv);  
    cpx wn = {cos(2 * pi / n), sin(2 * pi / n)}, w = {1, 0};  
    if (inv) wn = conj(wn);  
    for (int i = 0; i < n >> 1; i++, w *= wn) {  
        a[i] = a1[i] + w * a2[i], a[i + (n >> 1)] = a1[i] - w * a2[i];  
    }  
}

迭代法实现

位逆序置换

以 8 项多项式为例，模拟拆分的过程：

  • 初始序列为 {x0, x1, x2, x3, x4, x5, x6, x7}
  • 一次二分之后 {x0, x2, x4, x6}, {x1, x3, x5, x7}
  • 两次二分之后 {x0, x4}, {x2, x6}, {x1, x5}, {x3, x7}
  • 三次二分之后 {x0}, {x4}, {x2}, {x6}, {x1}, {x5}, {x3}, {x7}

注意到最终那个位置的下标是原坐标的二进制翻转对称，称这个变换为位逆序置换。

考虑从小到大求 R(x)，首先 R(0) = 0 。在求 R(x) 时，$R\left(\left\lfloor \frac{x}{2} \right\rfloor\right)$ 的值是已知的。

将 x 右移一位，然后翻转，再右移一位，就得到了 x 除了个位之外其它位的翻转结果。

考虑个位的翻转结果：如果个位是 0，翻转之后最高位就是 0。如果个位是 1，则翻转后最高位是 1，因此还要加上 $\frac{len}{2} = 2^{k-1}$。综上

$$ R(x) = \left\lfloor \frac{R\left(\left\lfloor \frac{x}{2} \right\rfloor\right)}{2} \right\rfloor + (x \bmod 2) \times \frac{len}{2} $$

蝶形运算优化

已知 G(ωn/2k) 和 H(ωn/2k) 后，需要使用下面两个式子求出 f(ωnk) 和 f(ωnk + n/2) ：

f(ωnk) = G(ωn/2k) + ωnk × H(ωn/2k)

f(ωnk + n/2) = G(ωn/2k) − ωnk × H(ωn/2k)

使用位逆序置换后，对于给定的 n, k ：

  • G(ωn/2k) 的值存储在数组下标为 k 的位置，H(ωn/2k) 的值存储在数组下标为 $k + \frac{n}{2}$ 的位置。
  • f(ωnk) 的值将存储在数组下标为 k 的位置，f(ωnk + n/2) 的值将存储在数组下标为 $k + \frac{n}{2}$ 的位置。

因此可以直接在数组下标为 k 和 $k + \frac{n}{2}$ 的位置进行覆写，而不用开额外的数组保存值。此方法即称为 蝶形运算 。

使用以上方法优化后，可以实现一个无需递归的 FFT 实现。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
void fft(vector<cpx> &a, bool inv = false) {  
    int n = a.size(), bit = 0;  
    while (1 << bit < n) bit++;  
    vector<int> rev(n);  
    for (int i = 1; i < n; i++) {  
        rev[i] = rev[i >> 1] >> 1 | (i & 1) << (bit - 1);  
        if (i < rev[i]) swap(a[i], a[rev[i]]);  
    }  
    for (int mid = 1; mid < n; mid <<= 1) {  
        cpx wn = {cos(pi / mid), sin(pi / mid)};  
        if (inv) wn = conj(wn);  
        for (int i = 0; i < n; i += mid << 1) {  
            cpx w = {1, 0};  
            for (int j = 0; j < mid; j++, w = w * wn) {  
                cpx u = a[i + j], v = a[i + mid + j] * w;  
                a[i + j] = u + v, a[i + mid + j] = u - v;  
            }  
        }  
    }  
}

预处理优化

预处理位逆序置换表 rev 和单位根表 wf（正变换）、wb（逆变换），通过当前长度 2 的幂次 len 调用。

直接查表可以优化 fft 函数常数，在多次调用卷积时有常数优势。

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
constexpr int maxl = 4e6; // 最大长度的两倍
  
vector<vector<int> > rev;  
vector<vector<cpx> > wf, wb;  
  
void init() {  
    int len = 1;  
    while (len < maxl) len <<= 1;  
    rev.resize(len + 1);  
    for (int n = 2, bit = 0; n < len; n <<= 1, bit++) {  
        rev[n].resize(n);  
        for (int i = 1; i < n; i++) rev[n][i] = rev[n][i >> 1] >> 1 | (i & 1) << bit;  
    }  
    wf.resize((len >> 1) + 1), wb.resize((len >> 1) + 1);  
    for (int mid = 1; mid < len; mid <<= 1) {  
        wf[mid].resize(mid), wb[mid].resize(mid);  
        cpx w = {cos(pi / mid), sin(pi / mid)}, iw = conj(w);  
        wf[mid][0] = wb[mid][0] = 1;  
        for (int i = 1; i < mid; i++) {  
            wf[mid][i] = wf[mid][i - 1] * w;  
            wb[mid][i] = wb[mid][i - 1] * iw;  
        }  
    }  
}  
  
void fft(vector<cpx> &a, bool inv = false) {  
    int n = a.size();  
    for (int i = 1; i < n; i++) {  
        if (i < rev[n][i]) swap(a[i], a[rev[n][i]]);  
    }  
    for (int mid = 1; mid < n; mid <<= 1) {  
        auto &w = inv ? wb[mid] : wf[mid];  
        for (int i = 0; i < n; i += mid << 1) {  
            for (int j = 0; j < mid; j++) {  
                cpx u = a[i + j], v = a[i + j + mid] * w[j];  
                a[i + j] = u + v, a[i + j + mid] = u - v;  
            }  
        }  
    }  
}

FFT 优化多项式乘法

多项式乘法即卷积。多项式的系数表示称时域，值表示称频域。时域上的卷积等同频域上的点乘。

以下实现使用 FFT 将系数表示转化为值表示，执行 O(n) 的点乘，再将结果转化为系数表示，总时间复杂度 O(nlog n) 。

1
2
3
4
5
6
7
8
9
10
poly conv(poly a, poly b) {  
    int n = 1;  
    while (n < a.size() + b.size() - 1) n <<= 1;  
    a.resize(n), b.resize(n);  
    fft(a), fft(b);  
    for (int i = 0; i < n; i++) a[i] = a[i] * b[i];  
    fft(a, true);  
    for (int i = 0; i < n; i++) a[i] /= n;  
    return a;  
}

快速数论变换

NTT（快速数论变换）是 FFT 在有限域（模数意义下）的推广，使用数论单位根代替复数单位根，在算法竞赛中用于替代 FFT ，其原理可以参考 FFT 部分模板介绍。

迭代法实现

以下代码使用迭代法实现 NTT ，依赖快速幂 qpow，复杂度 O(nlog ⁡n) 。

inv 传入 false 表示正变换，传入 true 表示逆变换，逆变换后需要将所有系数手动乘以 n 的逆元；

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
constexpr i64 p = 998244353, G = 3;  
  
void ntt(poly &a, bool inv = false) {  
    int n = a.size(), bit = 0;  
    while (1 << bit < n) bit++;  
    vector<int> rev(n);  
    for (int i = 1; i < n; i++) {  
        rev[i] = rev[i >> 1] >> 1 | (i & 1) << (bit - 1);  
        if (i < rev[i]) swap(a[i], a[rev[i]]);  
    }  
    for (int mid = 1; mid < n; mid <<= 1) {  
        i64 wn = qpow(3, (p - 1) / (2 * mid), p);  
        if (inv) wn = qpow(wn, p - 2, p);  
        for (int i = 0; i < n; i += mid << 1) {  
            i64 w = 1;  
            for (int j = 0; j < mid; j++, w = w * wn % p) {  
                i64 u = a[i + j], v = a[i + j + mid] * w % p;  
                a[i + j] = (u + v) % p, a[i + j + mid] = (u - v + p) % p;  
            }  
        }  
    }  
}

预处理优化

预处理位逆序置换表 rev 和单位根表 wf（正变换）、wb（逆变换），通过当前长度 2 的幂次 len 调用。

直接查表可以优化 ntt 函数常数，在多次调用卷积时有常数优势。

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
constexpr int maxl = 4e6; // 最大长度的两倍  
constexpr i64 p = 998244353, G = 3;  
  
vector<vector<int> > rev;  
vector<vector<i64> > wf, wb;  
  
void init() {  
    int len = 1;  
    while (len < maxl) len <<= 1;  
    rev.resize(len + 1);  
    for (int n = 2, bit = 0; n < len; n <<= 1, bit++) {  
        rev[n].resize(n);  
        for (int i = 1; i < n; i++) rev[n][i] = rev[n][i >> 1] >> 1 | (i & 1) << bit;  
    }  
    wf.resize((len >> 1) + 1), wb.resize((len >> 1) + 1);  
    for (int mid = 1; mid < len; mid <<= 1) {  
        wf[mid].resize(mid), wb[mid].resize(mid);  
        i64 w = qpow(G, (p - 1) / (mid << 1), p), iw = qpow(w, p - 2, p);  
        wf[mid][0] = 1, wb[mid][0] = 1;  
        for (int i = 1; i < mid; ++i) {  
            wf[mid][i] = wf[mid][i - 1] * w % p;  
            wb[mid][i] = wb[mid][i - 1] * iw % p;  
        }  
    }  
}  
  
void ntt(poly &a, bool inv = false) {  
    int n = a.size();  
    for (int i = 1; i < n; ++i) {  
        if (i < rev[n][i]) swap(a[i], a[rev[n][i]]);  
    }  
    for (int mid = 1; mid < n; mid <<= 1) {  
        auto &w = inv ? wb[mid] : wf[mid];  
        for (int i = 0; i < n; i += mid << 1) {  
            for (int j = 0; j < mid; ++j) {  
                i64 u = a[i + j], v = a[i + j + mid] * w[j] % p;  
                a[i + j] = (u + v) % p, a[i + j + mid] = (u - v + p) % p;  
            }  
        }  
    }  
}

NTT 优化多项式乘法

以下实现使用 NTT 将系数表示转化为值表示，执行 O(n) 的点乘，再将结果转化为系数表示，总时间复杂度 O(nlog n) 。

1
2
3
4
5
6
7
8
9
10
11
poly conv(poly a, poly b) {  
    int n = 1;  
    while (n < a.size() + b.size() - 1) n <<= 1;  
    a.resize(n), b.resize(n);  
    ntt(a), ntt(b);  
    for (int i = 0; i < n; i++) a[i] = a[i] * b[i] % p;  
    ntt(a, true);  
    i64 invn = qpow(n, p - 2, p);  
    for (int i = 0; i < n; i++) a[i] = a[i] * invn % p;  
    return a;  
}

分治 NTT

以下代码使用 CDQ 分治策略结合多项式乘法逐层计算递推式 $f_i = \sum_{j=1}^i f_{i-j} g_j$ ，时间复杂度为 O(nklog n) ，其中 k 取决于选择的多项式乘法算法。

1
2
3
4
5
6
7
8
9
10
void solve(int l, int r, poly &f, poly &g) {  
    if (r - l <= 1) return;  
    int m = l + r >> 1;  
    solve(l, m, f, g);  
    poly A(f.begin() + l, f.begin() + m);  
    poly B(g.begin() + 1, g.begin() + min(r - l + 1, (int) g.size()));  
    poly C = karatsuba(A, B);  
    for (int i = m; i < r; i++) f[i] = (f[i] + C[i - l - 1]) % mod;  
    solve(m, r, f, g);  
}

概率论

含绝对值的期望

计算两个整数区间上随机选取点的差的绝对值的期望。

关键公式：$\frac{1}{n^2}\sum_{i=1}^n\sum_{j=1}^n|i-j| = \frac{n^2-1}{3n}$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
i64 eabs(i64 a, i64 b, i64 c, i64 d) {  
    if (a > c) swap(a, c), swap(b, d);  
    if (b > d) {  
        i64 el = eabs(a, d, c, d);  
        i64 pl = qdiv(d - a + 1, b - a + 1);  
        i64 er = eabs(d + 1, b, c, d);  
        i64 pr = qdiv(b - d, b - a + 1);  
        return (el * pl + er * pr) % p;  
    }  
    i64 n = max(b - c + 1, 1ll);  
    i64 em = qdiv(n * n - 1, 3 * n);  
    i64 pm = qdiv(n * n, (b - a + 1) * (d - c + 1));  
    return (qdiv(c + d - a - b, 2) + em * pm) % p;  
}

博弈论

博弈论是经济学的一个分支，主要研究具有竞争或对抗性质的个体，在特定规则下所产生的各种行为。博弈论关注博弈中个体的预期行为与实际行为，并研究其最优策略。

公平组合游戏

SG 定理

Nim 游戏

巴什博弈

威佐夫博弈

斐波那契 -Nim 游戏

有一堆个数为 n 的石子，先手第一次可以取任意个石子，但不能取完，之后每次取的石子个数最多只能是上一次的两倍，不能不取，取走最后一个石子的人胜。

结论：当 n 是斐波那契数时，先手必败，否则先手必胜。

证明：待更新

线性代数

矩阵快速幂

以下代码实现了矩阵乘法和矩阵快速幂，时间复杂度分别为 O(n³) 和 O(n³log 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
using poly = vector<i64>;  
using mat = vector<vector<i64> >;

mat mati(int n) {  
    mat res(n, poly(n));  
    for (int i = 0; i < n; i++) res[i][i] = 1;  
    return res;  
}  
  
mat operator*(const mat &a, const mat &b) {  
    assert(!a.empty() && a[0].size() == b.size());  
    int m = a.size(), s = b.size(), n = b[0].size();  
    mat res(m, poly(n));  
    for (int i = 0; i < m; i++)  
        for (int k = 0; k < s; k++)  
            for (int j = 0; j < n; j++)  
                res[i][j] = (res[i][j] + a[i][k] * b[k][j]) % p;  
    return res;  
}  
  
mat mpow(mat a, i64 n) {  
    mat res = mati(a.size());  
    while (n) {  
        if (n & 1) res = res * a;  
        a = a * a;  
        n >>= 1;  
    }  
    return res;  
}

高斯消元

将方程写为系数矩阵形式，以下操作称为初等行变换，初等行变换不改变解：

  • 两方程互换，解不变；
  • 一方程乘以非零数 k，解不变；
  • 一方程加上另一方程，解不变。

高斯消元逐列消元，每次选定一行非零元素为主元，利用初等行变换消去其它行元素，以下实现将矩阵化为行最简型，时间复杂度 O(nmmin (n, m))。

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
constexpr f64 eps = 1e-8;  
  
int sgn(auto x) {  
    if (abs(x) < eps) return 0;  
    return x > 0 ? 1 : -1;  
}  
  
mat gauss(mat a) {  
    int n = a.size(), m = a[0].size(), rk = 0;  
    for (int c = 0; c < m && rk < n; c++) {  
        int r = rk;  
        for (int i = rk + 1; i < n; i++) {  
            if (abs(a[i][c]) > abs(a[r][c])) r = i;  
        }  
        if (!sgn(a[r][c])) continue;  
        swap(a[rk], a[r]);  
        for (int j = m - 1; j >= c; j--) a[rk][j] /= a[rk][c];  
        for (int i = 0; i < n; i++) {  
            if (i != rk && sgn(a[i][c])) {  
                f64 fac = a[i][c];  
                for (int j = c; j < m; j++) {  
                    a[i][j] -= fac * a[rk][j];  
                }  
            }  
        }  
        rk++;  
    }  
    return a;  
}

解线性方程组

解线性方程组 Ax = b ，增广矩阵为方程组系数矩阵 A 与常数列 b 的并生成的新矩阵，将增广矩阵化为行最简型，即可直接读出特解和基础解系，方程组的所有解均为特解加上基础解系的线性组合，存在基础解析说明方程组具有无穷解。

以下实现中，返回的矩阵第一行表示特解，其余行均为基础解系，若方程无解则返回空矩阵，复杂度瓶颈在消元法。

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
mat calc(const mat &aug) {  
    int n = aug.size(), m = aug[0].size() - 1;  
    mat rref = gauss(aug);  
  
    vector<int> piv(m, -1);  
    for (int i = 0; i < n; i++) {  
        int j = 0;  
        while (j < m && !sgn(rref[i][j])) j++;  
        if (j == m && sgn(rref[i][m])) return {};  
        if (j != m) piv[j] = i;  
    }  
  
    mat res(1, poly(m));  
    for (int j = 0; j < m; j++) {  
        if (~piv[j]) {  
            res[0][j] = rref[piv[j]][m];  
            continue;  
        }  
        res.emplace_back(m, 0.0);  
        poly &cur = res.back();  
        cur[j] = 1.0;  
        for (int k = 0; k < m; k++) {  
            if (~piv[k]) cur[k] = -rref[piv[k]][j];  
        }  
    }  
    return res;  
}

求逆矩阵

对于方阵 A，若存在方阵 A−1，使得 A × A−1 = A−1 × A = I，则称矩阵 A 可逆，A−1 被称为它的逆矩阵。

给出 n 阶方阵 A，求解其逆矩阵的方法如下：

  1. 构造 n × 2n 的增广矩阵 (A, In)；
  2. 用高斯消元法将其化简为最简形 (In, A−1)，即可得到 A 的逆矩阵 A−1。如果最终最简形的左半部分不是单位矩阵 In，则矩阵 A 不可逆。

以下代码实现了模意义下的高斯消元和矩阵求逆。

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
mat gauss(mat a) {  
    int n = a.size(), m = a[0].size(), rk = 0;  
    for (int c = 0; c < m && rk < n; c++) {  
        int r = rk;  
        for (int i = rk + 1; i < n; i++) {  
            if (abs(a[i][c]) > abs(a[r][c])) r = i;  
        }  
        if (!a[r][c]) continue;  
        swap(a[rk], a[r]);  
        i64 inv = qpow(a[rk][c], P - 2, P);  
        for (int j = c; j < m; j++) a[rk][j] = a[rk][j] * inv % P;  
        for (int i = 0; i < n; i++) {  
            if (i != rk && a[i][c]) {  
                i64 fac = a[i][c];  
                for (int j = c; j < m; j++) {  
                    a[i][j] = ((a[i][j] - fac * a[rk][j]) % P + P) % P;  
                }  
            }  
        }  
        rk++;  
    }  
    return a;  
}  
  
mat inv(const mat &a) {  
    int n = a.size();  
    if (n == 0 || a[0].size() != n) return {};  
    mat aug(n, poly(2 * n));  
    for (int i = 0; i < n; i++) {  
        for (int j = 0; j < n; j++) aug[i][j] = a[i][j];  
        aug[i][n + i] = 1;  
    }  
    aug = gauss(aug);  
    for (int i = 0; i < n; i++) {  
        for (int j = 0; j < n; j++) {  
            if (aug[i][j] != (i == j)) return {};  
        }  
    }  
    mat res(n, poly(n));  
    for (int i = 0; i < n; i++) {  
        for (int j = 0; j < n; j++) {  
            res[i][j] = aug[i][j + n];  
        }  
    }  
    return res;  
}

异或线性基

称线性空间 V  的一个极大线性无关组为 V 的一组 Hamel 基 或 线性基，简称 基 。

将 n 维实线性空间 Rn 下的线性基称为实数线性基，n 维布尔域线性空间 Z2n 下的线性基称为异或线性基。

构造方法

使用高斯消元化线性基矩阵为行最简形式，构造的异或线性基具有如下性质：

  • 线性基没有异或和为 0 的子集。
  • 线性基中各数二进制最高位不同。
  • 如果某一位是一个元素的最高有效位，则其他元素在这一位均为 0 。
1
2
3
4
5
6
7
vector<u64> base;
    
auto insert = [&](u64 x) {
    for (auto b: base) x = min(x, b ^ x);
    for (auto &b: base) b = min(b, b ^ x);
    if (x) base.push_back(x);
};

因为 bitset 不提供 < 操作符，可以使用用 base[n] 存最高有效位为 n 的元素。

1
2
3
4
5
6
7
8
9
10
11
12
bitset<N> base[N];

auto insert = [&](bitset<N> x) {
    for (int i = 0; i < N; i++)
        if (x[i]) x ^= base[i];
    if (x == 0) return;
    int msb = N - 1; // 最高有效位
    while (msb && x[msb] == 0) msb--;
    base[msb] = x;
    for (int i = msb + 1; i < N; ++i)
        if (base[i][msb]) base[i] ^= x;
};

前缀线性基

前缀线性基允许对于序列的每个前缀，都维护该前缀的所有后缀的线性基，这样就可以支持查询每个区间的线性基，甚至可以利用树上差分扩展到树上。

注意到序列的某个前缀 [1, i] 的所有后缀 [j, i] 的线性基是相互包含的，这些后缀的线性基中互不相同的至多只有 n 种。利用这个单调性，只需要为添加的每个向量 v，都标记它出现的最大下标 t，就可以在 O(n) 的空间内存储所有后缀的线性基。查询区间 [i, j] 对应的线性基时，只需要在 i 处的前缀线性基中仅保留标记 t ≥ j 的那些向量即可。

要维护线性基中每个向量 v 的时间戳，只需要贪心地选取尽可能新的向量替换掉旧的向量即可：

  • 从高位向低位扫，如果 v 的第 x 位是 1 ，就比较线性基中已有的向量 ax 的时间戳 tx 和当前时间 i ：
    • 如果 i > tx，即要添加的向量时间更晚，就将 ax 设为 v，并更新时间戳为 i，并将 ax ⊕ v 按照之前记录的时间 tx 继续添加过程；
    • 如果 i < tx，即要添加的向量时间更早，不更新 ax 和 tx，将 v 异或 ax 后继续添加。

注意在更新位置 x 的向量时，不能将异或的结果 ax ⊕ v 存入位置 x，因为异或的结果 ax ⊕ v 的时间戳为 min (t(ax), t(v)) = t(ax)，小于要添加的变量 v 的时间戳 t(v)。同样的原因，高斯消元法构造线性基的过程中向上更新时可能会破坏时间戳的性质，所以不再适用于构造前缀线性基。如果需要用到高斯消元法得到的线性基的性质，可以在查询时另行处理。

利用前缀线性基可以离线解决线性基问题，也可以将每个前缀处的前缀线性基都存下来再查询，这可以看作是一种「可持久化」线性基。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
struct PreBasis {  
    static constexpr int mxb = 20;  
    array<pair<u64, int>, mxb> base;  
  
    void insert(u64 x, int np) {  
        for (int i = mxb; i >= 0 && x; i--) {  
            if (~x >> i & 1) continue;  
            auto &[b, p] = base[i];  
            if (p < np) swap(b, x), swap(p, np);  
            x = min(x, b ^ x);  
        }  
    }  
};

vector<PreBasis> pb(n + 1);  
for (int i = 1; i <= n; i++) {  
    u64 v;  
    cin >> v;  
    pb[i] = pb[i - 1];  
    pb[i].insert(v, i);  
}

第 k 小异或和

对于某个元素，选后异或和总是比不选异或和更大，所以最大异或和就是整个集合的异或和。

进一步地，排序后还可求出第 k 小异或和，只需按照 k 的二进制分解去选择基数即可获得第 k 小异或和。

1
2
3
4
5
6
7
8
9
10
sort(base.begin(), base.end());

u64 ans = 0;
if (base.size() < n) k--; // 存在基被完全消元，即线性基非空子集可异或出 0
for (auto b : base) {
    if (k & 1) ans ^= b;
    k >>= 1;
}
if (k == 0) cout << ans << endl;
else cout << -1 << endl;

线性规划

线性规划研究线性约束条件下线性目标函数最值问题，下文中标准形式定义如下：

有 n 个实数非负变量和 m 条约束，其中第 i 条约束形如 $\sum_{j=1}^n​ a_{i,j}​x_{j}​\leq b_{i}​$ 。

在满足所有约束的情况下，求每个变量 xj​ 的取值，使得目标函数 $F=\sum_{j=1}^n​ c_{j}​x_{j}$​ 的值最大。

等价转换

通过以下方法，可以将任意线性规划问题转化为标准形式，规模不超过原问题的规模的二倍：

原问题中的情况 转换方法 转换后的标准型
目标函数为最小化：
$\min F = \sum_{j=1}^n c_j x_j$		 令 F′ = −F
转化为最大化问题		 $\max F' = \sum_{j=1}^n (-c_j) x_j$
约束为 ≥ 形式：
$\sum_{j=1}^n a_{i,j} x_j \geq b_i$		 两边乘以 −1 $\sum_{j=1}^n (-a_{i,j}) x_j \leq -b_i$
约束为 = 形式：
$\sum_{j=1}^n a_{i,j} x_j = b_i$		 拆分为两个约束：
$\sum_{j=1}^n a_{i,j} x_j \leq b_i$
$\sum_{j=1}^n a_{i,j} x_j \geq b_i$		 两个约束：
$\sum_{j=1}^n a_{i,j} x_j \leq b_i$
$\sum_{j=1}^n (-a_{i,j}) x_j \leq -b_i$
变量无符号限制：
xj ∈ ℝ		 引入两个非负变量：
xj = xj+ − xj−
其中 xj+ ≥ 0, xj− ≥ 0		 替换原变量：
xj → xj+ − xj−
新增变量 xj+, xj− 满足非负约束

单纯形法

单纯形法通过在可行域的顶点间迭代移动，每次选择使目标函数值增加最快的入基变量和保持可行性的出基变量进行主元旋转，直到找到最优解或判断问题无界。

以下代码实现了单纯形法求解标准形式线性规划问题， solve 函数传入目标函数系数向量 C、约束条件系数矩阵 A 和右侧向量 B，返回求解状态（−1 为不可行，0 为无界，1 为可行且有解）、最优目标函数值和最优解 。

单次转轴时间复杂度为 O(mn)，通常认为在 2m 次转轴内可以得到大多数问题的最优解，但最坏情况下需要指数次数的转轴。

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
constexpr f80 eps = 1e-8;  
  
int sgn(auto x) {  
    if (abs(x) < eps) return 0;  
    return x > 0 ? 1 : -1;  
}  
  
struct Simplex {  
    int cntr, cntc, neg;  
    vector<int> idxP, idxQ;  
    vector<vector<f80> > tbl;  
  
    f80 eval(int x, int y) {  
        if (sgn(tbl[x][y]) >= 0) return 0;  
        f80 s = 0;  
        for (int i = 0; i < cntc; ++i) s += tbl[x][i] * tbl[x][i];  
        return tbl[x][y] * tbl[x][y] / s;  
    }  
  
    void pivot(int x, int y) {  
        swap(idxP[x], idxQ[y]);  
        f80 v = -1.L / tbl[x][y];  
        for (int j = 0; j < cntc + 2; ++j)  
            tbl[x][j] = (j == y) ? -v : v * tbl[x][j];  
        for (int i = 0; i <= cntr; ++i) {  
            if (i == x) continue;  
            v = tbl[i][y];  
            tbl[i][y] = 0;  
            for (int j = 0; j < cntc + 2; ++j)  
                tbl[i][j] += v * tbl[x][j];  
        }  
    }  
  
    bool init() {  
        while (neg) {  
            int x = -1;  
            f80 ms = 0;  
            for (int i = 0; i < cntr; ++i) {  
                f80 s = eval(i, cntc + 1);  
                if (s > ms) ms = s, x = i;  
            }  
            if (x < 0) return false;  
            int y = -1, sy = -1;  
            for (int j = 0; j < cntc; ++j) {  
                int sj = sgn(tbl[x][j]);  
                if (sj == (idxQ[j] < 0 ? -1 : 1)) {  
                    if (y < 0) y = j, sy = sj;  
                    else {  
                        int s = sj * sy * sgn(tbl[cntr][j] * tbl[x][y] - tbl[cntr][y] * tbl[x][j]);  
                        if (!s) s = (idxQ[j] < idxQ[y]) ? -1 : 1;  
                        if (s < 0) y = j, sy = sj;  
                    }  
                }  
            }  
            if (idxQ[y] < 0) {  
                --neg;  
                idxQ[y] = ~idxQ[y];  
                pivot(x, y);  
                tbl[x][cntc + 1] += 1;  
            } else {  
                pivot(x, y);  
            }  
        }  
        return true;  
    }  
  
    bool simplex() {  
        while (true) {  
            int x = -1;  
            f80 ms = 0;  
            for (int i = 0; i < cntr; ++i) {  
                f80 s = eval(i, cntc);  
                if (s > ms) ms = s, x = i;  
            }  
            if (x < 0) return true;  
            int y = -1;  
            for (int j = 0; j < cntc; ++j) {  
                if (sgn(tbl[x][j]) > 0) {  
                    if (y < 0) y = j;  
                    else {  
                        int s = sgn(tbl[cntr][j] * tbl[x][y] - tbl[cntr][y] * tbl[x][j]);  
                        if (!s) s = (idxQ[j] < idxQ[y]) ? -1 : 1;  
                        if (s < 0) y = j;  
                    }  
                }  
            }  
            if (y < 0) return false;  
            pivot(x, y);  
        }  
    }  
  
    tuple<int, f80, vector<f80> > solve(const vector<f80> &C, const vector<vector<f80> > &A, const vector<f80> &B) {  
        cntr = C.size();  
        cntc = A.size();  
        tbl.assign(cntr + 1, vector<f80>(cntc + 2));  
        idxP.resize(cntr);  
        idxQ.resize(cntc);  
        neg = 0;  
  
        for (int i = 0; i < cntr; ++i) {  
            for (int j = 0; j < cntc; ++j)  
                tbl[i][j] = A[j][i];  
            tbl[i][cntc] = -C[i];  
        }  
        for (int j = 0; j < cntc; ++j)  
            tbl[cntr][j] = B[j];  
  
        iota(idxP.begin(), idxP.end(), 0);  
        iota(idxQ.begin(), idxQ.end(), cntr);  
  
        for (int j = 0; j < cntc; ++j) {  
            if (sgn(tbl[cntr][j]) < 0) {  
                for (int i = 0; i < cntr; ++i)  
                    tbl[i][cntc + 1] += tbl[i][j];  
                idxQ[j] = ~idxQ[j];  
                ++neg;  
            }  
        }  
  
        if (!init()) return make_tuple(-1, f80(), vector<f80>());  
        if (!simplex()) return make_tuple(0, f80(), vector<f80>());  
  
        vector<f80> z(cntr);  
        for (int i = 0; i < cntc; ++i)  
            if (idxQ[i] < cntr)  
                z[idxQ[i]] = tbl[cntr][i];  
  
        return make_tuple(1, tbl[cntr][cntc], z);  
    }  
};

数值算法

自适应辛普森算法

自适应辛普森算法的思想是将被积区间分为若干小段，每段套用二次函数的积分公式进行计算，通过选取合适的 eps 可以控制运行时间和精度的平衡。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
f80 a;  
  
f80 f(f80 x) {  
    return pow(x, a / x - x);  
}  
  
f80 sim(f80 l, f80 r) {  
    return (f(l) + 4 * f((l + r) / 2) + f(r)) * (r - l) / 6;  
}  
  
f80 asr(f80 l, f80 r, f80 e, f80 st) {  
    f80 m = (l + r) / 2, sl = sim(l, m), sr = sim(m, r);  
    if (abs(sl + sr - st) <= 15 * e) return sl + sr + (sl + sr - st) / 15;  
    return asr(l, m, e / 2, sl) + asr(m, r, e / 2, sr);  
}  
  
f80 asr(f80 l, f80 r) {  
    return asr(l, r, eps, sim(l, r));  
}

计算几何

二维几何

1
2
3
4
5
6
7
8
9
const f80 eps = 1e-8; // 根据题目精度要求进行修改  
const f80 pi = acos(-1.0L); // 3.1415926....  
  
// 取符号函数  
int sgn(auto x) {  
    // 精度范围内的近似相等  
    if (abs(x) < eps) return 0;  
    return x > 0 ? 1 : -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
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
template<typename T>  
struct pos {  
    T x, y;  
  
    pos operator+(const pos &b) const { return {x + b.x, y + b.y}; }  
  
    pos operator-(const pos &b) const { return {x - b.x, y - b.y}; }  
  
    T operator^(const pos &b) const { return x * b.y - y * b.x; }  
  
    T operator*(const pos &b) const { return x * b.x + y * b.y; }  
  
    pos operator*(const T &b) const { return {x * b, y * b}; }  
  
    pos operator/(const T &b) const { return {x / b, y / b}; }  
  
    bool operator==(const pos &b) const { return !sgn(x - b.x) && !sgn(y - b.y); }  
  
    [[nodiscard]] T abs2() const { return x * x + y * y; }  
};  
  
using Pi = pos<i64>;  
using Pf = pos<f80>;  
  
// 获得点 a 到点 b 的距离  
f80 dist(auto a, auto b) {  
    return sqrtl((a - b).abs2());  
}  
  
// 获得点 a 的单位化向量  
auto norm(auto a) {  
    return a / hypotl(a.x, a.y);  
}  
  
// 判断点 p 在向量 ab 的哪一侧  
int side(auto p, auto a, auto b) {  
    return sgn((a - p) ^ (b - p));  
}  
  
// 判断点 p 是否在线段 ab 上  
bool ons(auto p, auto a, auto b) {  
    return side(p, a, b) == 0 && sgn((a - p) * (b - p)) <= 0;  
}  
  
// 判断线段 ab 与线段 cd 是否存在交点  
bool iss(auto a, auto b, auto c, auto d) {  
    if (ons(a, c, d) || ons(b, c, d) || ons(c, a, b) || ons(d, a, b)) return true;  
    if (side(a, b, c) * side(a, b, d) >= 0) return false;  
    if (side(c, d, a) * side(c, d, b) >= 0) return false;  
    return true;  
}  
  
// 判断直线 ab 与线段 cd 是否存在交点  
bool ils(auto a, auto b, auto c, auto d) {  
    return side(a, b, c) * side(a, b, d) <= 0;  
}  
  
// 判断直线 ab 与直线 cd 是否平行  
bool para(auto a, auto b, auto c, auto d) {  
    return sgn(a - b ^ c - d) == 0;  
}  
  
// 获得直线 ab 与直线 cd 的交点  
Pf ill(auto a, auto b, auto c, auto d) {  
    auto u = b - a, v = d - c;  
    auto t = (a - c ^ v) / (v ^ u);  
    return a + u * t;  
}  
  
// 获取点 p 在直线 uv 上的垂足  
Pf foot(auto p, auto u, auto v) {  
    auto dir = norm(u - v);  
    return u + dir * (dir * (p - u));  
}  
  
// 获取向量 a 逆时针旋转角度 r 后的向量， r ∈ [0, PI]
Pf rot(auto a, f80 r) {  
    Pf b = {sinl(r), cosl(r)};  
    return {a ^ b, a * b};  
}  
  
// 获取三角形 ABC 的面积  
f80 tria(auto a, auto b, auto c) {  
    return abs(b - a ^ c - a) / 2;  
}

  
// 获取闭合多边形 p 的面积  
f80 polya(auto &p) {  
    f80 res = 0;  
    for (int i = 0; i < p.size(); ++i) res += p[i] ^ p[(i + 1) % p.size()];  
    return abs(res) / 2;  
}

直线的运算

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
struct line {  
    Pf s, e;  
  
    Pf dir() const { return e - s; }  
  
    f80 ang() const { return atan2((e - s).y, (e - s).x); }  
};  
  
// 判断点 p 在有向直线 a 的哪一侧  
int side(auto p, line a) {  
    return sgn((a.s - p) ^ (a.e - p));  
}  
  
// 判断直线 ab 与直线 cd 是否平行  
bool para(line a, line b) {  
    return sgn(a.dir() ^ b.dir()) == 0;  
}  
  
// 获得直线 a 与直线 b 的交点  
pos<f80> ill(line a, line b) {  
    return ill(a.s, a.e, b.s, b.e);  
}  
  
// 判断直线 l 与线段 cd 是否存在交点  
bool ils(auto l, auto c, auto d) {  
    return side(l.s, l.e, c) * side(l.s, l.e, d) <= 0;  
}  
  
// 获取点 p 在直线 l 上的垂足  
pos<f80> foot(auto p, line l) {  
    auto dir = norm(l.s - l.e);  
    return l.s + dir * (dir * (p - l.s));  
}

凸包算法

Andrew 算法

计算给定点集的凸包，返回按逆时针顺序排列的凸包顶点，时间复杂度为  O(nlog n) 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
template<typename T>
vector<T> hull(vector<T> p) {
    sort(p.begin(), p.end(), [&](auto &a, auto &b) {
        return a.x < b.x || a.x == b.x && a.y < b.y;
    });
    vector<T> h, l;
    for (auto a: p) {
        while (h.size() > 1 && side(a, h.back(), h[h.size() - 2]) <= 0) h.pop_back();
        while (l.size() > 1 && side(a, l.back(), l[l.size() - 2]) >= 0) l.pop_back();
        l.push_back(a), h.push_back(a);
    }
    reverse(h.begin(), h.end());
    l.insert(l.end(), h.begin() + 1, h.end() - 1);
    return l;
}
动态凸壳算法

基于 multiset 动态维护上凸壳周长，插入时均摊时间复杂度为 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
struct Hull : public multiset<Pi> {  
    f80 C = 0;  
  
    void modify(iterator t, int op) {  
        auto l = t == begin() ? end() : prev(t), r = next(t);  
        if (l != end() && r != end()) C += (dist(*l, *t) + dist(*t, *r) - dist(*l, *r)) * op;  
        else if (l != end()) C += dist(*l, *t) * op;  
        else if (r != end()) C += dist(*t, *r) * op;  
    }  
  
    void append(Pi p) {  
        auto t = insert(p);  
        if (t != begin() && next(t) != end() && side(*t, *prev(t), *next(t)) <= 0) {  
            erase(t);  
            return;  
        }  
        modify(t, 1);  
        while (t != begin() && prev(t) != begin() && side(*t, *prev(t), *prev(t, 2)) <= 0) {  
            modify(prev(t), -1), erase(prev(t));  
        }  
        while (next(t) != end() && next(t, 2) != end() && side(*t, *next(t), *next(t, 2)) >= 0) {  
            modify(next(t), -1), erase(next(t));  
        }  
    }  
};
半平面交算法

使用双端队列维护可行半平面，计算一组半平面的交集，返回按逆时针顺序排列的凸包或凸壳，时间复杂度为  O(nlog 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
vector<Pf> half(vector<line> ls) {  
    sort(ls.begin(), ls.end(), [](auto &x, auto &y) {  
        if (x.ang() != y.ang()) return x.ang() < y.ang();  
        return (x.s ^ x.e) < (y.s ^ y.e);  
    });  
  
    auto check = [](line a, line b, line c) {  
        return side(ill(b, c), a.s, a.e) > 0;  
    };  
  
    deque<line> dq;  
    for (auto &li: ls) {  
        while (dq.size() >= 2 && !check(li, dq[dq.size() - 2], dq.back())) dq.pop_back();  
        while (dq.size() >= 2 && !check(li, dq[0], dq[1])) dq.pop_front();  
        if (dq.empty() || sgn((dq.back().dir() ^ li.dir())) != 0) dq.push_back(li);  
    }  
    while (dq.size() >= 2 && !check(dq[0], dq[dq.size() - 2], dq.back())) dq.pop_back();  
    while (dq.size() >= 2 && !check(dq.back(), dq[0], dq[1])) dq.pop_front();  
  
    vector<Pf> poly;  
    if (dq.size() < 2) return poly;  
    for (int i = 0; i < dq.size(); ++i) poly.push_back(ill(dq[i], dq[(i + 1) % dq.size()]));  
    return poly;  
};

高精度计算

高精度无符号整型

p 进制压 18 位高精度无符号整型，部分中间运算使用 i128 保证不溢出。

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
struct wintu : private vector<i64> {  
    static constexpr i64 p = 1000000000000000000;  
    static constexpr int w = 18; // 压 18 位高精度，相当于以十进制位数为变量的复杂度除 18  
    wintu() : vector(1) {  
    }  
  
    explicit(false) wintu(const string &s) {  
        for (int i = (int) s.size(); i >= 1; i -= w) {  
            int begin = max(0, i - w);  
            push_back(stoll(s.substr(begin, i - begin)));  
        }  
    }  
  
    explicit(false) wintu(u64 n) {  
        if (n == 0) push_back(0);  
        while (n > 0) push_back((i64) (n % p)), n /= p;  
    }  
  
    friend ostream &operator<<(ostream &out, const wintu &self) {  
        out << self[self.size() - 1];  
        for (int i = (int) self.size() - 2; i >= 0; --i)  
            out << setw(w) << setfill('0') << self[i];  
        return out;  
    }  
  
    explicit operator bool() const {  
        return size() != 1 || operator[](0);  
    }  
  
    explicit operator string() const {  
        stringstream sio;  
        sio << *this;  
        return sio.str();  
    }  
  
    explicit operator u64() const {  
        return stoull((string) *this);  
    }  
  
    bool operator<(const wintu &other) const {  
        if (size() != other.size()) return size() < other.size();  
        for (int i = (int) size() - 1; i >= 0; --i) {  
            if (operator[](i) != other[i]) return operator[](i) < other[i];  
        }  
        return false;  
    }  
  
    bool operator>(const wintu &other) const { return other < *this; }  
    bool operator<=(const wintu &other) const { return !(*this > other); }  
    bool operator>=(const wintu &other) const { return !(*this < other); }  
    bool operator==(const wintu &other) const { return *this <= other && other <= *this; }  
  
    wintu operator+(const wintu &other) const {  
        wintu res("");  
        long long t = 0;  
        for (int i = 0; i < max(size(), other.size()); ++i) {  
            if (i < size()) t = t + operator[](i);  
            if (i < other.size()) t += other[i];  
            res.push_back(t % p);  
            t /= p;  
        }  
        if (t) res.push_back(t);  
        return res;  
    }  
  
    wintu operator-(const wintu &other) const {  
        wintu res("");  
        long long t = 0;  
        for (int i = 0; i < size(); ++i) {  
            t = operator[](i) - t;  
            if (i < other.size()) t -= other[i];  
            res.push_back((t + p) % p);  
            if (t < 0) t = 1;  
            else t = 0;  
        }  
        while (res.size() > 1 && res.back() == 0) res.pop_back();  
        return res;  
    }  
  
    wintu operator*(const wintu &other) const {  
        wintu res("");  
        res.resize(size() + other.size());  
        for (int i = 0; i < size(); ++i) {  
            for (int j = 0; j < other.size(); ++j) {  
                i128 tmp = res[i + j] + (i128) operator[](i) * other[j];  
                res[i + j + 1] += (i64) (tmp / p);  
                res[i + j] = (i64) (tmp % p);  
            }  
        }  
        while (res.size() > 1 && res.back() == 0) res.pop_back();  
        return res;  
    }  
  
    wintu operator/(int b) const {  
        wintu res("");  
        i64 r = 0;  
        for (int i = (int) size() - 1; i >= 0; --i) {  
            r = r * p + operator[](i);  
            res.push_back(r / b);  
            r %= b;  
        }  
        reverse(res.begin(), res.end());  
        while (res.size() > 1 && res.back() == 0) res.pop_back();  
        return res;  
    }  
  
    i64 operator%(i64 b) const {  
        i64 r = 0;  
        for (int i = (int) size() - 1; i >= 0; --i) r = (r * p + operator[](i)) % b;  
        return r;  
    }  
};

高精度有符号整型

基于高精度无符号整型 wintu 实现的高精度有符号整型，运算规则与 c++ 内建整型行为一致。

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
struct wint {  
    bool neg;  
    wintu val;  
  
    wint() : neg(false), val(0) {  
    }  
  
    explicit(false) wint(string s) {  
        if (s.front() == '-') {  
            neg = true;  
            s.erase(s.begin());  
        } else {  
            neg = false;  
        }  
        val = s;  
    }  
  
    explicit(false) wint(i64 n) {  
        if (n < 0) {  
            neg = true;  
            n = -n;  
        } else {  
            neg = false;  
        }  
        val = n;  
    }  
  
    friend ostream &operator<<(ostream &out, const wint &self) {  
        if (self.neg) out << '-';  
        out << self.val;  
        return out;  
    }  
  
    explicit operator bool() const {  
        return (bool) val;  
    }  
  
    explicit operator string() const {  
        stringstream sio;  
        sio << *this;  
        return sio.str();  
    }  
  
    explicit operator i64() const {  
        i64 res = (u64) val;  
        return neg ? -res : res;  
    }  
  
    bool operator<(const wint &other) const {  
        if (neg != other.neg) return neg;  
        if (neg) return other.val < val;  
        return val < other.val;  
    }  
  
    bool operator>(const wint &other) const { return other < *this; }  
    bool operator<=(const wint &other) const { return !(*this > other); }  
    bool operator>=(const wint &other) const { return !(*this < other); }  
    bool operator==(const wint &other) const { return neg == other.neg && val == other.val; }  
  
    wint operator+(const wint &other) const {  
        wint res;  
        if (neg == other.neg) {  
            res.neg = neg;  
            res.val = val + other.val;  
        } else {  
            if (val >= other.val) {  
                res.neg = neg;  
                res.val = val - other.val;  
            } else {  
                res.neg = other.neg;  
                res.val = other.val - val;  
            }  
        }  
        return res;  
    }  
  
    wint operator-(const wint &other) const {  
        wint res;  
        if (neg != other.neg) {  
            res.neg = neg;  
            res.val = val + other.val;  
        } else {  
            if (val >= other.val) {  
                res.neg = neg;  
                res.val = val - other.val;  
            } else {  
                res.neg = !neg;  
                res.val = other.val - val;  
            }  
        }  
        return res;  
    }  
  
    wint operator*(const wint &other) const {  
        wint res;  
        res.neg = neg != other.neg;  
        res.val = val * other.val;  
        return res;  
    }  
  
    wint operator/(int b) const {  
        wint res;  
        res.neg = neg != b < 0;  
        res.val = val / abs(b);  
        return res;  
    }  
  
    int operator%(int b) const {  
        int r = val % abs(b);  
        return neg ? -r : r;  
    }  
};

取模类

封装了加减乘除等基本运算的整数取模类，支持输入输出、欧几里得法求逆元和快速幂。

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
template<const int MOD>  
struct ModInt {  
    int x;  
  
    ModInt() : x(0) {  
    }  
  
    explicit(false) ModInt(i64 y) : x(y >= 0 ? y % MOD : (y % MOD + MOD) % MOD) {  
    }  
  
    ModInt &operator+=(const ModInt &p) {  
        if ((x += p.x) >= MOD) x -= MOD;  
        return *this;  
    }  
  
    ModInt &operator-=(const ModInt &p) {  
        if ((x += MOD - p.x) >= MOD) x -= MOD;  
        return *this;  
    }  
  
    ModInt &operator*=(const ModInt &p) {  
        x = (int) (1LL * x * p.x % MOD);  
        return *this;  
    }  
  
    ModInt &operator/=(const ModInt &p) {  
        *this *= p.inv();  
        return *this;  
    }  
  
    ModInt operator-() const { return ModInt(-x); }  
    ModInt operator+(const ModInt &p) const { return ModInt(*this) += p; }  
    ModInt operator-(const ModInt &p) const { return ModInt(*this) -= p; }  
    ModInt operator*(const ModInt &p) const { return ModInt(*this) *= p; }  
    ModInt operator/(const ModInt &p) const { return ModInt(*this) /= p; }  
  
    bool operator==(const ModInt &p) const { return x == p.x; }  
    bool operator!=(const ModInt &p) const { return x != p.x; }  
  
    ModInt inv() const {  
        int a = x, b = MOD, u = 1, v = 0;  
        while (b) {  
            int t = a / b;  
            swap(a -= t * b, b);  
            swap(u -= t * v, v);  
        }  
        return ModInt(u);  
    }  
  
    ModInt pow(i64 n) const {  
        ModInt res = 1, a = x;  
        while (n) {  
            if (n & 1) res *= a;  
            a *= a;  
            n >>= 1;  
        }  
        return res;  
    }  
  
    friend ostream &operator<<(ostream &os, const ModInt &p) {  
        return os << p.x;  
    }  
  
    friend istream &operator>>(istream &is, ModInt &a) {  
        i64 t;  
        is >> t;  
        a = ModInt(t);  
        return is;  
    }  
};

using mint = ModInt<mod>;