多项式与生成函数基础

多项式很有工程价值（不仅限于计算机科学），研究多项式计算及其性质，可以加速许多运算。而以 FFT 为基础的算法可以直接操纵生成函数。

概念

一个环 $R$ 上的多项式写作 $f(x)$ ，代表 $\sum\limits_{i=0} f_ix^{i}$ 。其中 $f_i\in R$ 。 $x$ 被称为这个多项式的自由元，所有 $R$ 上的多项式组成多项式环 $R[x]$ 。

如果允许项数无穷的，即 $f(x)=f_0+f_1x+f_2x^2+\cdots$ ，那么可以得到形式幂级数环 $R[[x]]$ ，其中的每个元素 $f$ 称为形式幂级数。

对于一个多项式 $f(x)$ ，称其系数非零的最高次项的次数为该多项式的度（Degree），记作 $\deg F$ 。
使得 $f(x)=0$ 的所有 $x$ 称为多项式的根。
如果 $f_i$ 均为实数，那么 $f$ 称为实系数多项式；若可以为复数，则成为复系数多项式。

表示方法

形式幂级数形式表示的多项式称为系数表示法，而将 $x=x_i$ 带入多项式，求出的 $y_i=f(x_i)$ ，称 $(x_i,y_i)$ 为 $f$ 在 $x_i$ 处的点值，使用若干点值表示多项式的方式称为点值表示法。

拉格朗日插值算法可以将点值表示法转化为系数表示法。

多项式加减

有：

$f(x)+g(x)=\sum_{i=0}^n (f_i+g_i)x^i$

时间复杂度为 $O(n)$ 。

卷积

设 $\textbf{f, g}$ 是两个数列，那么这两个数列的卷积 $\textbf{h}$ 定义为：

$h_k=\sum_{i+j=k}f_ig_j$

多项式乘法

给定多项式 $f(x),g(x)$ ，那么其乘积 $q(x)$ 满足：

$q(x)=\sum_{i=0}^{n}\sum_{j=0}^m a_ib_jx^{i+j}$

这个过程可以使用 FFT 加速，将在下文介绍。

普通生成函数

生成函数（generating function），又称母函数，是一种形式幂级数（它不是函数， $x$ 不是自变量而是自由元），其每一项的系数可以提供关于这个序列的信息^[1]。这是个很有用的东西，是用于刻画数列的组合工具，能够将很多复杂的组合问题赋以简洁的代数形式，考场上也会出现不需要多项式算法只用生成函数就能解决的问题。

大多数生成函数可以表示为：

$F(x)=\sum_{i} a_i k_i(x)$

其中 $k_i(x)$ 被称之为核函数，对于普通的生成函数， $k_i(x)=x^i$ 。注意这个 $a$ 既可以是有穷序列，也可以是无穷序列。

在本小节中我们只讨论普通生成函数（OGF），它的核函数为 $k_n(x)=x^n$ 。实际上，普通生成函数的系数就是序列 $a$ 的通项公式。

术语约定

对于 OGF $F(x)$ ， $x$ 是自由元，因此没有实际意义，只是一个占位符， $x_i$ 用来代表它是第 $i$ 个占位符， $F[i]$ 可以表示 $F(x)$ 的第 $i$ 次项系数，简称 $i$ 次项， $0$ 次项称为常数项。

为了方便提取其中某一项的系数，有 $[x^i]F(x)=F[i]$ 。

生成函数中的求和符号并不是无穷项求和，它只是一个记号，代表的是序列。两个生成函数 $F(x),G(x)$ 相等当且仅当 $\forall i\in\mathbb{N}, F[i]=G[i]$ 。

对于序列 $a$ ，如果从某一项开始它后面都是 $0$ ，那么可以认为它是一个有穷序列。如果它是有穷的，那么它的生成函数是一个有穷幂级数，否则是无穷幂级数。

基本运算

考虑两个序列 $a,b$ 的普通生成函数 $F(x),G(x)$ ，那么它们的加减法和乘法运算是跟多项式一样的。因此对于生成函数 $F(x),G(x)$ ，有：

$F(x)\pm G(x)=\sum_{n}(a_n\pm b_n)x^n,\\ \begin{aligned} F(x)G(x)&=\sum_{n}x^n\sum_{i=0}^n a_i b_{n-i}\\ &=\sum_{i}\sum_j F[i]G[j]x^{ij} \end{aligned}$

因此 $F(x)\pm G(x)$ 是序列 $\left<a_n\pm b_n\right>$ 的普通生成函数， $F(x)G(x)$ 是序列 $\left<\sum_{i=0}^n a_i b_{n-i}\right>$ 的普通生成函数。

