浅谈二项式反演

· · 算法·理论

摘自 OI 中的数学基础第 11 章。

11. 二项式反演

11.1 原理

设 f(n) 表示恰好使用 n 个不同元素形成特定结构的方案数，g(n) 表示从这 n 个不同元素中选出若干个元素形成特定结构的方案数。

已知 f(n) 求 g(n) 是简单的，枚举选出多少个元素有

g(n)=\sum_{i=0}^{n} C_n^i f(i)

反着做却是困难的，这个过程就叫二项式反演。有公式：

f(n)=\sum_{i=0}^{n} C_n^i (-1)^{n-i}g(i)

这是二项式反演的形式之一。二项式反演的作用就是把“恰好”转化为“钦定”。

*11.2 证明

我们先给出一个引理

C_n^r C_r^k=C^k_n C^{r-k}_{n-k}

组合意义和代数推导都易证。

把上述式子展开：

\begin{aligned} f(n)&=\sum_{i=0}^{n} C_n^i (-1)^{n-i} [\sum_{j=0}^{i} C_i^j f(j)]\\ &=\sum_{i=0}^{n} \sum_{j=0}^{i} C_n^i C_i^j (-1)^{n-i} f(j) \\ &= \sum_{j=0}^{n} \sum_{i=j}^{n} C_n^i C_i^j (-1)^{n-i} f(j) \\ &= \sum_{j=0}^{n} [f(j) \times \sum_{i=j}^{n} C_n^i C_i^j (-1)^{n-i} ] \end{aligned}

根据引理可得

\begin{aligned} f(n)&=\sum_{j=0}^{n} [f(j) \times \sum_{i=j}^{n} C_n^j C_{n-j}^{i-j} (-1)^{n-i} ]\\ &=\sum_{j=0}^{n} [C_n^j f(j) \times \sum_{i=j}^{n}C_{n-j}^{i-j} (-1)^{n-i} ]\\ \end{aligned}

作换元，令 k=i-j，则 i=k+j，即有

\begin{aligned} f(n)&=\sum_{j=0}^{n} [C_n^j f(j) \times \sum_{k=0}^{n-j}C_{n-j}^{k} (-1)^{n-j-k} ]\\ &=\sum_{j=0}^{n} [C_n^j f(j) \times \sum_{k=0}^{n-j}C_{n-j}^{k} (-1)^{n-j-k}1^k ] \end{aligned}

由组合数性质 6 可得当且仅当 n=j 时 \sum_{k=0}^{n-j}C_{n-j}^{k} (-1)^{n-j-k}1^k 为 1，故：

f(n)=f(n)

上述变换都是等价变换，证毕。

11.3 常用形式

下面我们给出二项式反演的全部形式：

形式一：

f(x)=\sum_{i=0}^x (-1)^iC_x^i g(n) \Leftrightarrow g(x)=\sum_{i=0}^x (-1)^i C_x^i f(i)

这是二项式反演的基本公式。

形式二：

f(x)=\sum_{i=x}^n C_i^x g(i) \Leftrightarrow g(x)=\sum_{i=x}^n (-1)^{i-x} C_i^x f(i)

这是钦定意义下的“至少”转“恰好”。

形式三：

f(x)=\sum_{i=0}^x C_{n-i}^{n-x} g(i) \Leftrightarrow g(x)=\sum_{i=0}^x (-1)^{x-i} C_{n-i}^{n-x} f(i)

这是钦定意义下的“至多”转“恰好”。

形式四：

f(x)=\sum_{i=0}^x C_{i}^{x} g(i) \Leftrightarrow g(x)=\sum_{i=0}^x (-1)^{x-i} C_{i}^{x} f(i)

这是选择意义下的“至多”转“恰好”。

形式五：

f(x)=\sum_{i=x}^n C_{n-x}^{n-i} g(i) \Leftrightarrow g(x)=\sum_{i=x}^n (-1)^{i-x} C_{n-x}^{n-i} f(i)

