基础算法
多项式乘法
多项式的基础就是多项式乘法,而且能够扩展出来的那些东西都是基于不损失精度的快速数论变换NTT,这是对于下面所有知识点的基础。
令(A(x) = sumlimits_{i = 0}^{n - 1}a_ix^i,B(x)=sumlimits_{i=0}^{n - 1}b_ix^i),当我们计算(A(x)*B(x))的时候,最朴素的方法是(O( ext{n}^2))的,但实际上我们可以做到在(O( ext{n log n}))的时间复杂度内解决这个问题,在这里就用到了快速数论变换NTT。
这个东西讲解比较复杂这里就不讲了,挂一篇很久之前写的学习笔记,静下心把所有式子一步一步推出来应该不用多久就可以学会了。
分治FFT
问题是给定(F(x)=sumlimits_{i=1}^xG(i)F(x-i)),已知(G(x)),求(F(x))。
首先这个东西我们直接做是(O(n^2)),我们考虑计算每个(F(x))对哪些状态有贡献,联想一下CDQ分治的过程,分治的时候每次计算左半部分区间对右半部分的影响,我们可以这样算:
那么我们在分治的时候做一遍多项式乘法就可以把对右半部分所有的贡献都求出来了,复杂度(O( ext{n log}^2 ext{n}))。
多项式求逆
多项式有乘法,那有没有乘法逆元呢?答案是肯定的,我们可以通过多项式求逆来解决这样一个问题,对于已知多项式(A(x))求解另一个多项式(B(x))满足(A(x)*B(x)=1( ext{mod x}^n)),那么多项式(B)就称为多项式(A)的逆元,也记作(A^{-1}(x))。
我们来推导一下多项式逆元的求解过程,首先有:
假设我们已经求出了(A)在模(x^{frac{n}{2}})的意义下的逆元(B'(x)),那么有:
两式相减后可得:
同时除去(A(x)):
两边同时平方:
同时乘上(A(x)):
那么最后可以得到递推式:
推式子的时候可能会遇到一个问题,那就是往下递归的时候可能会让一些奇数项没有计算进去,那么我们刚开始的时候把多项式转换为模(x^{2^p}(nle2^p<2n))意义下的就好了。
那么我们可以用这个式子递归往下求解就好了,由主定理可知复杂度是(O( ext{n log n}))的。
1.4 多项式除法
这里的除法指的是带余除法,也就是让你求(A(x)=B(x) imes C(x)+D(x)),给定(A(x)( ext{mod x}^n),B(x)( ext{mod x}^m)),求(C(x)( ext{mod x}^{n-m+1}),D(x)( ext{mod x}^{m - 1}))。我们令(A^R(x))为多项式(A)系数翻转后的结果,例如(A(x)=1+2x+3x^2,A^R(x)=3+2x+x^2),那么我们把原式的(x)换成(frac 1 x),两边同时乘上(x^n):
我们知道多项数(C)的次数是小于(n-m+1)的,而(D(x) imes x^{n-m+1})的最低次数是不会低于(n-m+1)的,我们把这个式子放到模(x^{n-m+1})意义下,那么(A^R(frac 1 x)=B^R(frac 1 x)C^R(frac 1 x)),求出(C)后回代就可以算出(D)了,复杂度(O( ext{n log n}))。
多项式求导与积分
这个东西应该不用怎么讲,高中数学书已经讲的很清楚了,这里就简单列一下式子。
多项式对数函数
问题是给定多项式(A(x)),求解多项式(B(x))满足(B(x)=ln A(x))。
这个东西推导很简单,直接把两边求导:
最后对(B'(x))进行积分就可以算出(B(x))了,用上文提到的多项式求逆和多项式求导与积分可以实现。
多项式牛顿迭代
问题是给定(G(x)),求解多项式(F(x))使得(G(F(x))=0( ext{mod x}^n))。
首先先引入一个知识点,一个函数(f(x))在某个点(x_0)的泰勒展开有这样一个式子:
在这里(f(x))的(n)阶导数记作(f^n(x))。
我们类似多项式求逆考虑把这个问题划分成子问题递归求解,假设我们已经知道了(G(F_0(x))=0( ext{mod x}^frac{n}{2})),我们在(F_0(x))处对这个函数进行泰勒展开:
我们知道(F(x))与(F_0(x))的前(frac n 2)项时相同的,那么(F(x)-F_0(x))对于(forall i<{frac{n}{2}}),都满足(a_i=0),那么对于展开后第三项以后的所有式子的次数最低项都不会小于(x^n),那么在模意义下它们的值都是(0),那么实际对我们有意义的只有式子的前两项,即:
化一下式子,就得到了:
那么我们就可以用之前学到的多项式求逆递归下去求解这个式子了,复杂度还是(O( ext{n log n} ))。
多项式开方
问题是给定(F(x)),求(G(x))满足(G(x)^2=F(x)( ext{mod x}^n))。
我们把之前的牛顿迭代的式子带进去就好了,设(H(G(x))=G(x)^2-F(x)),那么由牛顿迭代的式子有:
代入后可得:
如果常数项不是(1)或(0)的话还要计算二次剩余,就比较麻烦了,剩下的分治计算即可,复杂度(O( ext{n log n}))。
多项式指数函数
既然可以求对数,那么逆运算应该也是可以的吧,这个问题就是求给定(G(x)),求(F(x))满足:
先同时对两边取自然对数,移项可得:
我们令(H(F(x))=ln F(x)-G(x)),套用之前提到的牛顿迭代,可以得到:
对于(H(F(x))=ln F(x)-G(x)),求导可得(H'(F(x))=frac{1}{F(x)}),那么把这两个带进去:
整理得:
我们利用之前的多项式对数函数即可,复杂度(O( ext{n log n})),常数巨大。
常见应用
第一类斯特林数
第一类斯特林数(S(n,m))是表示把(n)个物品排成(m)个环的方案数,首先有递推式:
这个式子我们可以写出它的生成函数,也就是:
理解起来应该很容易,乘的这个东西分别对应两种转移,那么这个式子我们就可以用分治+NTT做到(O( ext{n log}^2 ext{n})) 的复杂度,
第二类斯特林数
第二类斯特林数(S(n, m))是表示把(n)个物品分成(m)个集合的方案数,首先有递推式:
但是这个式子我们不方便用多项式的运算来优化,我们把这个式子写成容斥的形式:
意义是我们每次枚举空集合的个数,计算至少有(k)个空集合的方案数,剩下的随便放,又因为集合是无序的,所以最后还要除去(m!),我们再把组合数给展开:
用NTT优化多项式乘法解决这个问题就可以做到(O( ext{n log n}))了。
例题讲解
[BZOJ3456] 城市规划
求(n)个点的带标号无向连通图个数,(nle10^5, ext{mod = 998244353})。
我们设(f(i))为(i)个点的带标号无向连通图个数,(g(i))表示(i)个点的带标号图的个数,那么(g(i)=2^{inom 2 i}),我们可以写出(f)的递推式:
把式子化开,我们可以得到:
把((i-1)!)除去:
我们发现把左边那个式子刚好是右边那个求和的第(i)项,合并这两个式子:
设:
那么(H(x)=G(x)*F(x)),所以(F(x)=H(x)*G^{-1}(x)),多项式求逆即可。
[HDU5322] Hope
对于1到n这n个数的任何一个排列A可以这样计算其价值:对所有下标(i)找到最小的(j)满足(j>i)且(A[j]>A[i]),然后(i)和(j)之间连边,最后所有连通块个数之积的平方就是该排列的价值,问所有排列的价值和是多少,(nle 10^5, ext{mod = 998244353})。
首先发现一个小性质,对于数列中最大的数一定和它前面的所有数在同一个集合,那么我们设(dp_i)为(i)个数组成的排列的价值,每次枚举最大的数放在第几位:
把这个式子拆一下:
就可以用分治FFT求解了,每次只要把前面的dp值除一个阶乘就好了。
[LOJ2541] 「PKUWC2018」猎人杀
有(n)个人,每个人有一个权值(w_i),每次随机杀一个人,杀第(i)个人的概率是(displaystylefrac{w_i}{sum_j[j is alive]}),求第一个人最后一个死的概率,对(998244353)取模。
第一个人最后一个死就代表恰好有(0)个人在第一个人之后死,算这个是个很套路的东西用容斥转化为至少,那么设(f(S))为至少有(S)集合的人在第一个人之后死的概率,那么(ans=sum(-1)^{|S|}f(S)),现在我们转化成了求(f(S)),实际上这个东西我们可以写出一个式子:
这个东西理解起来也不是很难,每次只有选择(S)集合内的人击杀或者(1)击杀的时候才会影响他们的相对死亡顺序,而选中这里面的人击杀时必须要击杀第一个人才满足(S)集合都在(1)之后死,那么我们发现(f(S))只与(sum_w[win S])有关,又因为题目中的条件有(sum wle10^5),我们可以考虑计算出每个值的容斥系数,可以直接设(f[i][j])表示考虑前(i)个人,权值和为(j)时的容斥系数,用背包的方法转移就是(f[i][j] = f[i - 1][j] - f[i - 1][j - w_i]),但这样子显然是无法通过的,我们需要优化这个方法。
我们设出这个东西的生成函数,第(i)项(a_ix^i)代表当前(dp_i=a_i),那么每次转移就是乘上((1-x^{w_i})),分别代表是否把当前这个人加入集合的决策,那么这个生成函数就是:
这个东西直接分治,合并两个区间的信息的时候用NTT计算就好了,这个东西开数组有点麻烦,那么我们可以类似于线段树动态开点回收空间的方法一样,用完以后利用以前的空间,分治出最深的一条链深度应该是(log n)的,那么我们开(log n)个数组就好了,复杂度(O(nlog ^2n))。
[AGC005F] Many Easy Problems
给你一棵(n)个点的树和一个整数(k),设(S)是树上某些点的集合,(f(S))是最小的包含(S)的联通子图的大小,(n)个点选(k)个点一共有(inom k n)种方案,你需要对于(kin[1,n])求出所有(k)个点的集合的答案,对924844033取模。
我们直接算(k)的所有答案是不好算的,这个时候考虑算每个点对于每个(k)的贡献,我们发现一个点产生贡献当且仅当选中的点不全在以它相邻的点为根的子树内,所以一个点(u)对一个(k)的贡献为:
观察这个式子发现,一个点的子树大小会在它的父亲统计它时记为(size_{u})统计一遍,在它统计它的父亲时作为(n-size_{u})统计一遍,那么我们设(cnt_x)为子树大小为(x)被统计到答案内的次数,那么:
前一部分很好算,后一部分把组合数拆开,就是:
这个东西就可以做卷积了,顺带提一句(924844033 = 2^{21} imes3^2 imes7^2+1),并且它的原根是(5)而不是(3)。
[LOJ2058] 「TJOI / HEOI2016」求和
设(S(i,j))表示第二类斯特林数,求(f(n)=sum_{i=0}^nsum_{j=0}^iS(i,j) imes2^j imes j!)。
我们首先发现(j>i)时,(S(i,j)=0),那么我们可以把(j)的枚举范围扩展到(n),那么就变成了求:
用第二类斯特林数的容斥公式:
那么要求的变成了:
(sum_{i = 0}^n(j-k)^i)用等比数列求和公式可以可以变成(frac{1-(j-k)^{n+1}}{1-(j-k)}),那么我们可以设:
那么:
把(f)和(g)卷起来之后再求和就可以了。
[CF438E] The Child and Binary Tree
求有多少种不同的二叉树,满足二叉树每个点的权值(in C),求出所有权值和(in[1,m])的二叉树种类数,对(998244353)取模。
设(f_i)表示权值为(i)的二叉树种类数,那么有如下递推式:
我们设出(f)答案与存在的物品的权值(c)的生成函数,那么:
最后那个加(1)是因为(f_0=1),但是在答案的卷积中不会计算(f_0)。
我们利用一元二次方程求根公式可得:
我们发现取正号的时候这个东西不收敛,所以舍去。
那么(f = frac{1-sqrt{1 - 4c}}{2c}),多项式开方+多项式求逆即可。
[LUOGU4389] 付公主的背包
你有(n)个物品,每个物品有无穷多个,求选择一些物品恰好体积为(v)的方案数,求出(vin[1,m])的所有答案,(n,mle10^5, mod = 998244353)。
我们还是按生成函数的套路,设(a_ix^i)这一项的系数表示体积为(i)的方案数,那么每加入一件物品就相当于给这个生成函数乘上了一个((1+x^{v_i}+x^{2v_i}+...)),也就是(frac{1}{1-x^{v_i}}),那么我们最后要求的就是:
我们发现直接算不好算,先对整个函数取(ln),那么我们来考虑快速计算每一项的(ln),也就是如何快速计算:
我们先对这个东西求导,也就是:
把(frac{1}{1-x^v})写成无穷级数求和的形式,那么:
把这个积分回去,有:
我们就可以用(O( ext{n ln n})),枚举倍数的方法计算出(ln(prod_{i = 1}^nfrac{1}{1-x^{v_i}}))的每一项系数了,我们最后再用多项式exp把这个转化回来就可以了,复杂度是(O( ext{n log n}))。
[LUOGU4705] 玩游戏
对于(kin[1,t]),求(frac{sum_{i = 1}^nsum_{j = 1}^m(a_i+b_j)^k}{nm})。
先不管分母,把分子用二项式定理展开:
这样子就写成了卷积的形式,那么我们只要能快速求出(sum_{k=1}^na^k),就可以用NTT计算卷积了。
我们写出这个东西的生成函数:
(x^k)的系数表示(k)次幂和,那么我们用等比数列求和公式解出来可以得到:
令(f(x))为这些生成函数的和,那么:
这个东西不好算,我们发现:
我们从对数函数角度考虑,又可以发现:
设(g(x)=displaystylesum_{i = 1}^ndisplaystylefrac{-a_i}{1-a_ix}),那么(f(x)=-x imes g(x)+n)。
(g(x))也很好算,化一下式子就变成了:
这样子(g)就可以用分治+NTT算,算出(g)后再推出(f),(f)中(x_i)的系数就是(displaystylesum_{j = 1}^na_j^i),那么我们代回原式再做一遍卷积就好了,总的复杂度为(O(nlog^2n))。
总结
多项式的题目其实有很多套路,而且多项式难的地方应该在推式子而不是模板上面,毕竟多项式的模板并没有什么可以改的地方,在OI中一般多项式不会直接拿来出题,可能有对(O(n^2))DP的优化,也有那种直接计算生成函数系数的,多项式可能比较烦的是打模板,因为这种东西根本不能调的,在后面我会附上包含我这上面所有操作的模板的,然后为了准备这篇学习笔记也花了蛮多时间整理博客,希望看完这篇学习笔记后你能对多项式相关内容有更深的理解:)。
多项式模板
#include <bits/stdc++.h>
using namespace std;
const int N = 1 << 18 | 1;
const int mod = 998244353;
inline int mul(int x, int y) { return 1ll * x * y % mod; }
inline int add(int x, int y) { return (x += y) < mod ? x : x - mod; }
inline int qpow(int _, int __) {
int ___ = 1;
for (; __; _ = 1ll * _ * _ % mod, __ >>= 1)
if (__ & 1) ___ = 1ll * ___ * _ % mod;
return ___;
}
namespace Poly {
static int rev[N];
int Get_Rev(int x) {
int limit = 1, k = 0;
while (limit < x) limit <<= 1, ++ k;
for (int i = 0; i < limit; ++ i)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
return limit;
}
void NTT(int *a, int n, int fh) {
for (int i = 0; i < n; ++ i)
if (i < rev[i]) swap(a[i], a[rev[i]]);
for (int Wn, limit = 2; limit <= n; limit <<= 1) {
Wn = qpow(fh ^ 1 ? qpow(3, mod - 2) : 3, (mod - 1) / limit);
for (int W = 1, j = 0; j < n; j += limit, W = 1)
for (int i = j; i < j + (limit >> 1); ++ i, W = mul(W, Wn)) {
int a1 = a[i], a2 = mul(W, a[i + (limit >> 1)]);
a[i] = add(a1, a2), a[i + (limit >> 1)] = add(a1, mod - a2);
}
}
if (fh ^ 1) for (int inv = qpow(n, mod - 2), i = 0; i < n; ++ i)
a[i] = mul(a[i], inv);
}
void Add(int *a, int *b, int *c, int len, int typ) {
for (int i = 0; i < len; ++ i)
c[i] = add(a[i], typ ? b[i] : mod - b[i]);
}
void Mul(int *a, int n, int *b, int m, int *c) {
int limit = Get_Rev(n + m), a_[N], b_[N];
for (int i = 0; i < limit; ++ i) {
a_[i] = i < n ? a[i] : 0;
b_[i] = i < m ? b[i] : 0;
}
NTT(a_, limit, 1), NTT(b_, limit, 1);
for (int i = 0; i < limit; ++ i)
c[i] = mul(a_[i], b_[i]);
NTT(c, limit, -1);
}
void Get_Inv(int *a, int *b, int len) {
if (len & 1) return (void)(b[0] = qpow(a[0], mod - 2));
Get_Inv(a, b, len >> 1);
static int A[N], B[N];
Mul(a, len, b, len, A);
Mul(A, len, b, len, B);
for (int i = 0; i < len; ++ i)
b[i] = add(b[i], add(b[i], mod - B[i]));
}
void Inv(int *a, int *b, int len) { Get_Inv(a, b, Get_Rev(len)); }
void Der(int *a, int *b, int len) {
for (int i = 0; i < len - 1; ++ i)
b[i] = mul(i + 1, a[i + 1]);
b[len - 1] = 0;
}
void Int(int *a, int *b, int len) {
for (int i = len - 1; i; -- i)
b[i] = mul(a[i - 1], qpow(i, mod - 2));
b[0] = 0;
}
void Ln(int *a, int *b, int len) {
static int A[N], B[N];
Inv(a, A, len), Der(a, B, len);
Mul(A, len, B, len, b), Int(b, b, len);
}
void Get_Exp(int *a, int *b, int len) {
if (len & 1) return (void)(b[0] = 1);
Get_Exp(a, b, len >> 1);
static int A[N], B[N];
Ln(b, A, len), Add(a, A, B, len, 0);
++ B[0], Mul(b, len, B, len, b);
}
void Exp(int *a, int *b, int len) { Get_Exp(a, b, Get_Rev(len)); }
void Get_Sqrt(int *a, int *b, int len) {
if (len & 1) return (void)(b[0] = calcsqrt(a[0]));
Get_Sqrt(a, b, len >> 1);
static int A[N], B[N];
Add(b, b, A, len, 1), Inv(A, B, len), Mul(b, len, b, len, A);
Add(A, a, A, len, 1), Mul(A, len, B, len, b);
}
void Sqrt(int *a, int *b, int len) { Get_Sqrt(a, b, Get_Rev(len)); }
void Div(int *a, int n, int *b, int m, int *c, int *d) {
static int A[N], B[N], C[N];
for (int i = 0; i < m; ++ i) A[m - i - 1] = b[i];
Inv(A, B, n - m + 1);
for (int i = 0; i < n; ++ i) A[n - i - 1] = a[i];
Mul(A, n, B, n - m + 1, C);
for (int i = 0; i < n - m + 1; ++ i) c[i] = C[n - m - i];
Mul(b, m, c, n - m + 1, A);
for (int i = 0; i < m; ++ i) d[i] = add(a[i], mod - A[i]);
}
}
int main() {
// 首先定义n - 1次多项式是长度为n的多项式。
// Mul(a, n, b, m, c) 表示把长度分别为n, m的两个多项式a, b乘起来的结果,答案存在c数组中。
// Inv(a, b, len) 表示把长度为len的多项式a在模x^len下的结果计算在c数组中。
// Int(a, b, len) 表示把长度为len的多项式a积分后的结果计算在数组b中。
// Der(a, b, len) 表示把长度为len的多项式a求导后的结果计算在数组b中。
// Ln(a, b, len) 表示把长度为len的多项式a取ln后的结果计算在数组b中。
// Sqrt(a, b, len) 表示把长度为len的多项式a开方后的结果计算在数组b中,但是在此处我没有写计算二次剩余的函数,需要自己实现。
// Exp(a, b, len) 表示把长度为len的多项式a求exp后的结果计算在数组b中。
// Div(a, n, b, m, c, d) 表示把长度为n,m的两个多项式a,b做带余除法的结果,商存在c数组,余数存在d数组中。
return 0;
}