竞赛算法(10) 组合数学

抽屉原理#

往3个抽屉里放一支铅笔, 必然会有一个抽屉中铅笔的数量大于等于2.

更一般地, 向n个抽屉中分配m支铅笔, 所有抽屉中铅笔的最大值大于等于m / n. 也叫鸽巢原理.

应用#

对于任意 $n$ 个整数 $a_1, a_2, \dots, a_n$ , 必然存在一个连续子段，其和是 $n$ 的倍数. 这是因为我们对这个数组模n再做前缀和, 前缀和有 $n + 1$ 种取值(注意 $S_0 = 0$ ), 所以必然有相同的两项, 则这两项之间的子段和模n为0, 即为n的倍数.
奇偶性: 在二维平面上给你 5 个整点(坐标都是整数), 证明其中必有一对点, 它们的中点也是整点. 注意到模 2 意义下的取值只有 4 种(奇偶, 偶奇, 奇奇, 偶偶), 却有 5 个点, 故一定有两个点的(各个维度的)奇偶性都相同, 故这两个点的中点必然是整点.
存在性证明: 证明某种解一定存在, 从而可以写暴力/贪心解决该问题.
优化搜索空间: 看似范围很大, 实则可以根据抽屉原理减小搜索范围.

容斥原理#

|A \cup B| = |A| + |B| - |A\cap B|

更一般地,

\left| \bigcup_{i=1}^n A_i \right| = \sum_{i=1}^n |A_i| - \sum_{1 \le i < j \le n} |A_i \cap A_j| + \sum_{1 \le i < j < k \le n} |A_i \cap A_j \cap A_k| - \dots + (-1)^{n-1} |A_1 \cap \dots \cap A_n|

即

\left| \bigcup_{i=1}^n A_i \right| = \sum_{\emptyset \neq J \subseteq \{1, \dots, n\}} (-1)^{|J|-1} \left| \bigcap_{j \in J} A_j \right|

示例: CSES-Prime Multiples

1
// CSES-Prime Multiples
2
void solution() {
3
    ll n;
4
    int k;
5
    cin >> n >> k;
6
    vector<ll> arr(k);
7
    for (int i = 0; i < k; i++) {
8
        cin >> arr[i];
9
    }
10
    int limit = 1 << k;
11
    ll ret = 0;
12
    for (int i = 1; i < limit; i++) {
13
        i128 now = 1;
14
        bool ok = 1;
15
        for (int j = 0; j < k; j++) {
16
            if ((i >> j) & 1) {
17
                now *= arr[j];
18
            }
19
            if (now > n) {
20
                ok = 0;
21
                break;
22
            }
23
        }
24
        int cnt = __builtin_popcount(i);
25
        int sig = (cnt % 2 == 0) ? -1 : 1;
26
        if (ok) {
27
            ret += sig * (n / now);
28
        }
29
    }
30
    cout << ret << '\n';
31
}

斐波那契数列#

F_0 = 0, F_1 = 1, F_i = F_{i - 1} + F_{i - 2}

常规求法 $O(n)$ , 利用倍增递推式或者矩阵快速幂可以优化到 $O(n\log n)$ . 倍增递推式:

\begin{aligned} F_{2k} &= F_k(2F_{k+1} - F_k) \\ F_{2k+1} &= F_{k+1}^2 + F_k^2 \end{aligned}

1
pair<int, int> fib(int n) {
2
    if (n == 0) {
3
        return { 0, 1 };
4
    }
5
    auto p = fib(n >> 1);
6
    int c = p.first * (2 * p.second - p.first);
7
    int d = p.first * p.first + p.second * p.second;
8
    if (n % 2 == 1) {
9
        return { d, c + d };
10
    } else {
11
        return { c, d };
12
    }
13
}

常见性质#

卡西尼性质:

F_{n-1}F_{n+1} - F_n^2 = (-1)^n

gcd性质:

\gcd(F_n, F_m) = F_{\gcd(n, m)}

相邻项互质, 斐波那契数列的相邻项是辗转相除法最差的起始情况.
前缀和性质:

\sum_{i=1}^n F_i = F_{n+2} - 1

奇偶项求和:

F_1 + F_3 + \dots + F_{2n-1} = F_{2n}\\ F_2 + F_4 + \dots + F_{2n} = F_{2n+1} - 1

加法规则:

F_{n + k} = F_kF_{n + 1} + F_nF_{k - 1}

这可以引出: 对于任意正整数k, $F_{nk}$ 是 $F_n$ 的倍数(反过来也成立, 若 $F_{m}$ 是 $F_n$ 的倍数, 则 $m$ 一定是 $n$ 的倍数).

齐肯多夫定理与斐波那契编码#

任意正整数n都能表示成若干不连续的斐波那契数的和. 于是我们可以使用这个性质给正整数编码. 即n表示为 $d_0d_1d_2...d_s1$ . 注意这里末尾手动添加一个1, 这样会强制出现两个连续的1, 表示编码结束. $d_i$ 若为1表示 $F_{i + 2}$ 被使用, 否则表示 $F_{i + 2}$ 未使用.