这是选择意义下的“至少”转“恰好”。

二项式反演的本质是系数特殊的容斥原理。

11.4 例题

二项式反演的题目中，最常用的是形式二。同时有些题转成“钦定”以后，还要配合计数 DP 求解。

在二项式反演时，首先要找到“恰好”，然后分清“钦定”和“选择”的区别，使用恰当的反演形式。

AT_abc423_f [ABC423F] Loud Cicada

设 g(x) 表示钦定至少 x 种蝉爆发的方案数。状压枚举集合 S，则

g(x)=\sum _{\operatorname{popcount}(S) = x} \lfloor \dfrac{Y}{\operatorname{lcm}_{i \in S} a_i} \rfloor

解释一下，对于集合 S 中的所有元素 i，求出 a_i 的最小公倍数，则集合 S 中的所有蝉爆发的充要条件就是年份为最小公倍数的倍数。在 Y 以内共有 \lfloor \dfrac{Y}{\operatorname{lcm}_{i \in S} a_i} \rfloor 个年份。

由二项式反演形式二得

f(m)=\sum_{i=m}^{n} (-1)^{i-m}C_m^i g(i)

复杂度 O(n2^n)，杨辉三角 O(n^2) 预处理组合数即可。注意 \operatorname{lcm} 最好开个 __int128。

#define popcount __builtin_popcount
const int N = 25;
int n, m;
long long y, a[N], g[N], c[N][N];

void _main() {
    cin >> n >> m >> y;
    for (int i = 1; i <= n; i++) cin >> a[i];
    for (int i = 0; i <= n; i++) c[i][0] = 1, c[i][i] = 1;
    for (int i = 0; i <= n; i++) {
        for (int j = 1; j < i; j++) c[i][j] = c[i - 1][j] + c[i - 1][j - 1];
    }
    for (int s = 0; s < (1 << n); s++) {
        __int128 x = 1;
        for (int i = 1; i <= n; i++) {
            if (s >> (i - 1) & 1) x = x / __gcd<__int128>(x, a[i]) * a[i];
            if (x > y) break;  // 注意这句
        } g[popcount(s)] += y / x;  
    }
    long long res = 0;
    for (int i = m; i <= n; i++) {
        if ((m - i) & 1) res -= c[i][m] * g[i];
        else res += c[i][m] * g[i];
    } cout << res;
}

P10596 BZOJ2839 集合计数

看到不好求的“恰好”，考虑二项式反演。设 f(i) 表示交集大小恰好为 i 的方案数，g(i) 表示钦定至少 i 个元素属于交集的方案数。

那么对于 g(i)，先从 n 个元素中选出 i 个，方案数 C_n^i。剩下 n-i 个元素可选可不选，就是求大小为 2^{n-i} 的集合的非空子集数，答案为 2^{2^{n-i}}-1，乘法原理：

g(i)=C_n^i (2^{2^{n-i}}-1)

由二项式反演形式二得

f(k)=\sum_{i=k}^{n}(-1)^{i-k} C_i^k g(i)

复杂度 O(n \log n)。注意求 g(i) 要用欧拉定理。

const int N = 1e6 + 5;
mint fac[N], ifac[N], g[N];
mint C(int n, int m) {
    return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];
}
int n, k;
void _main() {
    cin >> n >> k;
    fac[0] = ifac[0] = 1;
    for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
    for (int i = 0; i <= n; i++) g[i] = C(n, i) * (mint(2).pow(mod32<1000000006>(2).pow(n - i).val) - 1);
    mint res = 0;
    for (int i = k; i <= n; i++) {
        if ((i - k) & 1) res -= C(i, k) * g[i];
        else res += C(i, k) * g[i];
    } cout << res;
}

P6521 [CEOI 2010] pin (day2)

首先做一个等价变形，把不同变成相同，令 m \gets 4-m 即可。

将 i 位相同相同视为一个限制，共有 4 种限制，求的是恰好满足 m 个限制的方案数。

