📖 第 8.5 章 ⏱️ 约 55 分钟 🎯 Gold / 困难

第 8.5 章：组合数学与数论

📝 前置要求： 本章需要基础代数知识以及附录 E（数学基础）中的模运算内容。第 6.1 章（DP）也很有帮助，因为许多组合问题是通过 DP 解决的。

组合数学与数论题目在 USACO Gold 中定期出现——通常作为每套题的"数学题"。它们一般涉及计数（有多少种合法方案？）或整除性（哪些数满足某个性质？），答案通常需要对一个素数取模（通常为 10⁹+7）。

学习目标：

  • 正确实现模运算（加、乘、快速幂、逆元）
  • 高效计算 C(n, k) mod p
  • 应用容斥原理
  • 使用埃氏筛分解质因数
  • 识别并解决 USACO 计数问题

8.5.0 为什么需要模运算？

组合数学的答案可能极其庞大。C(100, 50) 有 30 位数字！USACO 题目始终要求输出答案对某个素数 p（通常 p = 10⁹+7 = 1,000,000,007）取模的结果。

模运算的核心规则：

  • (a + b) mod p = ((a mod p) + (b mod p)) mod p ✓
  • (a × b) mod p = ((a mod p) × (b mod p)) mod p ✓
  • (a − b) mod p = ((a mod p) − (b mod p) + p) mod p ✓（加 p 防止负数）
  • (a / b) mod p = (a mod p) × (b⁻¹ mod p) mod p ← 需要模逆元
const long long MOD = 1e9 + 7;

long long mod_add(long long a, long long b) { return (a + b) % MOD; }
long long mod_sub(long long a, long long b) { return (a - b + MOD) % MOD; }
long long mod_mul(long long a, long long b) { return (a % MOD) * (b % MOD) % MOD; }

8.5.1 快速幂（二进制取幂）

用反复平方法在 O(log n) 内计算 a^n mod p：

📄 用反复平方法在 O(log n) 内计算 a^n mod p： 
// 返回 a^n mod p
long long power(long long a, long long n, long long p = MOD) {
    a %= p;
    long long result = 1;
    while (n > 0) {
        if (n & 1)             // 若 n 当前位为 1
            result = result * a % p;
        a = a * a % p;         // 底数平方
        n >>= 1;               // 移动到下一位
    }
    return result;
}

示例： 2^10 = 2^(1010₂) = 2^8 × 2^2 = 256 × 4 = 1024

8.5.2 模逆元

要计算 a/b mod p，需要 b 的模逆元：一个值 b⁻¹ 满足 b × b⁻¹ ≡ 1 (mod p)。

何时存在？ 仅当 gcd(b, p) = 1 时。若 p 为素数且 0 < b < p，逆元始终存在。

方法一：费马小定理（p 必须为素数）

费马小定理：对素数 p 且 gcd(a,p)=1，有 a^(p−1) ≡ 1 (mod p)。

因此：a^(p−2) ≡ a⁻¹ (mod p)。

long long mod_inv(long long a, long long p = MOD) {
    return power(a, p - 2, p);  // O(log p)
}

// 模意义下的除法：
long long mod_div(long long a, long long b) {
    return mod_mul(a, mod_inv(b));
}

方法二：扩展欧几里得算法（适用于非素数模数）

