数学与计算

数论

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

素数筛

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

埃拉托斯特尼筛

如果我们从小到大考虑每个数，然后同时把当前这个数的所有（比自己大的）倍数记为合数，那么运行结束的时候没有被标记的数就是素数了，注意 0 和 1 要特殊处理，时间复杂度 $O(n\log{\log{n}})$。

注意到筛法求素数的同时也得到了每个数的最小质因子，因此筛法可以在 $O(n\log{n})$ 的时间复杂度内对 $[1,n]$ 区间所有数分解质因数，或者在 $O(k\log{n})$ 的时间复杂度内将 $[1,n]$ 区间的所有数分解出至多 $k$ 个质因数。

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
vector<int> p, f1, f2(r - l + 1);  
auto eratosthenes = [&](int maxn) {  
    f1.resize(maxn + 1);  
    for (i64 i = 2; i <= maxn; i++) {  
        if (f1[i]) continue;  
        p.push_back(i);  
        for (i64 j = i * i; j <= maxn; j += i) {  
            if (!f1[j]) f1[j] = i;  
        }  
        for (i64 j = max(i, (l + i - 1) / i) * i; j <= r; j += i) {  
            if (!f2[j - l]) f2[j - l] = i;  
        }  
    }  
};  
eratosthenes(ceil(sqrt(r)) + 1);

欧拉筛

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

欧拉筛不保证每个一定被合数被最小值质因数标记，因此不可以筛至平方根。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
vector<int> pr, vis;  
  
auto eular = [&](int maxn) {  
    vis.resize(maxn + 1);  
    for (i64 i = 2; i <= maxn; i++) {  
        if (!vis[i]) pr.push_back(i);  
        for (auto x: pr) {  
            if (i * x > maxn) break;  
            vis[i * x] = true;  
            if (i % x == 0) break;  
        }  
    }  
};  
eular(n);

分层 DP 求素数计数函数

每个合数的最小质因数不大于 $\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]$ ，则划去的数字个数可以表示为满足 $ti≤\frac{n}{j}$（即 $t<\frac{n}{ij}$）且 $t$ 的最小质因数不小于 $i$ 的 $t$ 的个数。注意到 $p[i][i*j]-v[i][i-1]$ 即为所求。

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

  • 有 $\pi(n^\frac{1}{4})$ 个素数 $p_k$ 满足 $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})$ 个素数 $p_k$ 满足 $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}})$ ，可以求 $10^{11}$ 量级的输入。

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

数论分块

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

1
2
3
4
5
6
7
8
9
i64 n;  
cin >> n;  
i64 l = 1, ans = 0;  
while (l <= n) {  
    i64 r = n / (n / l);  
    ans += n / l * (r - l + 1);  
    l = r + 1;  
}  
cout << ans << endl;

素数与质因数

Miller-Rabin 素性测试

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

可在底数列表中添加随机数以进行更多轮试验，常用于对高精度数进行测试，此时时间复杂度是 $O(k\log^{3}n)$ ，其中 $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(n^{1/4}logn)$ 。

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

扩展欧几里得算法

扩展欧几里得算法可以求得 $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\equiv{c}\pmod{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. 求得原方程的一组特解 $x_{0}=c’x’$, $y_{0}=c’y’$

原方程的所有解可以表示为 $x=x_{0}+tb$, $y=y_{0}-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 筛是一种用于快速计算积性函数前缀和的亚线性筛法，其在 $10^{13}$ 数据范围内的复杂度为常数很小的 $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;  
};

Meissel-Lehmer 算法

