多项式全家桶学习笔记(进阶)
好的继续一点一点啃吧
Todo List
- “鸽了”的实现
- 多项式复合
- 多项式复合逆
- 常系数非齐次线性递推
一.数学前置
1.1 二项式反演
等价于
Proof
证明:参考
显然只要证明前推后即可。
考虑二式的右边。
得证。
1.2 拉格朗日反演
则
Proof
证明:
设
所以
而当\(i\neq n\)时,显然对左式没有贡献,可得
注意到
所以
即
证毕。
1.3 单位根反演
设 \(\omega\) 表示 \(n\) 次本原单位根。
我们有
Proof
证明:
当 \(n|k\) 时,有 \(\omega^{ik}=1\),显然成立!
当 \(n\nmid k\) 时,有
得证!
1.4 斯特林数
先给出斯特林数的定义:
第一类斯特林数: \(\begin{bmatrix}n\\k\end{bmatrix}\) 表示将 \(n\) 个两两不同的元素,划分为 \(k\) 个非空圆排列的方案数。
第二类斯特林数: \(\begin{Bmatrix}n\\k\end{Bmatrix}\) 表示将 \(n\) 个两两不同的元素,划分为 \(k\) 个非空子集的方案数。
有如下递推式:
有用的式子:
斯特林反演:
等价于
通项公式:
这个可以用二项式反演或者直接组合意义容斥来证
二.进阶多项式算法
2.1 Bluestein's Algorithm
orz myy!
注:原论文似乎有点问题,如果发现这里也有问题欢迎指出。
\(\text{Bluestein's Algorithm}\) 支持任意长度 \(\text{CZT}\),即给出多项式 \(B(x)\),计算
考虑如何将其化为普通卷积的形式。
下面我们用一种膜法,注意到这么一个恒等式:
正确性显然。
通常称为 \(\text{Bluestein}\) 的第一个拆乘积的方法。
代入上面的式子:
发现是个卷积,于是直接 NTT 就好了……
吗?
我们发现指数不是整数,这意味着只有在 \(c\) 有二次剩余的时候才可以,但是通常不存在。
有没有不依赖二次剩余的方法呢?
我们发现了另一个式子:
通常称为 \(\text{Bluestein}\) 的第二个拆乘积的方法。
于是
这就是一个卷积的形式了,直接 NTT 即可。
复杂度 \(O(n\log n)\)。
Code
vec CZT(vec f,int len,int c,int m){
vec A(len),B(len+m),g(m);
for(int i=0;i<len;i++)A[i]=1LL*qpow(c,P-1-F(i))*f[i]%P;
for(int i=0;i<len+m;i++)B[i]=qpow(c,F(i));
reverse(A.begin(),A.end());
A=mul(A,B,len*2+m);
for(int i=0;i<m;i++)g[i]=1LL*qpow(c,P-1-F(i))*A[i+len-1]%P;
return g;
}
2.2 多项式多点求值
多项式多点求值,即求多项式 \(A(x)\) 在点 \(x_1,x_2,\cdots,x_n\) 上分别取到的值。
Sol 1
常用做法。
构造函数
设 \(A(x)=P_0(x)D(x)+R(x)\),且 \(\deg R<\deg P_0=\lfloor\frac{n}{2}\rfloor\)。
那么对于 \(\forall x\in \{x_1,x_2,\cdots,x_{\lfloor\frac{n}{2}\rfloor}\}\),都有 \(A(x)=R(x)\),于是我们就将一个 \(\deg =n\) ,有 \(n\) 个点的问题,转化为一个 \(\deg <\lfloor\frac{n}{2}\rfloor\) ,有 \(\lfloor\frac{n}{2}\rfloor\) 个点的问题。
另一半点值与上面同理。
因此分治 NTT 即可,注意计算 \(P_0(x)\) 和 \(P_1(x)\) 同样可以用一次分治 NTT 预处理。
设计算的复杂度为 \(T(x)\),则 \(T(x)=2T(\frac{x}{2})+O(n\log n)\),此处 \(O(n\log n)\) 的部分是算 \(D(x),R(x)\) 的复杂度,容易发现预计算复杂度也是 \(T(x)\)。
由主定理可得总复杂度为 \(O(n\log^2 n)\)。
实际应用时可以通过用暴力秦九韶进行剪枝来优化常数。
Code
inline vec mod(vec f,vec g,int n,int m){
f.resize(n),g.resize(m);
vec ff(n),gg(n),iv(n),t(n);
for(int i=0,t=n-1;i<n;i++,t--)ff[i]=f[t];
for(int i=0,t=m-1;i<m;i++,t--)gg[i]=g[t];
int len=n-m+1;
for(int i=len;i<n;i++)ff[i]=gg[i]=0;
getinv(gg,iv,len);
ff=mul(ff,iv,len<<1);
vec q(n);
for(int i=0,t=len-1;i<len;i++)q[i]=ff[t--];
for(int i=len;i<n;i++)q[i]=0;
t=mul(q,g,n<<1);
vec r(m-1);
for(int i=0;i<m-1;i++)r[i]=(f[i]-t[i]+P)%P;
return r;
}
void build(int k,int l,int r,const vec&a){
if(l==r){
Len[k]=1;
p[k].resize(2);
p[k][0]=P-a[l];
p[k][1]=1;
return;
}
int mid=l+r>>1;
build(k<<1,l,mid,a);
build(k<<1|1,mid+1,r,a);
Len[k]=Len[k<<1]+Len[k<<1|1];
p[k]=mul(p[k<<1],p[k<<1|1],Len[k]+1);
p[k].resize(Len[k]+1);
}
void solve(int k,int l,int r,const vec&a,const vec&A,vec&ans){
if(Len[k]<=500){
int m=Len[k]-1;
for(int i=l;i<=r;i++)
for(int j=m;j>=0;j--)
ans[i]=(1LL*ans[i]*a[i]+A[j])%P;
return;
}
if(l==r){ans[l]=A[0];return;}
int mid=l+r>>1;
vec R;
R=mod(A,p[k<<1],Len[k],Len[k<<1]+1);
solve(k<<1,l,mid,a,R,ans);
R=mod(A,p[k<<1|1],Len[k],Len[k<<1|1]+1);
solve(k<<1|1,mid+1,r,a,R,ans);
}
vec evaluation(vec f,vec a,int n,int m){
build(1,1,m,a);
if(n>m){
f=mod(f,p[1],n,Len[1]+1);
}
vec ans(m+1);
solve(1,1,m,a,f,ans);
return ans;
}
Sol 2
转置做法。
先简单解释一下 转置原理 的核心思想:
对矩阵 \(\mathbf M\) 和向量 \(\mathbf v\),为了快速计算 \(\mathbf {Mv}\),我们可以先考察怎么计算 \(\mathbf M^\mathsf T\mathbf v\),也就是矩阵的转置。设 \(\mathbf M\) 可以分解为 \(\mathbf M=\mathbf E_1 \mathbf E_2\cdots \mathbf E_k\),那么 \(\mathbf M^\mathsf T=\mathbf E_1^\mathsf T \mathbf E_2^\mathsf T\cdots \mathbf E_k^\mathsf T\)。在这些初等矩阵中,它们的功能可以分为两类:
- \(\mathbf {Ev}\) 的效果是让向量 \(\mathbf v\) 的第 \(i\) 位乘上 \(c\),那么其转置效果不变。
- \(\mathbf{Ev}\) 的效果是让向量 \(\mathbf v\) 的第 \(i\) 位乘上 \(c\) 然后加到第 \(j\) 位,那么它的转置的效果就是让向量的第 \(j\) 位乘上 \(c\) 然后加到第 \(i\) 位。
接下来我们考虑多点求值的过程,其本质是 VanderMonde 矩阵
对应的线性变换 \(\mathbf{V}(x_0,x_1,\dots,x_{n-1})\mathbf v\),其中 \(\mathbf v\) 是系数向量 \(\begin{bmatrix}v_0\\v_1\\\vdots\\v_{n-1}\end{bmatrix}\)。
根据转置原理,我们考虑其转置 \(\mathbf{V}^\mathsf T(x_0,x_1,\dots,x_{n-1})\mathbf v\):
考虑这个转置所求的是什么。
第 \(k\) 项是
而为了求这个,我们可以这么做:
- 分治,维护左右两部分的 \((u_1,L)\) 和 \((u_2,R)\),也就是答案的分子和分母的部分,然后合并为 \((u_1R+u_2L,LR=P)\)。
- 答案就是 \(\frac{u}{P}=uP^{-1}\),这里 \(P\) 指的是分母,即 \(\prod_i(1-x_ix)\)。
我们只要算出其转置,就可以得到新的多点求值做法。在此之前,我们还要考虑一下多项式乘法 \(\mathbf{MUL}(A)\mathbf v\) 的转置 \(\mathbf {MUL}(A)^\mathsf T\mathbf v\)。
对多项式 \(A(x)=\sum_{i}a_ix^i\),在多项式乘法中会把每个 \(v_j\) 乘上 \(a_i\) 然后贡献给 \(i+j\),因此根据转置原理,转置的多项式乘法就是把每个 \(v_{i+j}\) 乘上 \(a_i\) 然后贡献给 \(i\)。可以看到,这就是一个标准的减法卷积。此外,如果原来的变换没有溢出部分,那么转置后的溢出部分我们也要扔掉。
于是,我们就可以得到新的多点求值做法:
- 整体分治(其实就是线段树的 build),算出每个线段树节点对应的区间里 \((1-x_ix)\) 的乘积。
- 对要求的 \(m-1\) 次多项式,我们对它计算 \(\mathbf{MUL}(P^{-1}\bmod x^m)^\mathsf T\) 并保留前 \(n\) 位。
- 从线段树自顶向下递归,令下传的两个向量分别为 \(\mathbf l=\mathbf {MUL}(R)^\mathsf T\mathbf v\) 和 \(\mathbf r=\mathbf{MUL}(L)^\mathsf T\mathbf v\)。
- 最后得到的叶节点的向量长度均为 \(1\),对应于该处的点值。
新做法不需要多项式取模,而且常数也有显著提升。
Code
vec mulT(vec f,vec g){
int n=f.size(),m=g.size()-1;
init(n);
f.resize(lmt),g.resize(lmt);
reverse(g.begin()+1,g.end());
vec ret(lmt);
NTT(f,lmt,1);NTT(g,lmt,1);
for(int i=0;i<lmt;i++)ret[i]=1LL*f[i]*g[i]%P;
NTT(ret,lmt,-1);
ret.resize(n-m+1);
return ret;
}
pair<vec,vec> mul2(vec a,vec b1,vec b2){
int n=a.size(),m1=b1.size()-1,m2=b2.size()-1;
init(n);
a.resize(lmt),b1.resize(lmt),b2.resize(lmt);
reverse(b1.begin()+1,b1.end());
reverse(b2.begin()+1,b2.end());
vec ret1(lmt),ret2(lmt);
NTT(a,lmt,1);NTT(b1,lmt,1);NTT(b2,lmt,1);
for(int i=0;i<lmt;i++)ret1[i]=1LL*a[i]*b1[i]%P,ret2[i]=1ll*a[i]*b2[i]%P;
NTT(ret1,lmt,-1);NTT(ret2,lmt,-1);
ret1.resize(n-m1+1);
ret2.resize(n-m2+1);
return make_pair(ret1,ret2);
}
void build(int n){
lc.assign(n*2,0);
rc.assign(n*2,0);
deg.assign(n*2,0);
p.assign(n*2,{0,0});
q.assign(n*2,{0,0});
pos.resize(n);
int cnt=0;
function<int(int,int)> dfs=[&](int l,int r){
if(l==r){
pos[l]=cnt;
deg[cnt]=1;
return cnt++;
}
int ret=cnt++,mid=l+r>>1;
lc[ret]=dfs(l,mid);
rc[ret]=dfs(mid+1,r);
deg[ret]=deg[lc[ret]]+deg[rc[ret]];
return ret;
};
dfs(0,n-1);
}
vec eval(vec&f,const vec&x){
int n=f.size();
build(n);
rep(i,0,n-1)q[pos[i]]={1,P-x[i]},deg[pos[i]]=1;
Rep(i,n*2-2,0)if(lc[i])q[i]=mul(q[lc[i]],q[rc[i]],deg[i]+1);
f.resize(n<<1);
vec g;
int len=q[0].size()+1;
getinv(q[0],g,len);
g.resize(len);
p[0]=mulT(f,g);
rep(i,0,n*2-1)if(lc[i])tie(p[lc[i]],p[rc[i]])=mul2(p[i],q[rc[i]],q[lc[i]]);
vec ans(n);
rep(i,0,n-1)ans[i]=p[pos[i]][0];
return ans;
}
vec evaluation(vec f,vec x){
int n=f.size(),m=x.size();
f.resize(max(n,m));
x.resize(max(n,m));
vec ret=eval(f,x);
ret.resize(m);
return ret;
}
2.3 多项式快速插值
多项式快速插值,即给定 \(n\) 个点 \((x_i,y_i)\),你需要求出一个多项式 \(f(x)\),使得 \(\deg f<n\) 且 \(\forall i,f(x_i)=y_i\)。
首先根据拉格朗日插值公式可得
这时如果你直接模拟那么复杂度是 \(O(n^2)\) 的,我们考虑如何优化。
先考虑系数部分的
不难发现如果设
那么
(最后一步是导数的定义)
那么使用多点求值就可以做到 \(O(n\log^2n)\)。
下设 \(\frac{y_i}{P'(x_i)}=v_i\)。
考虑和多点求值类似的做法,我们还是构造函数
设
那么容易发现
分治 NTT 即可。
设这部分复杂度为 \(T(x)\),则 \(T(x)=2T(\frac{x}{2})+O(n\log n)\),此处 \(O(n\log n)\) 的部分是计算 \(f(x)\) 的复杂度。由于 \(P(x)\) 在多点求值时已经算过,所以不需再算。
由主定理可得总复杂度为 \(O(n\log^2 n)\)。
讲一些卡常的技巧:注意到在递推的时候大部分都是先 \(n\) 大小 IDFT 再 \(2n\) 大小 DFT,可以考虑直接维护 DFT 结果,可以显著优化常数。
Code
鸽了
2.4 普通多项式转下降幂多项式
给定普通多项式
求下降幂多项式
使得 \(F(x)=G(x)\)。
首先有
所以
所以
多点求值预处理 \(F(t)\) ,然后卷积即可。
复杂度 \(O(n\log ^2n)\)。
Code
鸽了
2.5 下降幂多项式乘法
给定 \(n\) 次下降幂多项式 \(A(x)\) 和 \(m\) 次下降幂多项式 \(B(x)\),求一个 \(n+m\) 次的下降幂多项式 \(F(x)\),使 \(F(x)=A(x)B(x)\)。
对于下降幂多项式 \(F(x)\),考虑它的点值的 EGF
当 \(F(x)=x^\underline c\) 时
设
则
其中
也就是说,\(F^\#\) 就是把 \(F\) 当成普通多项式,乘上 \(e^x\) 的结果。
同理可得其逆运算就是 \(e^{-x}\)。
点值 EGF 点乘即可,\(e^x\) 和 \(e^{-x}\) 可以手算。
复杂度 \(O(n\log n)\)。
Code
void FDT(vec&A,int len,int tp){
init(len<<1);
vec ee(lmt);
if(A.size()<lmt)A.resize(lmt);
if(tp==-1)for(int i=0;i<lmt;i++)A[i]=1LL*A[i]*invfac[i]%P;
for(int i=0;i<len;i++){
if(tp==-1&&i&1)ee[i]=P-invfac[i];
else ee[i]=invfac[i];
}
for(int i=len;i<lmt;i++)ee[i]=A[i]=0;
NTT(A,lmt,1);NTT(ee,lmt,1);
for(int i=0;i<lmt;i++)A[i]=1LL*A[i]*ee[i]%P;
NTT(A,lmt,-1);
if(tp==1)for(int i=0;i<lmt;i++)A[i]=1LL*A[i]*fac[i]%P;
for(int i=0;i<lmt;i++)ee[i]=0;
}
vec mulDown(vec f,vec g,int len){
FDT(f,len,1);FDT(g,len,1);
for(int i=0;i<len;i++)f[i]=1LL*f[i]*g[i]%P;
FDT(f,len,-1);
return f;
}
2.6 下降幂多项式转普通多项式
给定下降幂多项式
求普通多项式
使得 \(F(x)=G(x)\)。
根据上题结论,下降幂多项式的点值是与 \(e^x\) 的卷积。
那我们只要求出 \(F(x)\) 在 \(0,1,\dots,n-1\) 的点值,然后用多点求值就可以得到原多项式。
复杂度 \(O(n\log^2 n)\)。
2.7 多项式复合
2.8 多项式复合逆
2.9 快速沃尔什变换
给定序列 \(A\) 和 \(B\),求
其中 \(\oplus\) 是位运算
尝试采取和 DFT 相类似的手法,考虑另一个数列变换 \(\operatorname{FWT}(A)\),称之为沃尔什变换。
为了能做到快速进行乘法,我们需要使得
其中 \(\cdot\) 表示点乘。
下面分与、或、异或三种进行讨论。
2.9.1 或
显然可得
构造
所以
下面考虑如何快速计算和逆运算。
和 FFT 一样,将 \(a\) 拆成下标最高位为 \(0\) 与下标最高位为 \(1\) 两个序列,分别记为 \(a_0\) 和 \(a_1\)。
所以可得
这里 \(\operatorname{merge}\) 表示拼接,\(+\) 表示对应位置相加。
分治就好了。
逆运算:根据上面那个式子我们同样可以得到 \(a=\operatorname{merge}(a_0,a_1-a_0)\),就可以了。
2.9.2 与
完全同理地,我们定义
推理过程略。
2.9.3 异或
这里我们定义运算 \(\otimes\),使得 \(x\otimes y=\operatorname{popcount}(x\& y)\bmod 2\)。
因此我们有 \((i\otimes j)\oplus(i\otimes k)=i\otimes(j\oplus k)\)。
那么定义
所以
所以
以上复杂度均为 \(O(n\log n)\)。
2.10 常系数齐次线性递推
常系数齐次线性递推,是指给定 \(f[1,2,\dots,k]\) 和 \(a[0,1,\dots,k-1]\),且对 \(\forall m\ge k\),都有
求 \(a_n\) 的值。
先考虑一般是如何处理这个问题的。
考虑矩阵
那么
所以
使用矩阵快速幂就可以做到 \(O(k^3\log n)\) 的复杂度。
考虑如何优化这个过程,由于我们只要求 \(a_n\),所以我们只关心右式的第一项即可。
下记 \(m=n-k+1\),且
现在假设我们不知道从哪找来一个数组 \(c\),使得
则
下考虑如何构造 \(c\)。
考虑构造一个
其中 \(Q,G,R\) 是一些多项式,且使得 \(\deg R<\deg G\)。
令 \(\deg G=k\) 且 \(G(A)=0\),有:
因此我们如果构造出了 \(G\),就直接求出 \(x^n\bmod G\) 就可以了,使用快速幂即可,复杂度为 \(O(k\log k\log n)\) 的。
因此只要取
就行。
结论:如果递推系数为 \(\{a\}\),那么 \(g_{k-i}=-a_i,g_k=1\)。
Proof
先定义特征值:如果 $$ (\lambda I-A)v=0 $$ 那么 $\lambda$ 是 $A$ 的特征值,$v$ 是 $A$ 的特征向量。再定义特征多项式:
是特征多项式。
根据 Cayley-Hamilton 定理可得,\(G(A)=0\),因此 \(G\) 就是我们要求的。
得证。
复杂度 \(O(k\log k\log n)\)。