套路地，钦定至少 x 个限制必须满足，设其为 S，对于选出的两个串使用字符串哈希，且仅对 S 中的位置哈希，这样就能求出 x 位必须相同的方案数。根据二项式反演形式二求得答案即可。

const int N = 5e4 + 5;
int n, m;
long long g[5];
char a[N][5];
const int fac[] = {1, 1, 2, 6, 24};
int C(int n, int m) {return n < m ? 0 : fac[n] / fac[m] / fac[n - m];}

void _main() {
    cin >> n >> m, m = 4 - m;
    for (int i = 1; i <= n; i++) cin >> (a[i] + 1);
    for (int s = 0; s < (1 << 4); s++) {
        umap<uint64_t, int> mp;
        for (int i = 1; i <= n; i++) {
            uint64_t h = 0;
            for (int j = 1; j <= 4; j++) {
                h *= 13331;
                if (s >> (j - 1) & 1) h += a[i][j];
            }
            g[popcount(s)] += mp[h], mp[h]++;
        }
    }
    long long res = 0;
    for (int k = m; k <= 4; k++) {
        if ((k - m) & 1) res -= C(k, m) * g[k];
        else res += C(k, m) * g[k];
    } cout << res;
}

P6076 [JSOI2015] 染色问题

将“每一”变成“恰好”，用二项式反演的思维来思考。

先处理颜色的限制，不妨计算选择至多出现 x 种颜色的方案数。这样每个格子要不然不染色，要不然染 x 种种的一种。

类似地处理掉行列限制后：问题变成：一个 n\times m 的棋盘，至多 n 行被染色，至多 m 行被染色，至多出现 c 种颜色，基础乘法原理得到 (c+1)^{nm}。

重复应用三次二项式反演形式四，得到答案为

\sum_{k=0}^c (-1)^{c-k}C_c^k\sum_{j=0}^m (-1)^{m-j} C_m^j \sum_{i=0}^n (-1)^{n-i} C_n^i (k+1)^{ij}

复杂度 O(nmc)，可以通过。

const int N = 405;
mint fac[N], ifac[N];
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];}

int n, m, c;
void _main() {
    fac[0] = ifac[0] = 1;
    for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
    cin >> n >> m >> c;
    mint res = 0;
    for (int k = 0; k <= c; k++) {
        mint a1 = 0;
        for (int j = 0; j <= m; j++) {
            mint a2 = 0;
            for (int i = 0; i <= n; i++) {
                if ((n - i) & 1) a2 -= C(n, i) * mint(k + 1).pow(i * j);
                else a2 += C(n, i) * mint(k + 1).pow(i * j);
            } 
            if ((m - j) & 1) a1 -= C(m, j) * a2;
            else a1 += C(m, j) * a2;
        }
        if ((c - k) & 1) res -= C(c, k) * a1;
        else res += C(c, k) * a1;
    } cout << res;
}

P5505 [JSOI2011] 分特产

正难则反，记 f(x) 为钦定至少 x 人没有分到的方案数，由二项式反演形式二得到答案为

\sum_{i=0}^{n-1} (-1)^i C_n^if(i)

需要求出 f(x)。考虑插板法，第 k 个特产没有分到的方案是 C_{a_k+n-x-1}^{n-x-1}，根据乘法原理合并：

f(x)=\prod_{k=1}^m C_{a_k+n-x-1}^{n-x-1}

逆元法预处理组合数即可，复杂度 O(nm)。

const int N = 2005;
mint fac[N], ifac[N];
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m]; }
int n, m, a[N];
mint f(int x) {
    mint res = 1;
    for (int k = 1; k <= m; k++) res *= C(a[k] + n - x - 1, n - x - 1);
    return res;
}

void _main() {
    fac[0] = ifac[0] = 1;
    for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
    cin >> n >> m;
    for (int i = 1; i <= m; i++) cin >> a[i];
    mint res = 0;
    for (int i = 0; i < n; i++) {
        if (i & 1) res -= C(n, i) * f(i);
        else res += C(n, i) * f(i);
    } cout << res;
}