📄 C++ 完整代码 
// 返回 x 使得 a*x ≡ 1 (mod m)
// 用扩展 GCD：找 x, y 满足 a*x + m*y = gcd(a, m)
long long ext_gcd(long long a, long long b, long long& x, long long& y) {
    if (b == 0) { x = 1; y = 0; return a; }
    long long x1, y1;
    long long g = ext_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

long long mod_inv_general(long long a, long long m) {
    long long x, y;
    long long g = ext_gcd(a, m, x, y);
    if (g != 1) return -1;  // 逆元不存在
    return (x % m + m) % m;
}

8.5.3 C(n, k) mod p 的计算

二项式系数 C(n, k) = n! / (k! × (n−k)!) 计算从 n 个中选取 k 个的方案数。

预处理阶乘（适用于多次查询，n ≤ 10⁶）

📄 查看代码：预处理阶乘（适用于多次查询，n ≤ 10⁶） 
const int MAXN = 1000001;
const long long MOD = 1e9 + 7;

long long fact[MAXN], inv_fact[MAXN];

void precompute_factorials(int n) {
    fact[0] = 1;
    for (int i = 1; i <= n; i++)
        fact[i] = fact[i-1] * i % MOD;

    inv_fact[n] = power(fact[n], MOD - 2);  // 费马小定理
    for (int i = n - 1; i >= 0; i--)
        inv_fact[i] = inv_fact[i+1] * (i+1) % MOD;
    // inv_fact[i] = 1/i! mod p，倒序计算
}

long long C(int n, int k) {
    if (k < 0 || k > n) return 0;
    return fact[n] % MOD * inv_fact[k] % MOD * inv_fact[n-k] % MOD;
}

帕斯卡三角 DP（适用于 n, k ≤ 2000 的小规模情形）

long long dp[2001][2001];  // dp[n][k] = C(n, k) mod p

void precompute_pascal(int maxn) {
    for (int i = 0; i <= maxn; i++) {
        dp[i][0] = 1;
        for (int j = 1; j <= i; j++)
            dp[i][j] = (dp[i-1][j-1] + dp[i-1][j]) % MOD;
    }
}
// C(n, k) = dp[n][k]

8.5.4 常用组合公式

公式数值含义
C(n, k)n! / (k!(n-k)!)从 n 个中无序选 k 个（不重复）
P(n, k)n! / (n-k)!从 n 个中有序选 k 个（不重复）
n^kn^k将 k 个不同物品放入 n 个不同盒子
C(n+k-1, k)(n+k-1)! / (k!(n-1)!)隔板法：k 个物品放入 n 个盒子（可重复）
n! / (a! b! c! ...)多项式：排列 a 种类型的物品
C(2n, n) / (n+1)卡特兰数二叉树、合法括号序列

卡特兰数（在 USACO 中出现频率出人意料地高）：

long long catalan(int n) {
    // C_n = C(2n, n) / (n+1)
    return C(2*n, n) % MOD * mod_inv(n+1) % MOD;
}
// C_0=1, C_1=1, C_2=2, C_3=5, C_4=14, C_5=42, ...

8.5.5 容斥原理

容斥原理通过交替加减计算集合并集的大小：

|A₁ ∪ A₂ ∪ ... ∪ Aₙ| = Σ|Aᵢ| − Σ|Aᵢ ∩ Aⱼ| + Σ|Aᵢ ∩ Aⱼ ∩ Aₖ| − ...

2～3 个集合的模板：

|A ∪ B| = |A| + |B| − |A ∩ B|
|A ∪ B ∪ C| = |A| + |B| + |C| − |A∩B| − |A∩C| − |B∩C| + |A∩B∩C|

USACO 题型：统计"至少满足一个条件"的序列

"统计长度为 N 的序列（每个元素取 1..M），使得 1..K 中每个值至少出现一次的方案数。"

对"缺失值"做容斥：

总数 = Σ_{j=0}^{K} (-1)^j × C(K, j) × (M-j)^N
  • 选 j 个值排除（C(K, j) 种方案）
  • 用剩余 M-j 个值填满 N 个位置：(M-j)^N 种序列
  • 容斥的交替符号
long long count_surjective(int n, int m, int k) {
    // 统计每个 K 值都至少出现一次的 N 元序列数
    long long ans = 0;
    for (int j = 0; j <= k; j++) {
        long long term = C(k, j) * power(m - j, n) % MOD;
        if (j % 2 == 0) ans = (ans + term) % MOD;
        else             ans = (ans - term + MOD) % MOD;
    }
    return ans;
}

8.5.6 埃氏筛

O(N log log N) 时间内找出 N 以内所有质数：

📄 O(N log log N) 时间内找出 N 以内所有质数： 
const int MAXN = 1000001;
bool is_prime[MAXN];
vector<int> primes;

void sieve(int n) {
    fill(is_prime, is_prime + n + 1, true);
    is_prime[0] = is_prime[1] = false;

    for (int i = 2; i <= n; i++) {
        if (is_prime[i]) {
            primes.push_back(i);
            for (long long j = (long long)i * i; j <= n; j += i)
                is_prime[j] = false;
        }
    }
}

线性筛 O(N)——用于质因数分解

📄 查看代码：线性筛 O(N)——用于质因数分解 
int min_prime[MAXN];  // 每个数的最小质因子

void linear_sieve(int n) {
    for (int i = 2; i <= n; i++) {
        if (min_prime[i] == 0) {  // i 是质数
            min_prime[i] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            if (p > min_prime[i] || (long long)i * p > n) break;
            min_prime[i * p] = p;
        }
    }
}

// 用 min_prime[] 在 O(log n) 内分解质因数
vector<pair<int,int>> factorize(int n) {
    vector<pair<int,int>> factors;
    while (n > 1) {
        int p = min_prime[n], cnt = 0;
        while (n % p == 0) { n /= p; cnt++; }
        factors.push_back({p, cnt});
    }
    return factors;
}

因子个数

若 n = p₁^a₁ × p₂^a₂ × ... × pₖ^aₖ，则因子个数为 (a₁+1)(a₂+1)...(aₖ+1)。

int count_divisors(int n) {
    auto factors = factorize(n);
    int cnt = 1;
    for (auto [p, e] : factors)
        cnt *= (e + 1);
    return cnt;
}

8.5.7 欧拉 φ 函数

φ(n)（欧拉 φ 函数）统计 [1, n] 中与 n 互质（即 gcd(k, n) = 1）的整数个数。

φ(1) = 1
φ(2) = 1  （只有 1 与 2 互质）
φ(6) = 2  （1 和 5 与 6 互质）
φ(12) = 4 （1, 5, 7, 11）
φ(p) = p-1，对任意质数 p（1..p-1 均与 p 互质）

计算公式

若 n = p₁^a₁ × p₂^a₂ × ... × pₖ^aₖ，则：

φ(n) = n × (1 - 1/p₁) × (1 - 1/p₂) × ... × (1 - 1/pₖ)

单值计算

📄 查看代码：单值计算 
int euler_phi(int n) {
    int result = n;
    for (int p = 2; (long long)p * p <= n; p++) {
        if (n % p == 0) {
            while (n % p == 0) n /= p;  // 去掉所有 p 因子
            result -= result / p;        // result *= (1 - 1/p)
        }
    }
    if (n > 1) result -= result / n;     // n 本身是剩余的质因子
    return result;
}

筛求 φ(1..N)——O(N log log N)

📄 查看代码：筛求 φ(1..N)——O(N log log N) 
const int MAXN = 1000001;
int phi[MAXN];

void phi_sieve(int n) {
    // 初始化 phi[i] = i（乘法单位元步骤）
    iota(phi, phi + n + 1, 0);

    for (int p = 2; p <= n; p++) {
        if (phi[p] == p) {   // p 是质数（尚未被修改）
            for (int j = p; j <= n; j += p) {
                phi[j] -= phi[j] / p;  // phi[j] *= (1 - 1/p)
            }
        }
    }
}
// 调用 phi_sieve(n) 后，phi[i] = φ(i) 对所有 i ∈ [1, n]

φ 函数在 USACO/组合数学中的应用

  1. 费马小定理的推广： 对满足 gcd(a, n) = 1 的任意 a：a^φ(n) ≡ 1 (mod n)。这就是欧拉定理。

  2. 原根/乘法阶： a 的阶整除 φ(n)。

  3. 项链计数（Burnside）： 公式中用到 N 的每个因子 d 对应的 φ(d)。

  4. φ 的求和： Σ_{d|n} φ(d) = n。在对因子做容斥时非常有用。

// 示例：统计满足 1<=a<=b<=n 且 gcd(a,b)=1 的对数
// 答案 = 1 + Σ_{i=2}^{n} φ(i)   （"+1" 对应 (1,1)）
phi_sieve(n);
long long count = 1;
for (int i = 2; i <= n; i++)
    count += phi[i];

8.5.8 中国剩余定理（CRT）

中国剩余定理指出：若有两两互质模数的同余方程组：

x ≡ r₁ (mod m₁)
x ≡ r₂ (mod m₂)
...
x ≡ rₖ (mod mₖ)

则在模 M = m₁ × m₂ × ... × mₖ 的意义下存在唯一解 x。

两方程 CRT

对 x ≡ r₁ (mod m₁) 和 x ≡ r₂ (mod m₂)（其中 gcd(m₁, m₂) = 1）：

📄 对 `x ≡ r₁ (mod m₁)` 和 `x ≡ r₂ (mod m₂)`（其中 gcd(m₁, m₂) = 1）： 
// 返回 x 使得 x ≡ r1 (mod m1) 且 x ≡ r2 (mod m2)
// 要求 gcd(m1, m2) = 1
// 解在 mod (m1 * m2) 意义下唯一
long long crt(long long r1, long long m1, long long r2, long long m2) {
    // x = r1 + m1 * k，对某个 k
    // r1 + m1 * k ≡ r2 (mod m2)
    // m1 * k ≡ r2 - r1 (mod m2)
    // k ≡ (r2 - r1) * inv(m1) (mod m2)
    long long k = (r2 - r1 % m2 + m2) % m2 * mod_inv(m1 % m2, m2) % m2;
    return r1 + m1 * k;
    // 结果在 [0, m1*m2) 范围内，若 m1*m2 > 10^18 可能溢出
    // 必要时使用 __int128
}

广义 CRT（模数非互质）

模数非互质时，解可能不存在，用扩展 GCD 判断：

📄 模数**非互质**时，解可能不存在，用扩展 GCD 判断： 
// 返回 {x, lcm(m1,m2)} 使得 x ≡ r1 (mod m1) 且 x ≡ r2 (mod m2)
// 无解时返回 {-1, -1}
// 即使 gcd(m1, m2) > 1 也适用
pair<long long, long long> crt_general(long long r1, long long m1, long long r2, long long m2) {
    long long g = __gcd(m1, m2);
    if ((r2 - r1) % g != 0) return {-1, -1};  // 无解

    long long lcm = m1 / g * m2;
    long long diff = (r2 - r1) / g;
    long long m2g = m2 / g;

    // k ≡ diff * inv(m1/g) (mod m2/g)
    long long k = diff % m2g * mod_inv(m1 / g % m2g, m2g) % m2g;
    long long x = (r1 + m1 * k) % lcm;
    if (x < 0) x += lcm;
    return {x, lcm};
}

多方程 CRT（迭代求解）

📄 查看代码：多方程 CRT（迭代求解） 
// 求解方程组：x ≡ r[i] (mod m[i])，i = 0..k-1
// 返回 {x, M}，其中 M = 所有模数的 lcm
// 无解时返回 {-1, -1}
pair<long long, long long> crt_multi(vector<long long>& r, vector<long long>& m) {
    long long cur_r = r[0], cur_m = m[0];
    for (int i = 1; i < (int)r.size(); i++) {
        auto [x, M] = crt_general(cur_r, cur_m, r[i], m[i]);
        if (x == -1) return {-1, -1};
        cur_r = x;
        cur_m = M;
    }
    return {cur_r, cur_m};
}

USACO CRT 题型模式

"分别每 A₁ 步、B₁ 天、C₁ 小时发生一次的事件，何时三者同时发生？"

x ≡ r₁ (mod A₁)
x ≡ r₂ (mod B₁)
x ≡ r₃ (mod C₁)
→ 用 crt_multi 迭代求解

8.5.9 USACO Gold 数学题型模式

模式 1：用 DP 计数

"统计长度为 N 的合法序列数，每个元素从 1..M 中选取并满足约束。"

建模为 DP：dp[i][状态] = 以某个状态结尾的长度为 i 的序列数。答案通常需要取模。

模式 2：整除约束

"1 到 N 中有多少数至少能被 {a₁, a₂, ..., aₖ} 之一整除？"

容斥：Σ|aᵢ 的倍数| − Σ|lcm(aᵢ, aⱼ) 的倍数| + ...

// 统计 n 以内 m 的倍数：
int count_multiples(long long n, long long m) {
    return n / m;
}

模式 3：隔板法（Stars and Bars）

"将 N 个相同小球分配到 K 个盒子，满足某些约束。"

无约束：C(N+K-1, K-1)。有"每盒至多 X 个"的约束：使用容斥。

模式 4：对称性 / Burnside 引理

"统计不同的项链/着色方案数（考虑旋转/翻转的等价性）。"

Burnside 引理：所有群元素的不动点数的平均值。在 USACO 中不常见但印象深刻。

💡 思路陷阱

陷阱 1：对非素数模数使用费马小定理求逆元

错误判断： "求 a 的逆元就是 power(a, MOD-2, MOD)"
实际情况： 费马小定理要求 MOD 是素数，且 gcd(a, MOD)=1；若 MOD 不是素数（如 MOD=10⁶）则结果错误

// 错误：MOD = 10^6（不是素数）
long long inv = power(6, 1e6 - 2, 1e6);  // 6^(10^6-2) mod 10^6 ≠ 6⁻¹

// 正确：MOD 不是素数时，用扩展欧几里得
long long inv = mod_inv_general(6, 1000000);  // ext_gcd 方法
// 大多数 USACO 题的 MOD=10⁹+7（素数），直接用 Fermat 没问题

识别信号： 题目给的模数不是 10⁹+7 或 998244353 → 先验证是否为素数再决定求逆方法

陷阱 2：容斥原理的符号方向搞反

错误判断： "偶数大小子集加，奇数大小子集减"（或者反过来）
实际情况： 容斥公式：|A₁∪...∪Aₙ| = Σ|单集合| - Σ|二元交| + Σ|三元交| - ...
奇数大小子集加，偶数大小子集减（包含-排除交替）

📄 奇数大小子集**加**，偶数大小子集**减**（包含-排除交替） 
// 常见题：N 个元素，统计"至少满足 k 个条件之一"
// 错误：把 + - 方向搞反
long long ans = 0;
for (int mask = 1; mask < (1<<k); mask++) {
    int bits = __builtin_popcount(mask);
    long long term = compute_intersection(mask);
    if (bits % 2 == 0) ans += term;  // ← 错误！偶数应该减
    else               ans -= term;  // ← 错误！奇数应该加
}

// 正确
    if (bits % 2 == 1) ans += term;  // 奇数大小交集：加
    else               ans -= term;  // 偶数大小交集：减

识别信号： 容斥答案出现负数或明显偏大 → 检查 +/- 符号是否与 popcount % 2 对应正确

陷阱 3：C(n, k) 中 k < 0 或 k > n 时未做边界检查

错误判断： "公式直接套，fact[n] * inv_fact[k] * inv_fact[n-k]"
实际情况： 当 k < 0 或 k > n 时，inv_fact 数组越界或数学上 C(n,k) 应为 0

📄 C++ 完整代码 
// 错误：没有边界检查
long long C(int n, int k) {
    return fact[n] * inv_fact[k] % MOD * inv_fact[n-k] % MOD;
    // 若 k=-1 或 k=n+1，访问 inv_fact[-1] 是野指针访问
}

// 正确：加边界保护
long long C(int n, int k) {
    if (k < 0 || k > n || n < 0) return 0;  // ← 必须有！
    return fact[n] * inv_fact[k] % MOD * inv_fact[n-k] % MOD;
}

识别信号： 组合数计算在某些特殊输入下崩溃或返回极大值 → 检查边界条件

⚠️ 常见错误

  1. a * b % MOD 中的整数溢出： 若 a, b ≈ 10⁹，则 a * b 可能溢出 int 甚至 long long。务必先转换：(long long)a * b % MOD。

  2. 减法结果为负： (a - b) % MOD 在 C++ 中可能为负数。始终写成 (a - b + MOD) % MOD。

  3. inv_fact[0] = 1： 确保 inv_fact[0] = 1（因为 0! = 1）。precompute_factorials 中的倒序循环会处理此问题。

  4. C(n, k) 当 k > n 或 k < 0 时为 0： 始终检查这些边界情况。

  5. MOD 不是素数： 费马小定理要求 p 为素数。若题目使用非素数模数（罕见），用 ext_gcd 求模逆元。

  6. 大 n 下的 Lucas 定理： 当 n 极大（10¹²+）但素数模数 p 较小（< 10⁶）时，使用 Lucas 定理：C(n, k) mod p = C(n mod p, k mod p) × C(n/p, k/p) mod p。在 USACO Gold 中罕见，但在 Platinum 中会出现。

📋 章节小结

📌 核心要点

概念说明
模逆元a⁻¹ mod p = a^(p-2) mod p（费马，p 为素数）；O(log p)
阶乘表预处理 fact[], inv_fact[] 到 10⁶；O(N) 空间
C(n, k) mod pfact[n] * inv_fact[k] * inv_fact[n-k] mod p；每次查询 O(1)
容斥原理对约束的子集交替求和
筛法O(N log log N) 找出 N 以内所有质数；O(log N) 分解质因数
卡特兰数C(2n,n)/(n+1)；统计二叉树、合法括号序列数

❓ 常见问题

Q：USACO 中最常见的模数是什么？
A：10⁹+7（1,000,000,007），是素数。偶尔用 998,244,353（也是素数，用于 NTT）。

Q：怎么判断一道题需要组合数学还是 DP？
A：若问题有"漂亮的"封闭形式答案（如 C(n,k)），用组合数学。若约束有复杂依赖，可能需要 DP。通常两者结合：先用 DP 建立表格，再用组合数学求和。

Q：GCD 是什么？什么时候需要用它？
A：gcd(a, b) = 最大公因数。C++ 中用 __gcd(a, b)。用途：化简分数、检验整除性、计算 lcm = a*b/gcd(a,b)。

Q：什么时候用 Lucas 定理？
A：当 n 极大（10¹²+）但素数模数 p 较小（< 10⁶）时。在 USACO Gold 中罕见，但在 Platinum 中出现。

🔗 与后续章节的联系

  • 附录 E（数学基础）： 本章是附录 E 的延伸——模运算、质数和组合数学都在那里介绍。
  • 第 6.3 章（进阶 DP）： 数位 DP（统计满足数位约束的整数）结合了数论与 DP。
  • 第 8.2 章（DAG DP）： DAG 中的路径计数通常需要对结果取模 p。

🏋️ 练习题

🟢 简单

8.5-E1. 模意义快速幂
给定 a, n, p（p 为素数，n ≤ 10¹⁸），计算 a^n mod p。

提示

二进制取幂（power(a, n, p)）。注意溢出：若 a, p ≈ 10¹⁸，使用 __int128 或谨慎处理乘法。

8.5-E2. 网格路径计数
统计在 N×M 网格中从 (0,0) 到 (n,m) 的单调路径数（只能向右或向下）。输出 mod 10⁹+7 的结果。

提示

答案为 C(n+m, n) = (n+m)! / (n! × m!)。用预处理的阶乘和模逆元计算。

🟡 中等

8.5-M1. 序列计数 （USACO 风格）
统计长度为 N 的序列数，每个元素取自 {1, 2, ..., M}，且 K 个"特殊值"都至少出现一次。输出 mod 10⁹+7 的结果。

提示

容斥：使用 8.5.5 节中的 count_surjective(N, M, K)。枚举 j 个被排除的值。

8.5-M2. 因子和
给定 N 个数 a₁, a₂, ..., aₙ，对每个 aᵢ 输出其所有因子之和对 10⁹+7 取模的结果。

提示

用线性筛预计算最小质因子并分解每个 aᵢ。若 aᵢ = p₁^e₁ × ... × pₖ^eₖ，因子和 = 各 (1 + pᵢ + pᵢ² + ... + pᵢ^eᵢ) 的乘积 = 各 (pᵢ^(eᵢ+1) - 1) / (pᵢ - 1) 的乘积。用模逆元处理除法。

🔴 困难

8.5-H1. 项链计数 （Burnside 引理）
统计由 N 颗珠子组成、每颗用 K 种颜色之一着色的不同项链数（旋转等价的视为相同）。

提示

Burnside 引理：答案 = (1/N) × Σ_{d|N} φ(N/d) × K^d，其中 φ 为欧拉 φ 函数，求和对 N 的因子 d 进行。需要 GCD、模逆元和欧拉 φ 函数的计算。

🏆 挑战

8.5-C1. 树上期望值 （接近 USACO Platinum 难度）
给定 N 个顶点的树，每个顶点初始无色。以 1/2 的概率将每个顶点染成红色，以 1/2 的概率染成蓝色。求两端颜色相同的边的期望数。将答案表示为最简分数 p/q，输出 p × q⁻¹ mod 10⁹+7。

提示

由期望的线性性：E[同色边数] = 边数 × P(两端同色)。对每条边，P(同色) = 1/4（均红）+ 1/4（均蓝）= 1/2。所以答案为 (N-1)/2。输出 (N-1) × mod_inv(2) % MOD。

（有趣的变体是给每个顶点不同的着色概率——用相同思路推广即可。）