一个连模板题都过不去的 Meissel-Lehmer ，也许只能求 $n=10^{12}$ 量级的素数计数函数 $\pi(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
17
constexpr int maxa = 100000; // 预处理范围  
  
vector<i64> A(maxa + 1);  
vector<i64> IA(maxa + 1);  
vector<i64> inv(maxa + 1);  
  
void init() {  
    A[0] = 1, inv[1] = 1, IA[0] = 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 = n_0 + n_1 p + n_2 p^2 + \dots + n_m p^m
$$

$$
k = k_0 + k_1 p + k_2 p^2 + \dots + k_m p^m
$$

其中 $0 \leq n_i, k_i < p$ ，则组合数模 $p$ 满足：

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

此过程可以使用递归或循环实现，至多进行 $\log_p{n}$ 次。

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

第二类斯特林数

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

多项式与生成函数

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

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

拉格朗日插值

通过已知的 $n$ 个点 $(x_i,y_i)$，构造一个 $n-1$ 次多项式 $L(x)$，使得 $L(x_i)=y_i$ 。

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
i64 lagrange(const vector<i64> &x, const vector<i64> &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(n^2)$ ：

  1. 分子多项式 $\text{num}(X) = \prod_{i}(X - x_i)$ ：逐步乘入每个 $(X - x_i)$，循环递推计算多项式系数；
  2. 分母 $\text{den}[i] = y_i \cdot \prod_{j \neq i}(x_i - x_j)^{-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 vector<i64> &x, const vector<i64> &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;  
    }  
  
    vector<i64> 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;  
}

快速数论变换

快速数论变换（NTT）是快速傅里叶变换（FFT）在有限域（模数意义下）的推广，使用数论单位根代替复数单位根，适用于模数满足特定条件的整数多项式乘法计算。

普通实现

以下代码实现快速数论变换，用于模数 p 下的多项式乘法，依赖快速幂 qpow，复杂度 $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
void ntt(poly &a, bool inv) {  
    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 root = qpow(3, (p - 1) / (2 * mid), p);  
        if (inv) root = qpow(root, p - 2, p);  
        for (int i = 0; i < n; i += mid << 1) {  
            i64 w = 1;  
            for (int j = 0; j < mid; ++j) {  
                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;  
                w = w * root % p;  
            }  
        }  
    }  
    if (!inv) return;  
    i64 invn = qpow(n, p - 2, p);  
    for (auto &x: a) x = x * invn % p;  
}

预处理优化

预处理逆序置换表 rev 和单位根表 wf（正变换）、wb（逆变换），优化 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
42
43
44
constexpr i64 maxl = 600; // 最大长度的两倍
constexpr i64 p = 998244353;  
constexpr i64 G = 3;  
  
vector<vector<i64> > rev, 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) {  
    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;  
            }  
        }  
    }  
    if (!inv) return;  
    i64 invn = qpow(n, p - 2, p);  
    for (auto &x: a) x = x * invn % p;  
}

形式幂级数

形式幂级数是以无穷级数形式表示的代数对象，其系数序列对应生成函数的离散信息，生成函数本质上是将组合序列编码为形式幂级数的系数。

形式幂级数乘法

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

稀疏多项式传入形参 $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;  
}

NTT 优化

以下代码实现基于 快速数论变换（NTT） 的多项式乘法，包含两个函数：

  1. conv(a, b)：计算两个多项式 aa 和 bb 的标准循环卷积（结果长度为 $len(a)+len(b)−1$）。
  2. pmul(a, b, m)：计算 $a$ 和 $b$ 的模 $m$ 循环卷积（即多项式乘法后系数对 $m$ 取模求和）。

依赖 ntt（快速数论变换）进行加速，时间复杂度均为 $O(n\log{⁡n})$ ，其中 $n$ 是变换长度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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, false), ntt(b, false);  
    for (int i = 0; i < n; ++i) a[i] = a[i] * b[i] % p;  
    ntt(a, true);  
    return a;             
}  
  
poly pmul(const poly &a, const poly &b, int m) {  
    poly c = conv(a, b), res(m);  
    for (int i = 0; i < c.size(); ++i) {  
        res[i % m] = (res[i % m] + c[i]) % p;  
    }  
    return res;  
}

概率论

含绝对值的期望

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

关键公式：$\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;  
}

线性代数

异或线性基

构造方法

本方法构造的异或线性基具有如下性质：

  • 线性基没有异或和为 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;
};

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

计算几何

二维几何

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

凸包算法

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

半平面交算法

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 进制压位高精度整型，默认压九位，部分中间运算应使用 i64 保证不溢出。

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 : private vector<int> {
    static constexpr int p = 1000000000, w = 9;

    wint() : vector(1) {
    }

    explicit(false) wint(const string &s) {
        for (int i = (int) s.size(); i >= 1; i -= w) {
            int begin = max(0, i - w);
            push_back(stoi(s.substr(begin, i - begin)));
        }
    }

    explicit(false) wint(u64 n) {
        if (n == 0) push_back(0);
        while (n > 0) push_back((int) (n % p)), n /= p;
    }

    friend ostream &operator<<(ostream &out, const wint &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 stoi((string) *this);
    }

    bool operator<(const wint &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 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 *this <= other && other <= *this; }

    wint operator+(const wint &other) const {
        wint res("");
        int 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;
    }

    wint operator-(const wint &other) const {
        wint res("");
        int 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;
    }

    wint operator*(const wint &other) const {
        wint res("");
        res.resize(size() + other.size());
        for (int i = 0; i < size(); ++i) {
            for (int j = 0; j < other.size(); ++j) {
                i64 tmp = res[i + j] + (i64) operator[](i) * other[j];
                res[i + j + 1] += (int) (tmp / p);
                res[i + j] = (int) (tmp % p);
            }
        }
        while (res.size() > 1 && res.back() == 0) res.pop_back();
        return res;
    }

    wint operator/(int b) const {
        wint res("");
        i64 r = 0;
        for (int i = (int) size() - 1; i >= 0; --i) {
            r = r * p + operator[](i);
            res.push_back((int) (r / b));
            r %= b;
        }
        reverse(res.begin(), res.end());
        while (res.size() > 1 && res.back() == 0) res.pop_back();
        return res;
    }

    int operator%(int b) const {
        i64 r = 0;
        for (int i = (int) size() - 1; i >= 0; --i) r = (r * p + operator[](i)) % b;
        return (int) r;
    }
};