CF1228E Another Filling the Grid

我们给出二项式反演形式二的二维形式：

f(i,j)=\sum_{x=i}^n \sum_{y=j}^n C_x^i C_y^j g(x,y) \Leftrightarrow g(i,j)=\sum_{x=i}^n \sum_{y=j}^n (-1)^{x-i+y-j} C_x^i C_y^j f(x,y)

特判 k=1。设 f(i,j) 表示钦定至少 i 行 j 列不满足要求的方案数，g(i,j) 表示恰好 i 行 j 列不满足要求的方案数。容易发现 f,g 满足上述公式。只需求出 f(i,j)。基础乘法原理得到

f(i,j)=(k-1)^{n^2-(n-i)(n-j)} k^{(n-i)(n-j)}

直接做即可。复杂度 O(n^2)。

const int N = 255;
int n, k;
mint fac[N], ifac[N];
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];}

mint f(int x, int y) {
    return mint(k - 1).pow(n * n - (n - x) * (n - y)) * mint(k).pow((n - x) * (n - y));
} 

void _main() {
    cin >> n >> k;
    if (k == 1) return cout << 1, void();
    fac[0] = ifac[0] = 1;
    for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
    mint res = 0;
    for (int i = 0; i <= n; i++) {
        for (int j = 0; j <= n; j++) {
            if ((i + j) & 1) res -= C(n, i) * C(n, j) * f(i, j);
            else res += C(n, i) * C(n, j) * f(i, j);
        }
    } cout << res;
}

AT_abc172_e [ABC172E] NEQ

其实可以用第 12 章错排的方法解决，放在这里有点大炮打蚊子了。

设 g(x) 表示钦定 x 位相同的方案数，f(x) 为钦定至少 x 位相同的方案数，根据二项式反演形式二

g(x)=\sum_{i=x}^n (-1)^{i-x}C_x^i f(i)

答案即为

g(0)=\sum_{i=0}^n (-1)^n f(i)

此时二项式反演退化为普通容斥原理。容易得到

f(x)=C_n^x A_m^x (A_{n-x}^{m-x})^2

复杂度 O(n)。

const int N = 5e5 + 5;
int n, m;
mint fac[N], ifac[N];
mint A(int n, int m) {return n < m ? 0 : fac[n] * ifac[n - m];}
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];}
mint f(int x) {
    return C(n, x) * A(m, x) * A(m - x, n - x) * A(m - x, n - x);
}