对于幂次运算， $F(x)$ 的 $n$ 次幂记作 $F^n(x)$ ，特别地， $F^0(x)=1$ 。

对 $F(x)$ 的 $n$ 次幂提取 $m$ 次项系数，可以记 $F^n[m]$ 为：

$F^n[m]=[x^m]F^n(x)$

封闭形式

形式幂级数的生成函数不一定好算，有时会将其转换为封闭形式。

比如说 $\left<1,1,\cdots\right>$ 的普通生成函数为 $F(x)=\sum_{n\ge 0}x^n$ ，发现 $F(x)x+1=F(x)$ ，因此得到 $F(x)=\cfrac{1}{1-x}$ ，这就是生成函数的封闭形式。

作为练习，可以看一下将这几个题：

写出 $a=\left<0,1,1,\cdots \right>$ 的生成函数。

有 $F(x)=\sum_{n\ge 1}x^n$ ，因此 $F(x)x+x=F(x)\Rightarrow F(x)=\cfrac{x}{1-x}$ 。

写出 $a=\left<1,0,1,0,1,0\cdots \right>$ 的生成函数。

有 $F(x)=\sum_{n\ge 0} x^{2n}=\sum_{n\ge 0} (x^{2})^n$ ，因此 $F(x)x^2+1=F(x)\Rightarrow \cfrac{1}{1-x^2}$ 。

写出序列 $a_n=\dbinom{m}{n}$ 的生成函数（ $m$ 是常数， $n\ge 0$ ）。

$\begin{aligned} F(x)&=\sum_{n\ge 0}\binom{m}{n}x^n\\ &=(1+x)^m \end{aligned}$

写出 $a=\left<1,2,3,4,5,6\cdots \right>$ 的生成函数。

$\begin{aligned} F(x)&=\sum_{n\ge 0}nx^{n-1}\\ &=\sum_{n\ge 0} (x^n)'\\ &=\frac{1}{(1-x)^2} \end{aligned}$

可以发现求一个序列的前缀和的生成函数只需要乘上 $\cfrac{1}{1-x}$ ，求差分只需要乘上 $1-x$ 。

实际上我们有：

$\frac{1}{(1-x)^{n+1}}=\sum_{i}\binom{n+i}{n}x^i$

为什么呢？当 $n=1$ 时，原式成立；当 $n>1$ 时，有：

$\begin{aligned} \frac{1}{(1-x)^{n+1}}&=\frac{1}{1-x}\frac{1}{(1-x)^{n}}\\ &=\sum_{i}x^i\sum_{i}\binom{n+i-1}{n-1}x^i\\ &=\sum_{i}x^i\sum_{j=0}^i\binom{n+j-1}{j}\\ &=\sum_{i}\binom{n+i}{n}x^i \end{aligned}$

它还可以说明一个问题：如果 $f(x)$ 是关于 $x$ 的 $k$ 次多项式，那么 $f(x)$ 的差分是 $k-1$ 次多项式，前缀和是 $k+1$ 次多项式。

简单题

说了这么多感觉很绕？我们通过几道简单的题目来见证生成函数的强大威力，即使多项式算法还没有登场。

生成函数能简单的将序列的运算转换为式子的运算，并通过系数来还原序列。对于无穷项的序列，可以试将其转化为封闭形式进行计算；当封闭形式无法直接得出结果，也可以将其转化会形式幂级数。

[六省联考 2017] 组合数问题

Portal.

令 $\displaystyle F(x)=\sum_{i} \binom{nk}{i}x^i=(1+x)^{nk}$ ，那么要求的就是 $\displaystyle\sum_{i\bmod k=r}F[i]$ 。

前者直接暴力计算生成函数的卷积，然后快速幂即可。底数是什么？千万记住， $x$ 只是自由元，因此 $1+x\Rightarrow a[0]=a[1]=1$ 。由于是求对 $k$ 取模的结果，因此幂次可以在模 $k$ 的意义下进行，最终答案就是 $ans[r]$ 。

查看代码

#include <bits/stdc++.h>
using namespace std;
typedef vector<int> VI; 
typedef long long i64; 

int n, p, k, r; 
VI operator * (const VI &a, const VI &b) {
    VI ans(k);
    for (int i = 0; i < k; ++i)
        for (int j = 0; j < k; ++j)
            ans[(i + j) % k] = (ans[(i + j) % k] + 1ll * a[i] * b[j]) % p; 
    return ans; 
}