给n编码的过程可以使用贪心算法:

从大到小找到最大的 $F_i <= n$ , 添加1, 从n中减去 $F_i$ .
继续往下找, $F_i > n$ 添加0, 否则添加1并减去 $F_i$ . 直到 $n$ 为0.
在末尾添加1表示结束.

解码过程同理.

Pisano周期#

斐波那契数列模任意正整数m下会形成周期数列. 证明如下: 我们知道斐波那契数列只要确定了相邻两项, 就能确定后续所有项. 模m意义下任意相邻两项只会有 $m^2$ 个取值, 所以对于前 $m^2 + 1$ 个相邻项来说, 必然会有两个相同的相邻项. 那么从这两项出发形成的斐波那契数列相同, 即形成周期, 最小正周期长度小于等于 $m^2$ .

Pisano证明了这个周期数列的最小正周期长度不会超过6m, 称为Pisano周期. OEIS A001175

二项式系数(组合数)的计算#

杨辉三角#

C(n, m) % p = (C(n - 1, m) + C(n - 1, m - 1)) % p, 没什么好说的, 复杂度 $O(n^2)$ .

1
mint _C[N][N];
2
constexpr int N = 2026;
3
void init_pascal() {
4
    static bool inited = 0;
5
    if (inited) {
6
        return;
7
    }
8
    for (int i = 0; i < N; i++) {
9
        _C[i][0] = _C[i][i] = 1;
10
        for (int j = 1; j < i; j++) {
11
            _C[i][j] = _C[i - 1][j] + _C[i - 1][j - 1];
12
        }
13
    }
14
    inited = 1;
15
}
16
mint nCr(ll n, ll m) {
17
    if (m < 0 || m > n) {
18
        return 0;
19
    }
20
    init_pascal();
21
    return _C[n][m];
22
}

线性预处理逆元#

求出前n个整数的阶乘和阶乘逆元fact[n + 1], factinv[n + 1].
使用定义C(n, m) = fact[n] * factinv[m] * factinv[n - m]即可.

模数应为大于n的素数, 否则无法保证逆元存在.

复杂度 $O(n), O(1)$

1
constexpr int N = 2e5 + 50;
2
ll MOD = 998244353;
3
vector<ll> fact, factinv;
4
void init_fact() {
5
    static bool inited = 0;
6
    if (inited) {
7
        return;
8
    }
9
    int n = N;
10
    inited = 1;
11
    ll p = MOD;
12
    fact.resize(n + 1);
13
    factinv.resize(n + 1);
14
    fact[0] = 1;
15
    for (int i = 1; i <= n; i++) {
16
        fact[i] = mul(fact[i - 1], i, p);
17
    }
18
    factinv[n] = inv(fact[n], p);
19
    for (int i = n - 1; i >= 0; i--) {
20
        factinv[i] = mul(factinv[i + 1], i + 1, p);
21
    }
22
    inited = 1;
23
}
24

25
ll nCr(ll n, ll m) {
26
    if (m < 0 || m > n) {
27
        return 0;
28
    }
29
    ll p = MOD;
30
    init_fact();
31
    ll ret = fact[n];
32
    ret = mul(ret, factinv[m], p);
33
    ret = mul(ret, factinv[n - m], p);
34
    return ret;
35
}

直接计算#

注意到C(n, m) = _fact[n] / (_fact[m] * _factinv[n - m]) = (n * (n - 1) * ... * (n - m + 1)) / (m * (m - 1) * ... * 1)(上面m项, 下面m项). 所以只需要 $O(m + log P)$ 的复杂度即可直接计算. 在n很大, m很小的时候很好用.

1
ll MOD = 998244353;
2
ll nCr(ll n, ll m, ll p = MOD) {
3
    if (m < 0 || m > n) {
4
        return 0;
5
    }
6
    if (m > n / 2) {
7
        m = n - m;
8
    }
9
    ll A = 1;
10
    ll B = 1;
11
    for (ll i = 0; i < m; i++) {
12
        A = mul(A, n - i, p);
13
        B = mul(B, m - i, p);
14
    }
15
    return mul(A, inv(B, p), p);
16
}

Lucas定理#

C(n, m) % p = C(n / p, m / p) * C(n % p, m % p). 在p较小的时候有用.

1
ll Lucas(ll n, ll m, ll p) {
2
    if (m == 0) {
3
        return 1;
4
    }
5
    return mul(Lucas(n / p, m / p, p), nCr(n % p, m % p, p), p); // 通常使用预处理版本的nCr, 注意这里的nCr要求在m > n时返回0, 别忘记了.
6
}

高精度求组合数#

直接用BigInt乘超级慢. 这里介绍一种使用Legendre定理优化的方法.