void _main() {
    cin >> n >> m;
    fac[0] = ifac[0] = 1;
    for (int i = 1; i <= max(n, m); i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
    mint res = 0;
    for (int i = 0; i <= n; i++) {
        if (i & 1) res -= f(i);
        else res += f(i);
    } cout << res;
}

*P4859 已经没有什么好害怕的了

形式化题意：给出长为 n 的序列 a,b，将其两两配对使得 a_i>b_i 的数目减去 a_i<b_i 的数目恰为 k，求方案数。

设 a_i>b_i 的数量为 m，则 m-(n-m)=k，解得 m=\dfrac{n+k}{2}。若 2 \nmid (n+k)，直接判定无解。此时问题变成 a_i>b_i 的数目恰好为 m。使用二项式反演化恰好为钦定。设 g(i) 为钦定至少 a_i>b_i 的组数不小于 i 的方案数，根据二项式反演所求为

\sum_{i=m}^{n} (-1)^{i-m} C_{i}^k g(i)

现在的问题就是求 g(i)，这玩意很 DP，所以设 dp_{i,j} 表示前 i 个数中选出 j 组 a_i>b_i 的方案数。DP 题套路对 a,b 排序，然后可以推出转移方程

dp_{i,j}=dp_{i-1,j}+(last_i-j+1)dp_{i-1,j-1}

可以看看下文的排列计数 DP，从插入角度考虑 (a_i,b_i) 的贡献，一种情况是不作为新的 a_i>b_i，还有一种就是考虑有多少个取代位置。设 last_i 表示 b 序列中第一个小于 a_i 的数的下标，则插入位置就有 last_i-j+1 个，加法原理合并答案。

于是我们 O(n^2) 求出 dp_{i,j}。考察 g(i) 的定义可得

g(i)=dp_{n,i}\times (n-i)!

就是在有序的 DP 基础上考虑剩下 n-i 个全排列。

二项式反演中组合数可以逆元法来求，总复杂度 O(n^2)。

const int N = 2005;
int n, m, k, a[N], b[N], last[N];
mint fac[N], ifac[N], dp[N][N], g[N];
mint C(int n, int m) {
    if (n < m) return 0;
    return fac[n] * ifac[m] * ifac[n - m]; 
}

void _main() {
    cin >> n >> k;
    if ((n + k) & 1) return cout << 0, void();
    fac[0] = ifac[0] = 1;
    for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];

    m = (n + k) / 2;
    for (int i = 1; i <= n; i++) cin >> a[i];
    for (int i = 1; i <= n; i++) cin >> b[i];
    sort(a + 1, a + n + 1), sort(b + 1, b + n + 1);
    for (int i = 1; i <= n; i++) last[i] = lower_bound(b + 1, b + n + 1, a[i]) - b - 1;
    dp[0][0] = 1;
    for (int i = 1; i <= n; i++) {
        dp[i][0] = dp[i - 1][0];
        for (int j = 1; j <= i; j++) {
            dp[i][j] = dp[i - 1][j] + dp[i - 1][j - 1] * max(0, last[i] - j + 1);
        }
    }
    for (int i = 0; i <= n; i++) g[i] = dp[n][i] * fac[n - i];
    mint res = 0;
    for (int i = m; i <= n; i++) {
        if ((i - m) & 1) res -= C(i, m) * g[i];
        else res += C(i, m) * g[i];
    } cout << res;
}

P10597 BZOJ4665 小 w 的喜糖

和 P4859 很像。设 f(x) 表示恰好 x 个位置不同的方案数，g(x) 表示钦定至少 x 个位置的不同方案数，由二项式反演形式二得

f(x)=\sum_{i=x}^n (-1)^{i-x} C_x^i g(i)

所求为

f(0)=\sum_{i=0}^n (-1)^i g(i)

我们再一次看到了二项式反演退化为普通容斥原理。考虑求 g(x)。

记第 i 种糖果的个数有 c_i 个，先假设同种糖果不同，根据多重集排列数，除掉 c_i! 即可。

设 dp_{i,j} 表示前 i 种糖果中至少 j 个人的糖果与原来相同的方案数，枚举 k 颗被钦定的第 i 种糖果，容易得到

dp_{i,j}=\sum_{k=0}^{\min(c_i,j)} \dfrac{dp_{i-1,j-k} \times C_{c_i}^k}{(c_i-k)!}

则 g(x)=dp_{n,x} \times (n-x)!。复杂度 O(n^3)。

const int N = 2005;  
mint fac[N], ifac[N], dp[N][N];  
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];}
int n, a[N], c[N];

void _main() {
    fac[0] = ifac[0] = 1;
    for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
    cin >> n; 
    for (int i = 1; i <= n; i++) cin >> a[i], c[a[i]]++;
    dp[0][0] = 1;
    for (int i = 1; i <= n; i++) {
        for (int j = 0; j <= n; j++) {
            for (int k = 0; k <= min(c[i], j); k++) dp[i][j] += dp[i - 1][j - k] * C(c[i], k) * ifac[c[i] - k];
        }
    }
    mint res = 0;
    for (int i = 0; i <= n; i++) {
        if (i & 1) res -= dp[n][i] * fac[n - i];
        else res += dp[n][i] * fac[n - i];
    } cout << res;
}