int main(void) {
    cin >> n >> p >> k >> r; 
    vector<int> a(k), ans(k); ans[0] = 1;
    if (k == 1) a[0] = 2; else a[0] = a[1] = 1; 
    i64 e = 1ll * n * k;
    for (; e; e >>= 1, a = a * a) if (e & 1) ans = ans * a; 
    cout << ans[r] << "\n";
    return 0;
}

[Luogu P2000] 拯救世界

Portal.

对于一种召唤方式，设 $a_n$ 代表这种召唤方式选择 $n$ 个的方案数，求出序列 $a$ 的生成函数。那么两种召唤方式一共有 $n$ 块石头的选择方式就是这两个生成函数的卷积 $F$ 的第 $n$ 项系数。

一共有 $10$ 个限制条件，依次写出它们的生成函数（封闭形式），然后把它们都乘起来得到 $\cfrac{1}{(1-x)^5}$ 。转化成形式幂级数就是 $\displaystyle\sum_{i}\binom{i+4}{4}x^i$ 。因此答案是 $\dbinom{n+4}{4}$ 。

查看代码

# decimal 采用压位 NTT 实现，足够快
from decimal import *
getcontext().prec = 1000000
n = Decimal(input())
print((n + 1) * (n + 2) * (n + 3) * (n + 4) // 24)

拉格朗日插值

拉格朗日插值是沟通多项式的系数形式与点值形式的重要公式。

事实上， $n$ 个点确定唯一的多项式，其最高项次数为 $n-1$ 。

从系数表示法转化为点值表示法可以考虑直接带入，时间复杂度为 $O(n^2)$ 。可以使用多项式算法加速到 $O(n\log^2 n)$ 。

可以使用拉格朗日插值法将点值表示法转化为系数表示法，通用情况时时间复杂度为 $O(n^2)$ ，可以用多项式算法加速到 $O(n\log^2 n)$ 。

一般形式

模板。

$n$ 个点可以确定一个多项式 $f(x)$ ，请求出 $f(k)$ 。

核心思想是利用点值的可加性。什么意思？

构造一个函数 $f(x)$ 经过 $P(x_1,y_1),\cdots P(x_n,y_n)$ ，设第 $i$ 个点在 $x$ 轴上的投影为 $P_i^{\prime}(x_i,0)$ 。考虑构造 $n$ 个函数，使得 $f_i(x)$ 的图像经过 $\begin{cases}P_j^{\prime}(x_j,0),(j\neq i)\\P_i(x_i,y_i)\end{cases}$ ，则 $f(x)=\sum_{i=1}^{n}f_i(x)$ 。

设 $f_i(x)=a\prod_{j\ne i}(x-x_j)$ ，然后将 $P_i$ 的坐标代入求出 $a=\cfrac{y_i}{\prod_{j\ne i}(x_i-x_j)}$ ，进而导出：

$f(k)=\sum_{i=1}^n y_i\prod_{i\ne j}\frac{k-x_j}{x_i-x_j}$

这就是拉格朗日插值表达式。

如果要求出系数，用 $O(n^2)$ 卷出 $F(x)=\prod_{i=1}^n (x-x_i)$ ，然后对于每个 $i$ 暴力将 $F$ 除掉 $x-x_i$ 算出 $\cfrac{F(x)}{x-x_i}$ 的各项系数，再乘 $\cfrac{y_i}{\prod_{j\ne i}{x_i-x_j}}$ 即可得到 $f_i$ 的各项系数，加起来即可。时间复杂度为 $O(n^2)$ 。

查看代码

#include <iostream>
#include <cstdio>

using namespace std;
const int MOD = 998244353;

int poww(int a, int b) {
    int res = 1;
    for (; b; b >>= 1, a = 1ll * a * a % MOD)
        if (b & 1) res = 1ll * res * a % MOD;
    return res;
}

int n, k;
int x[2005], y[2005];

int main(void) {
    scanf("%d%d", &n, &k);
    for (int i = 1; i <= n; ++i) scanf("%d%d", x + i, y + i);
    int ans = 0;
    for (int i = 1; i <= n; ++i) {
        int s1 = y[i], s2 = 1;
        for (int j = 1; j <= n; ++j)
            if (i != j) s1 = 1ll * s1 * (k - x[j]) % MOD, s2 = 1ll * s2 * (x[i] - x[j]) % MOD;
        ans = (ans + 1ll * s1 * poww(s2, MOD - 2)) % MOD;
    }
    printf("%d\n", (ans % MOD + MOD) % MOD);
    return 0;
}

连续情况

当已知点的坐标是连续整数时，可以 $O(n)$ 完成拉格朗日插值。以 $0,1\cdots n-1$ 为例：

有：

$f(x)=\sum_{i=0}^{n-1} y_i\prod_{i\ne j}\frac{x-j}{i-j}$

分子是一堆数的乘积挖掉了一个项，典型维护前后缀积。

对于 $i>j$ 的分母，有 $\displaystyle \prod_{j=0}^{i-1}(i-j)=i!$ ，对于 $i<j$ ，有 $\displaystyle \prod_{j=i+1}^{n-1}(i-j)=(-1)(-2)\cdots(i-n+1)$ ，也就是 $(-1)^{n-i+1}(n-i-1)!$ 。

所以：

$f(x)=\sum\limits_{i=0}^{n-1} y_i \frac{p_is_i}{i!(-1)^{n - i + 1}(n-i-1)!}$

集合幂级数

形式幂级数是刻画数列的重要方式，而集合幂级数是处理集合的特殊形式。

由于笔者现在不会 FFT，因此会避开其与 FFT 相关的内容。

定义

形式化地，域 $F$ 上的集合幂级数 $f$ 是 $2^U\rightarrow F$ 的函数，对于每个 $S\subset U$ ，函数 $f$ 均有对应的取值 $f_S\in F$ 。

集合幂级数的加法非常容易，按位相加即可。而乘法比较复杂，有很多种，但是都满足分配律（包括逆运算）。

FMT & FWT

给定序列 ${a_n},{b_n}$ ，求出序列 ${c_n}$ ，满足 $c_i=\sum_{j\oplus k=i}a_j b_k$ ，这样的操作称之为卷积。当其为加号时是我们熟悉的多项式乘法（和卷积），而其为乘号时为狄利克雷卷积。

模板。

用 FMT 解决与卷积和或卷积（或卷积被称为并卷积），用 FWT 来解决异或卷积（异或卷积为对称差卷积）。

我们首先来介绍 FMT（快速莫比乌斯变换）。以或卷积为例，定义序列 $\{a_0,\cdots a_{2^n-1}\}$ 的 FMT 变换为： $\operatorname{FMT}(a)_i=\sum_{j\subseteq i}a_j$ ，IFMT 为上述变换的逆变换。那么或卷积等价于 $c=\operatorname{IFMT}(\operatorname{FMT}(a)\cdot \operatorname{FMT}(b))$ 。

然后我们来介绍 FWT（快速沃尔什变换），其本质就是高维不进位的加法卷积。其定义式为： $\operatorname{FWT}(a)_i = \sum_{j\subseteq i} (-1)^{|i-j|}a_j$ 。IFWT 为上述变换的逆变换。

时间复杂度均为 $O(n\times 2^n)$ 。

性质

FMT 和 FWT 以及它们的逆运算都是一种线性变换。

以 FMT 或卷积为例（FWT 同理），设 $f_i$ 代表第 $i$ 个集合幂级数，那么：

$\prod f_i=\operatorname{IFMT}(\prod\operatorname{FMT}(f_i))$

左边的求积符号代表集合并卷积，右边的代表对应系数相乘。

这样我们就可以解决模板题了：

查看代码

#include <bits/stdc++.h>
using namespace std;
const int N = 1 << 17;
const int P = 998244353;

int n, a[N], b[N], A[N], B[N];
void in(void) {
    memcpy(A, a, N << 2), memcpy(B, b, N << 2);
}
void get(void) {
    for (int i = 0; i < 1 << n; ++i)
        A[i] = 1ll * A[i] * B[i] % P; 
}
void print(void) {
    for (int i = 0; i < 1 << n; ++i) printf("%d ", A[i]);  
    putchar('\n'); 
}
void OR(int *f, int c) {
    for (int i = 1; i < 1 << n; i <<= 1)
        for (int j = 0; j < 1 << n; ++j)
            if (i & j) f[j] = (f[j] + 1ll * c * f[j ^ i] % P + P) % P;
}
void AND(int *f, int c) {
    for (int i = 1; i < 1 << n; i <<= 1)
        for (int j = 0; j < 1 << n; ++j)
            if (i & j) f[j ^ i] = (f[j ^ i] + 1ll * c * f[j] + P) % P;
}
void XOR(int *f, int c) {
    for (int k = 1; k < 1 << n; k <<= 1)
        for (int i = 0; i < 1 << n; i += k << 1)
            for (int j = 0, x, y; j < k; ++j)
                x = f[i | j], y = f[i | j | k], 
                f[i | j] = (x + y) % P, 
                f[i | j | k] = (x - y + P) % P; 
    for (int i = 0; i < 1 << n; ++i) f[i] = 1ll * f[i] * c % P; 
}

int main(void) {
    scanf("%d", &n); 
    for(int i = 0; i < 1 << n; ++i) scanf("%d", a + i);
    for(int i = 0; i < 1 << n; ++i) scanf("%d", b + i);
    
    in(); OR(A, 1); OR(B, 1); get(); OR(A, P - 1); print(); 
    in(); AND(A, 1); AND(B, 1); get(); AND(A, P - 1); print(); 
    int inv = (P + 1) >> 1;
    for (int i = 1; i < n; ++i) inv = 1ll * inv * (P + 1 >> 1) % P; 
    in(); XOR(A, 1); XOR(B, 1); get(); XOR(A, inv); print(); 

    return 0;
}

高维前缀和

集合幂级数 $f$ 经过 FMT 得到的莫比乌斯变换形式与高维前缀和有着密不可分的关系。

FMT 在令第 $i$ 位为 $1$ 开始操作时可以看作在第 $i$ 个维度进行高维前缀和。

子集卷积

定义子集卷积为：

$c_i=\sum_{j\operatorname{or} k = i,j\operatorname{and}k=0}a_j b_k$

注意当 $j\operatorname{or}k = i,j\operatorname{and}k=0$ 的充要条件是 $|j|+|k|=i$ 。设 $f_{i,S}=[|S|=i]f_S$ ， $g$ 同理。令 $H_i=\sum_{j+k=i}f_jg_k$ ，这里的 $H_i$ 代表一个集合幂级数。

等式两边同时 FMT，计算出 $H$ 之后 IFMT 回来即可。时间复杂度 $O(n^2 2^n)$ 。

查看代码

#include <bits/stdc++.h>
using namespace std;
const int P = 1000000009, N = 20; 

int n; 
int f[21][1 << N], g[21][1 << N], h[1 << N], ans[1 << N], popcnt[1 << N]; 

void FMT(int *f, int c) {
    for (int i = 1; i < 1 << n; i <<= 1)
        for (int j = 0; j < 1 << n; ++j)
            if (i & j) f[j] = (f[j] + c * f[j ^ i]) % P; 
    for (int i = 0; i < 1 << n; ++i) f[i] = (f[i] % P + P) % P;
}

int main(void) {
    scanf("%d", &n); 
    for (int i = 0; i < 1 << n; ++i) popcnt[i] = popcnt[i >> 1] + (i & 1); 
    for (int i = 0; i < 1 << n; ++i) scanf("%d", &f[popcnt[i]][i]); 
    for (int i = 0; i < 1 << n; ++i) scanf("%d", &g[popcnt[i]][i]); 
    for (int i = 0; i <= n; ++i) FMT(f[i], 1), FMT(g[i], 1); 
    ans[0] = 1ll * f[0][0] * g[0][0] % P; 
    for (int i = 1; i <= n; ++i) {
        memset(h, 0, sizeof h); 
        for (int j = i; j >= 0; --j)
            for (int k = 0; k < 1 << n; ++k) 
                h[k] = (h[k] + 1ll * f[j][k] * g[i - j][k] % P) % P; 
        FMT(h, -1); 
        for (int j = 0; j < 1 << n; ++j) if (popcnt[j] == i) ans[j] = h[j]; 
    }
    for (int i = 0; i < 1 << n; ++i) printf("%d ", (ans[i] % P + P) % P); 
    return putchar('\n'), 0;
}

快速傅里叶变换（FFT）

万恶之源，快速傅里叶变换，FFT，它可以快速计算多项式的卷积。

模板，代码如下：

快速数论变换（NTT）

Problemset

前 $2$ 个小节没有使用多项式算法。

拉格朗日插值

很有用。

[CF622F] The Sum of the k-th Powers

Portal.

将其写成一个序列： $\sum_{i=1}^1 i^k,\sum_{i=1}^2 i^k,\sum_{i=1}^3 i^k,\cdots$ 。

差分后得到： $1^k,2^k,3^k\cdots$ ，这个数列的通项公式是一个 $k$ 次多项式，那么取前缀和之后是一个 $k+1$ 次多项式。

使用连续情况的拉格朗日插值，不预处理逆元的话是线性对数的。

查看代码

#include <bits/stdc++.h>
using namespace std;
const int P = 1000000007; 

int poww(int a, int b) {
    int res = 1; 
    for (; b; b >>= 1, a = 1ll * a * a % P) if (b & 1) res = 1ll * res * a % P;
    return res; 
}

int n, k, pre[1000005], suf[1000005], fac[1000005]; 

int main(void) {
    scanf("%d%d", &n, &k); fac[0] = pre[0] = suf[k + 3] = 1; 
    for (int i = 1; i <= k + 2; ++i) fac[i] = 1ll * fac[i - 1] * i % P;
    for (int i = 1; i <= k + 2; ++i) pre[i] = 1ll * pre[i - 1] * (n - i) % P; 
    for (int i = k + 2; i >= 1; --i) suf[i] = 1ll * suf[i + 1] * (n - i) % P; 
    int ans = 0, y = 0; 
    for (int i = 1; i <= k + 2;  ++i) {
        y = (y + poww(i, k)) % P;
        int a = 1ll * pre[i - 1] * suf[i + 1] % P; 
        int b = 1ll * fac[i - 1] * (k - i & 1 ? -1 : 1) * fac[k + 2 - i] % P;
        ans = (ans + 1ll * y * a % P * poww(b, P - 2) % P) % P; 
    }
    printf("%d\n", (ans + P) % P); return 0;
}

[集训队互测 2012] calc

Portal.

设 $f_{i,j}$ 表示考虑前 $i$ 个数，值域为 $[1,j]$ ，而且限制取的要上升，这样：

$f_{i,j}=f_{i-1,j-1}\times j + f_{i,j-1}$

最终答案为 $f_{n,k}\times n!$ 。

推测 $f_{n,i}$ 是关于 $i$ 的 $g(n)$ 次多项式。这个转移方程有差分的形式：

$f_{n,i}-f_{n,i-1}=f_{n-1,i-1}\times i$

左边是关于 $i$ 的 $g(n)-1$ 次多项式，右边是关于 $i$ 的 $g(n-1)+1$ 次多项式，因此：

$g(n)-1=g(n-1)+1$

而 $g(0)=0$ ，这样可以得出 $g(n)=2n$ 。

那么我们只需要求出 $f(n,1)\cdots f(n,2n+1)$ ，然后就可以使用拉格朗日插值求出 $f(n,k)$ 。

查看代码

#include <bits/stdc++.h>
using namespace std;

int P; 

int poww(int a, int b) {
    int res = 1; 
    for (; b; b >>= 1, a = 1ll * a * a % P) if (b & 1) res = 1ll * res * a % P;
    return res;
}

int k, n, fac = 1; 
int f[505][1005];

int main(void) {
    scanf("%d%d%d", &k, &n, &P);
    for (int i = 1; i <= n; ++i) fac = 1ll * fac * i % P;
    for (int j = 0; j <= n * 2 + 1; ++j) f[0][j] = 1; 
    for (int i = 1; i <= n; ++i)
		for (int j = 1; j <= n * 2 + 1; ++j)
			f[i][j] = (1ll * f[i - 1][j - 1] * j % P + f[i][j - 1]) % P;
    int ans = 0; 
    for (int i = 1; i <= n * 2 + 1; ++i) {
        int a = 1ll * f[n][i] * fac % P, b = 1;
        for (int j = 1; j <= n * 2 + 1; ++j)
            if (i != j) a = 1ll * a * (k - j) % P, b = 1ll * b * (i - j) % P; 
        ans = (ans + 1ll * a * poww(b, P - 2)) % P; 
    }
    printf("%d\n", (ans + P) % P); return 0;
}

[省选联考 2022] 填树

Portal.

考虑钦定一个最小值 $w$ ，选择的范围就是 $[w,w+K]$ ，然后容斥减去没有取到最小值的 $[w+1,w+K]$ 的部分。

确定一个 $w$ 时答案可以直接进行一次树形 DP 求出。不同取值范围可以划分成不同段，在同一段内的答案是一个关于 $w$ 的多项式（一个点的方案数是一次的，权值是二次的，链求和后依然是多项式），拉格朗日插值维护每一段即可。

查看代码

#include <bits/stdc++.h>
using namespace std;
const int P = 1000000007; 
const int sgn[2] = {1, P - 1}; 

inline int poww(int a, int b) {
    int r = 1; 
    for (; b; b >>= 1, a = 1ll * a * a % P) if (b & 1) r = 1ll * r * a % P; 
    return r; 
}

int n, K, l[205], r[205]; 
vector<int> G[205]; 
int f[205], g[205], d1[205], d2[205]; // f 经过，d 从开始
int ans1, ans2; 

void dfs(int x, int fa, int L, int R) {
    f[x] = g[x] = d1[x] = d2[x] = 0; 
    int ql = max(L, l[x]), qr = min(R, r[x]); 
    int cntx = (ql <= qr) ? qr - ql + 1 : 0; 
    int sumx = (ql <= qr) ? 1ll * (ql + qr) * (qr - ql + 1) / 2 % P : 0; 
    f[x] = d1[x] = cntx; g[x] = d2[x] = sumx; 
    for (int y : G[x]) if (y != fa) {
        dfs(y, x, L, R); f[x] = (f[x] + f[y]) % P; g[x] = (g[x] + g[y]) % P; 
        g[x] = (g[x] + 1ll * d2[x] * d1[y] + 1ll * d1[x] * d2[y]) % P; 
        d2[x] = (d2[x] + 1ll * d2[y] * cntx + 1ll * d1[y] * sumx) % P; 

        f[x] = (f[x] + 1ll * d1[x] * d1[y]) % P; 
        d1[x] = (d1[x] + 1ll * d1[y] * cntx) % P; 
    }
}

int m, a[805]; 
int Y1[805], Y2[805], pre[805], suf[805], inv[405]; 
void solve(int w) {
    m = 0; a[++m] = w + 1; a[++m] = 1000000001; 
    for (int i = 1; i <= n; ++i) {
        if (l[i] > w) a[++m] = l[i]; 
        if (r[i] > w) a[++m] = r[i]; 
        if (l[i] > K) a[++m] = l[i] - K + w; 
        if (r[i] > K) a[++m] = r[i] - K + w; 
    } 
    sort(a + 1, a + m + 1); m = unique(a + 1, a + m + 1) - (a + 1); 
    for (int i = 1; i < m; ++i) {
        int k = min(a[i + 1] - a[i] + 3, n + 5); 
        for (int j = 1; j <= k; ++j) dfs(1, 0, a[i] + (j - 1), a[i] + (j - 1) + K - w), Y1[j] = f[1], Y2[j] = g[1];  
        for (int j = 1; j <= k; ++j) Y1[j] = (Y1[j] + Y1[j - 1]) % P, Y2[j] = (Y2[j] + Y2[j - 1]) % P; 
        int x = a[i + 1] - a[i]; pre[0] = suf[k + 1] = 1; 
        for (int j = 1; j <= k; ++j) pre[j] = 1ll * pre[j - 1] * (x - j + P) % P; 
        for (int j = k; j >= 1; --j) suf[j] = 1ll * suf[j + 1] * (x - j + P) % P; 
        for (int j = 1; j <= k; ++j) {
            int t = 1ll * pre[j - 1] * suf[j + 1] % P * inv[j - 1] % P * inv[k - j] % P * sgn[(k - j) & 1] % P * sgn[w] % P; 
            ans1 = (ans1 + 1ll * t * Y1[j]) % P; ans2 = (ans2 + 1ll * t * Y2[j]) % P;
        }
    }
}

int main(void) {
    scanf("%d%d", &n, &K); 
    for (int i = 1; i <= n; ++i) scanf("%d%d", l + i, r + i); 
    inv[0] = inv[1] = 1; 
    for (int i = 2; i <= n + 5; ++i) inv[i] = 1ll * (P - P / i) * inv[P % i] % P; 
    for (int i = 3; i <= n + 5; ++i) inv[i] = 1ll * inv[i - 1] * inv[i] % P; 
    for (int i = 1; i < n; ++i) {
        int u, v; scanf("%d%d", &u, &v); 
        G[u].emplace_back(v); G[v].emplace_back(u); 
    } solve(0); solve(1); 
    return printf("%d\n%d\n", ans1, ans2), 0;
}

集合幂级数

[ABC212H] Nim Counting

Portal.

给定两个数 $N,K$ ，以及一个长度为 $K$ 的整数数组 $(A_i,A_2...A_K)$ 。

现在通过以下方式生成一个 Nim 游戏：

任意选择一个 $1\le M\le N$ ， $M$ 表示石子堆数。
对于每一堆，其石子数是 $A$ 中任意一个数。

对于 $\sum_{i=1}^N K^i$ 种游戏，求先手获胜的游戏数，答案对 $998244353$ 取模。

$1\le N\le 2\times 10^5$ ，
$1\le A_i,K<2^{16}$ 。
$A_i$ 两两不同。

相当于是给定了一个集合幂级数。 $n$ 堆式子就是 $n$ 个这样的集合幂级数的异或卷积。

利用 FWT 的性质就可以先将集合幂级数们加起来（等比数列求和），再 IFWT 一次。

查看代码

#include <bits/stdc++.h>
using namespace std;
const int N = 16, P = 998244353;

inline int poww(int a, int b) {
    a %= P; int r = 1; 
    for (; b; b >>= 1, a = 1ll * a * a % P) if (b & 1) r = 1ll * r * a % P; 
    return r; 
}

int n, k, f[1 << N];
void FWT(int *f, int c) {
    for (int k = 1; k < 1 << N; k <<= 1)
        for (int i = 0; i < 1 << N; i += k << 1)
            for (int j = 0, x, y; j < k; ++j)
                x = f[i | j], y = f[i | j | k], 
                f[i | j] = (x + y) % P, 
                f[i | j | k] = (x - y + P) % P; 
    for (int i = 0; i < 1 << N; ++i) f[i] = 1ll * f[i] * c % P; 
}

int main(void) {
    scanf("%d%d", &n, &k); 
    for (int i = 1, x; i <= k; ++i) scanf("%d", &x), ++f[x]; 
    FWT(f, 1); 
    for (int i = 0; i < 1 << N; ++i)
        if (f[i] == 1) f[i] = n; 
        else f[i] = 1ll * (poww(f[i], n + 1) - f[i] + P) * poww(f[i] - 1 + P, P - 2) % P; 
    FWT(f, poww(P + 1 >> 1, N)); 
    int ans = 0; 
    for (int i = 1; i < 1 << N; ++i) ans = (ans + f[i]) % P; 
    return printf("%d\n", ans), 0; 
}

[省选联考 2022] 卡牌

Portal.

根号分治， $\le \sqrt{2000}$ 的质数只有 $13$ 个（到 $41$ ）。考虑维护 $f_i$ 代表有多少种选择方式使得恰好覆盖 $i$ 的高维前缀和。预处理时 $f$ 只维护幂次（这样可以减少取模）。

对于含有大质数的卡牌，直接每个都乘上当前质数的选择即可。

最后直接 IFMT 求出真实值。

查看代码

#include <bits/stdc++.h>
using namespace std;
const int MOD = 998244353; 
const int P[13] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41}; 
const int V = 1 << 13; 