先筛n以下的质数.
注意到C(n, m) = fact[n] / (fact[m] * fact[n - m]. 且fact[n]中一定包含且只包含所有小于等于n的质数作为质因数. 我们枚举小于等于n的质数p, 处理出fact[n]中p的幂次, 减去fact[m]中的幂次, 再减去fact[n - m]中的幂次, 即为C(n, m)中质数p的幂次.
快速得到fact[i]的幂次可以使用Legendre定理, 即fact[i]中p的幂次为i / p + i / p^2 + ... + i / p^k. 这是个有限项的求和.

1
vector<int> primes;
2
int N = 2e6 + 50;
3
void sieve(int n = N) {
4
    // 埃氏筛或线性筛
5
}
6
ll legendre(ll n, ll p) {
7
    // 返回 n! 中, p的幂次
8
    ll ret = 0;
9
    while (n > 0) {
10
        ret += n / p;
11
        n /= p;
12
    }
13
    return ret;
14
}
15
BigInt nCr(ll n, ll m) {
16
    if (m < 0 || m > n) {
17
        return 0;
18
    }
19
    BigInt ret = 1;
20
    for (int i = 0; primes[i] <= n; i++) {
21
        int p = primes[i];
22
        int exp = legendre(n, p) - legendre(m, p) - legendre(n - m, p);
23
        ret *= BigInt(p).pow(exp);
24
    }
25
    return ret;
26
}

概率和期望#

本章涵盖概率dp和期望dp的内容.

期望的线性性#

计算期望最有用的工具, 比如计算某集合 $A$ 中满足某性质的子集 $B$ 的元素个数的期望是多少. 可以设定指示器变量 $I_i$ (若 $a[i]$ 在 $B$ 中为1, 否则为0). 求 $B$ 元素个数期望即转化为求 $E[\sum_i I_i]$ . 由于期望的线性型, 我们可以把 $\sum_i I_i$ 拆开, 转化成求 $\sum_iE[I_i]$ , 设 $p_i$ 为 $a[i]$ 在 $B$ 中的概率, $E[I_i] = 1 \cdot p_i + 0 \cdot (1 - p_i) = p_i$ . 这时候对 $p_i$ 求和即为所求结果.

例题: 牛客126319F

简单来记的话就是”概率dp正推, 期望dp倒推”. 这是因为算概率的时候通常需要计算从起点到某状态的可能性, 而算期望的时候通常计算从当前状态到终点的期望代价.

概率的转移(概率dp)#

我们把状态看作点(每个状态拥有一个状态概率 $p_x$ ), 状态的转移看作有向边, 终点为 $y$ , 每一条有向边赋予一个转移概率(流量) $p_{xy}$ . 所有初始状态的概率和为1. 若存在终止状态, 则所有终止状态的概率和也为1. 每个节点的所有转移概率的和为1. 大致可以分成以下三类.

DAG上#

假设当前在 $x$ 点, 状态概率为 $p_x$ . 对于每条出边概率 $p_{xy}$ 有 $p_y += p_{xy} \times p_x$ . 同理可以倒过来得到 $p_x$ 等于所有入边概率乘以前驱的状态概率.

带自环的DAG上#

假设当前 $x$ 点状态概率为 $p_x$ , 自环转移概率为 $p'$ , 对于除自环外的出边 $y$ , $p_{xy}$ 来说, 我们需要对所有概率进行归一化, 即所有出边的概率都要除以 $1 - p'$ , 这时候就消除了自环. (很好理解, 因为如果 $p'$ 小于1, 那就不可能一直走自环, 总有一天会走到另外的出边上). 所以 $p_y += \frac{p_xp_{xy}}{1 - p'}$ . 注意到这种情况是不能倒过来用前驱来计算当前节点的, 因为每条入边归一化的分母 $p'$ 都不一样(应为前驱节点的 $p'$ , 而非本节点的 $p'$ ).

可能带环的有向图#

需要比较复杂的算法, 待补充.

期望的转移(期望dp)#

DAG上#

先计算终点期望, 然后反推到起点.

对于每条出边y, p:

$E[x] = \sum_{(x, y) \in E} p_{xy} \cdot (E[y] + w_{xy})$

其中 $p_{xy}$ 是从 $x$ 转移到 $y$ 的概率, $w_{xy}$ 是这次转移的代价.

带自环的DAG上#

假设点 $x$ 有 $p'$ 的概率走自环(代价 $w'$ ).有 $p_{xi}$ 的概率走到 $y_i$ (代价 $w_{xi}$ ):

$E[x] = p'(E[x] + w') + \sum p_{xi}(E[y_i] + w_{xi})$

移项化简可得: $E[x] = \frac{p' w' + \sum p_{xi}(E[y_i] + w_{xi})}{1 - p'}$ 注: 若每步代价均为 1, 则公式简化为 $E[x] = \frac{1 + \sum p_{xi} E[y_i]}{1 - p'}$ .

可能带环的有向图#

也就是马尔科夫链. 需要高斯消元, 待补充.

trick合集#

任何一个排列, 都可以被唯一地分解为若干个”极大连续递增块”的拼接.