*P6478 [NOI Online #2 提高组] 游戏

“恰好 k 次非平局”很难处理，不妨转为“钦定 k 次非平局”。设恰好 k 次非平局的方案数为 g(k)，钦定至少 k 次非平局的方案数为 f(k)，显然有 f(n)=\sum_{i=n}^m C_i^n g(i)，根据二项式反演形式二得

g(n)=\sum_{i=n}^m (-1)^{i-n} C_i^n f(i)

考虑求 f(i) 即可。考虑树形 DP，设 dp_{u,i} 表示在以 u 为根的子树中钦定 i 个点且必须有胜负的方案数。

考虑对于儿子 v 枚举在 v 中选择 j 个点，有转移：

dp_{u,i}=\sum_{(u,v)} \sum_{j =0}^i dp_{v,j}dp_{u,i-j}

还可以选择一个点与 u 配对，于是：

dp_{u,i} \gets dp_{u,i}+dp_{u,i-1}+cnt_{u,x \oplus 1}

其中 x 为点 u 属于哪位玩家，cnt_{u,0/1} 表示以 u 为根的子树中有多少个属于 0/1 玩家。

这是一个树形背包，且应该跑 01 背包。根据树形背包的复杂度证明，复杂度是严格 O(n^2) 的。

这样有 g(i)=dp_{1,i}\times (m-i)!。求出来以后套进二项式反演即可。

const int N = 5005;
int n, m, u, v, belong[N];
char c;
int tot = 0, head[N];
struct Edge {
    int next, to;
} edge[N << 1];
inline void add_edge(int u, int v) {
    edge[++tot].next = head[u], edge[tot].to = v, head[u] = tot;
}
mint f[N], g[N], fac[N], ifac[N];
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];}

int sz[N], cnt[2][N];
mint dp[N][N];
void dfs(int u, int fa) {
    sz[u] = 1, cnt[belong[u]][u] = 1, dp[u][0] = 1;
    for (int j = head[u]; j != 0; j = edge[j].next) {
        int v = edge[j].to;
        if (v == fa) continue;
        dfs(v, u);
        vector<mint> f(min(sz[u] + sz[v], m), 0);
        for (int j = 0; j <= min(sz[u], m); j++) {
            for (int k = 0; k <= min(sz[v], m - j); k++) {
                f[j + k] += dp[u][j] * dp[v][k];
            }
        }
        sz[u] += sz[v], cnt[0][u] += cnt[0][v], cnt[1][u] += cnt[1][v];
        copy(f.begin(), f.end(), dp[u]);
    }
    for (int i = cnt[belong[u] ^ 1][u]; i >= 0; i--) {
        dp[u][i + 1] += dp[u][i] * (cnt[belong[u] ^ 1][u] - i);
    }
}

void _main() {
    fac[0] = ifac[0] = 1;
    for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
    cin >> n, m = n / 2;
    for (int i = 1; i <= n; i++) cin >> c, belong[i] = c ^ 48;
    for (int i = 1; i < n; i++) {
        cin >> u >> v;
        add_edge(u, v), add_edge(v, u);
    }
    dfs(1, -1);
    for (int i = 0; i <= m; i++) g[i] = dp[1][i] * fac[m - i];
    for (int i = 0; i <= m; i++) {
        for (int j = i; j <= m; j++) {
            if ((j - i) & 1) f[i] -= C(j, i) * g[j];
            else f[i] += C(j, i) * g[j];
        }
        cout << f[i] << '\n';
    }
}

*P3270 [JLOI2016] 成绩比较

前置知识：第 24 章 Lagrange 插值求自然数幂和。

设 f(k) 表示恰好 k 个人被碾压的方案数，g(k) 表示钦定至少 k 个人被碾压的方案数。由二项式反演形式二得

f(k)=\sum_{i=k}^n (-1)^{i-k} C_i^k g(i)

考虑求 g(k)。分步地，从 n-1 个人中钦定 k 个人被碾压。对于每一科 i，枚举 D 神的分数 j，选出 r_i-1 个人高于 D 神，再乘法原理做一做得到