inline void add(int &x, int t) { x += t; if (x >= MOD) x -= MOD; }
inline void sub(int &x, int t) { x -= t; if (x < 0) x += MOD; }

int n, m; 
int a[2005], pw[1000005]; 
int c, p[18005], f[V], g[2005][V], h[V]; 

int main(void) {
    scanf("%d", &n); 
    for (int i = 1, x; i <= n; ++i) scanf("%d", &x), ++a[x]; 
    for (int i = pw[0] = 1; i <= n; ++i) pw[i] = (pw[i - 1] << 1) % MOD; 

    for (int i = 1; i <= 2000; ++i) if (a[i]) {
        int tmp = i, msk = 0; 
        for (int j = 0; j < 13; ++j) while (tmp % P[j] == 0) tmp /= P[j], msk |= 1 << j; 
        for (int j = 0; j < V; ++j) if ((j & msk) == msk) f[j] += a[i];
        if (tmp == 43 * 43) tmp = 43; // 43 * 43 < 2000, 43 * 47 > 2000
        if (tmp > 1) {
            for (int j = 0; j < V; ++j)
                if ((j & msk) == msk) g[tmp][j] += a[i];
        }
    }
    
    scanf("%d", &m); 
    while (m--) {
        scanf("%d", &c); int ans = 0, msk = 0; 
        for (int i = 1; i <= c; ++i) scanf("%d", p + i); 
        sort(p + 1, p + c + 1); c = unique(p + 1, p + c + 1) - (p + 1); 

        memcpy(h, f, sizeof h); 
        for (int i = 1; i <= c; ++i) 
            if (p[i] <= 41) {
                for (int j = 0; j < 13; ++j)
                    if (p[i] == P[j]) msk |= 1 << j; 
            } else {
                // 将这个大质数的贡献减去，到时候重新统计
                for (int j = 0; j < V; ++j) 
                    h[j] -= g[p[i]][j];
            }
        
        for (int i = 0; i < V; ++i) h[i] = pw[h[i]]; 
        for (int i = 1; i <= c; ++i) if (p[i] > 41)
            for (int j = 0; j < V; ++j)
                h[j] = 1ll * h[j] * (pw[g[p[i]][j]] - 1) % MOD; // 排掉一个都不选
        
        // IFMT 或卷积逆运算，高维差分
        for (int i = 1; i < V; i <<= 1)
            for (int j = 0; j < V; ++j)
                if (i & j) sub(h[j], h[j ^ i]); 

        for (int i = 0; i < V; ++i) if ((i & msk) == msk) add(ans, h[i]); 
        printf("%d\n", ans); 
    }
    return 0;
}

具体的各类生成函数请参考《组合计数进阶》，这里仅介绍基本内容。 ↩︎