g(k)=C_{n-1}^k \prod_{i=1}^m \sum_{j=1}^{U_i} C_{n-k-1}^{r_i-1} j^{n-r_i}(U_j-1)^{n-r_i}

复杂度 O(nmV \log n)，无法通过。

考虑消去值域限制。开始大力推柿子：

\begin{aligned} g(k)&=C_{n-1}^k \prod_{i=1}^m \sum_{j=1}^{U_i} C_{n-k-1}^{r_i-1} j^{n-r_i}(U_j-1)^{n-r_i}\\ &=C_{n-1}^k \prod_{i=1}^m \sum_{j=1}^{U_i} C_{n-k-1}^{r_i-1} j^{n-r_i}\sum_{k=0}^{r_i-1} (-1)^{r_i-k-1}C_{r_i-1}^k {U_i}^k j^{r_i-k-1}\\ &=C_{n-1}^k \prod_{i=1}^m C_{n-k-1}^{r_i-1} \sum_{k=0}^{r_i-1} (-1)^{r_i-k-1} C_{r_i-1}^{k} {U_i}^k \sum_{j=1}^{U_i} j^{n-k-1} \end{aligned}

设

h(i)=\sum_{k=0}^{r_i-1} (-1)^{r_i-k-1} C_{r_i-1}^{k} {U_i}^k \sum_{j=1}^{U_i} j^{n-k-1}

则

g(k)=C_{n-1}^k \prod_{i=1}^m C_{n-k-1}^{r_i-1} h(i)

只要处理出 \sum_{j=1}^{U_i} j^{n-k-1}，h(i) 的预处理是简单的。使用 Lagrange 插值在 O(n) 内求出自然数幂和，然后预处理 h(i)，最后套回去即可。总复杂度 O(n^2m)。

const int N = 105;
int n, m, k, u[N], r[N];
mint h[N], fac[N], ifac[N];
mint C(int n, int m) {return n < m ? 0 : fac[n] * ifac[m] * ifac[n - m];}
mint g(int k) {
    mint res = C(n - 1, k);
    for (int i = 1; i <= m; i++) res *= C(n - k - 1, r[i] - 1) * h[i];
    return res;
}
namespace Lagrange {
    mint a[N], p[N], s[N];
    mint solve(int n, int k) {
        p[0] = 1, s[k + 3] = 1;
        for (int i = 1; i <= k + 2; i++) a[i] = a[i - 1] + mint(i).pow(k), p[i] = p[i - 1] * (n - i);
        for (int i = k + 2; i >= 1; i--) s[i] = s[i + 1] * (n - i);
        if (n <= k + 2) return a[n];
        mint res = 0;
        for (int i = 1; i <= k + 2; i++) {
            mint cur = a[i] * p[i - 1] * s[i + 1] * ifac[i - 1] * ifac[k + 2 - i];
            if ((k + 2 - i) & 1) res -= cur;
            else res += cur;
        } return res;
    }
} using Lagrange::solve;

void _main() {
    fac[0] = ifac[0] = 1;
    for (int i = 1; i < N; i++) fac[i] = fac[i - 1] * i, ifac[i] = ~fac[i];
    cin >> n >> m >> k;
    for (int i = 1; i <= m; i++) cin >> u[i];
    for (int i = 1; i <= m; i++) cin >> r[i];
    for (int i = 1; i <= m; i++) {
        for (int k = 0; k < r[i]; k++) {
            mint cur = C(r[i] - 1, k) * mint(u[i]).pow(k) * solve(u[i], n - k - 1);
            if ((r[i] - k - 1) & 1) h[i] -= cur;
            else h[i] += cur;
        }
    }
    mint res = 0;
    for (int i = k; i <= n; i++) {
        if ((i - k) & 1) res -= C(i, k) * g(i);
        else res += C(i, k) * g(i);
    } cout << res;
}