这篇文章会从零开始证一遍多项式全家桶中的基础部分,初中水平就能理解(因为我能理解qaq)。
讲了泰勒展开,多项式牛顿迭代,多项式求逆,多项式开根,多项式ln,多项式exp,多项式快速幂这些非常重要基础的东西。
如果想要例题可以到我博客里面,多项式 和 生成函数 两个标签里面找,都配有题解。
【upd 2021.1.7】增加了 (exp) 的 (O(n^2)) 做法,这东西是有用处的。
【upd 2021.1.7】发现牛迭式子锅了,已经修改。
【upd 2021.1.18】增加了一个小Trick:CZT变换。因为多项式笔记(二)需要CZT变换这个前置知识。原本最好丢到笔记三的,但是怕鸽子太久,所以放到这里来了。
【upd 2021.1.19】从后面的笔记拉了多项式平移过来,觉得还是放在这里合适。
【upd 2021.1.19】NTT改成了预处理单位根的写法,实测洛谷多项式快速幂可以快1.2s。
【upd 2021.1.22】增加了多项式除法
【upd 2021.1.23】写多点求值的时候发现求逆写法过于不优良,已经更新(大概从直接NTT答案数组改成了copy过去,不然多点求值会直接RE)
泰勒展开
多项式求逆
一个很非常重要的思想。
已知 (F(x)*G(x)equiv 1pmod {x^n}) 和 (F(x)) ,求 (G(x))
假设我们已经求出了 (H(x)) 使得 (F(x)H(x)equiv 1 pmod{x^{lceilfrac{n}{2} ceil}})
又因为 (F(x)G(x)equiv 1pmod{x^{lceilfrac{n}{2} ceil}})
所以 (F(x)(G(x)-H(x))equiv 0pmod{x^{lceilfrac{n}{2} ceil}})
所以 (G(x)-H(x)equiv 0pmod{x^{lceilfrac{n}{2} ceil}})
两边平方可得 (G^2(x)+H^2(x)-2G(x)H(x)equiv 0pmod{x^n})
注意后面变成 (mod x^n) 了,这一步特别牛逼。
其实就是最高次变成了 (le n) 罢了,所以要 (mod x^n)
接着往下推
(G^2(x)+H^2(x)-2G(x)H(x)equiv 0pmod {x^n})
两边同乘 (F(x)) ,利用 (F(x)G(x)equiv 1pmod {x^n}) 这个性质
(G(x)+H^2(x)F(x)-2H(x)equiv 0pmod {x^n})
移项可得
(G(x)equiv 2H(x)-H^2(x)F(x)pmod {x^n})
得到了一个可以递归处理的东西,每次先把 (H(x)) 递归处理出来,在这一层FFT合并。
边界为 (G(0)=F(0)^{-1})
复杂度 (T(n)=T(dfrac{n}{2})+O(nlog n)=O(nlog n))
多项式牛顿迭代
考虑多项式求逆的过程,假设已经有 (F(x)*H(x)equiv 1pmod{x^{lceil frac{n}{2} ceil}}) ,有一个函数满足 (G(F(x))equiv 0pmod {x^n}) ,请记住这个设定
把函数 (G(F(x))) 在 (H(x)) 泰勒展开
从第三项开始,由于 (F(x),H(x)) 后 (lceildfrac{n}{2} ceil) 项相同,((F(x)-H(x))^2) 在 (mod x^n) 意义下为 (0) ,往后也都是 (0)
所以
多项式ln
两边求导
把 (A) 求逆,求导,得到的两个多项式卷起来,积分回去,得到 (B)
多项式exp
已知 (A(x)) ,求 (F(x)=e^{A(x)}) ,保证 (A(0)=0)
好像保证 (A(0)=0) 是为了保证exp有意义,具体原因等我知道了再补。
两边求导
令 (G(F(x))=ln (F(x))-A(x)) ,那么 (G'(F(x))=dfrac{1}{F(x)}) 。如果你不知道为什么 (A(x)) 可以当常数项消掉可以看看 某个洛谷帖子
对其牛顿迭代
用多项式求逆的思想,每次递归计算 (H(x)) 即可。
边界为 (G(0)=1)
复杂度 (T(n)=T(dfrac{n}{2})+O(nlog n)=O(nlog n))
【upd】(O(n^2)) 的做法,好写好推。有些题目 (n) 很小并且模数不一定为NTT模数就可以这么干。
提取 ([x^n]) 系数得
直接递推。
多项式开根
给定 (A(x)) ,求 (F^2(x)equiv A(x)pmod {x^n})
设 (G(F(x))=F^2(x)-A(x)) ,那么 (G'(F(x))=2F(x))
对其牛顿迭代
用多项式求逆的思想递归求 (H(x)) 即可。
(F(0) ot=1) 要二次剩余,不会,鸽子了。
所以边界(默认 (F(0)=1))为 (G(0)=1) 。
复杂度 (T(n)=T(dfrac{n}{2})+O(nlog n)=O(nlog n))
多项式快速幂
(ln + exp) 是非常强大的工具,结合起来用非常无敌。
这题是个经典例子,(ln +exp) 可以把幂运算转化成多项式的 +-*
运算。
求 (B(x)equiv A^k(x)pmod {x^n})
两边取 (ln) ,(ln (B(x))=kln (A(x)))
然后把 (ln (B(x))) (exp) 回去就可以得到 (B(x)) 了。
多项式除法
给定 (n) 次多项式 (F(x)) ,(m) 次多项式 (G(x)) ,求 (Q(x),R(x)) ,使得 (F(x)=Q(x)G(x)+R(x)) ,(Q(x)) 的次数为 (n-m) ,(R(x)) 的次数小于 (m) 。
如果我们想要翻转一个多项式的系数(reverse),可以采用如下方法:
把 (dfrac{1}{x}) 带到题目给出的等式里((F_r(x)) 表示把 (F(x)) 系数翻转的结果)
这样子就求出 (Q(x)) 了。(R(x)) 把 (Q(x)) 带回去就好了。
代码实现
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp(x,y) make_pair(x,y)
#define pb(x) push_back(x)
#define sz(v) (int)v.size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
const int N=100005;
const int M=N<<2;
#define mod 998244353
namespace math{
int inv[N],fac[N],ifc[N];
inline int qpow(int n,int k){int res=1;for(;k;k>>=1,n=1ll*n*n%mod)if(k&1)res=1ll*n*res%mod;return res;}
inline void fmod(int&x){x-=mod,x+=x>>31&mod;}
void initmath(const int&n=N-1){
inv[1]=1;for(int i=2;i<=n;++i)inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
fac[0]=1;for(int i=1;i<=n;++i)fac[i]=1ll*i*fac[i-1]%mod;
ifc[n]=qpow(fac[n],mod-2);for(int i=n-1;i>=0;--i)ifc[i]=1ll*(i+1)*ifc[i+1]%mod;
}
}
using namespace math;
namespace poly{
int rev[M],lg,lim,rt[M];
void initpoly(const int&n){
for(lg=0,lim=1;lim<=n;++lg,lim<<=1);
for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
for(int i=1;i<lim;i<<=1){
const int w=qpow(3,(mod-1)/(i<<1));
rt[i]=1;
for(int j=1;j<i;++j)rt[i+j]=1ll*rt[i+j-1]*w%mod;
}
}
void NTT(int*a,int op){
if(!op)reverse(a+1,a+lim);
for(int i=0;i<lim;++i)
if(i>rev[i])swap(a[i],a[rev[i]]);
for(int i=1;i<lim;i<<=1){
for(int j=0;j<lim;j+=i<<1){
for(int k=0;k<i;++k){
const int X=a[j+k],Y=1ll*rt[i+k]*a[i+j+k]%mod;
fmod(a[j+k]=X+Y),fmod(a[i+j+k]=X-Y+mod);
}
}
}
if(op)return;const int ilim=qpow(lim,mod-2);
for(int i=0;i<lim;++i)a[i]=1ll*a[i]*ilim%mod;
}
#define clr(a,n) memset(a,0,sizeof(int)*(n))
#define cpy(a,b,n) memcpy(a,b,sizeof(int)*(n))
void poly_mul(int*f,int*g,int*ans,int n,int m){
static int A[M],B[M];initpoly(n+m);
cpy(A,f,n),clr(A+n,lim-n),NTT(A,1);
cpy(B,g,m),clr(B+m,lim-m),NTT(B,1);
for(int i=0;i<lim;++i)A[i]=1ll*A[i]*B[i]%mod;
NTT(A,0),cpy(ans,A,n+m-1);
}
void poly_inv(int*g,int*f,int n){
static int A[M],B[M];
if(n==1)return g[0]=qpow(f[0],mod-2),void();
int m=(n+1)>>1;
poly_inv(g,f,m),initpoly(n<<1);
cpy(A,f,n),clr(A+n,lim-n),cpy(B,g,m),clr(B+m,lim-m);
NTT(A,1),NTT(B,1);
for(int i=0;i<lim;++i)A[i]=1ll*B[i]*(2-1ll*A[i]*B[i]%mod+mod)%mod;
NTT(A,0),cpy(g,A,n);
}
void poly_sqrt(int*g,int*f,int n){
static int A[M];
if(n==1)return g[0]=1,void();
poly_sqrt(g,f,(n+1)>>1);
clr(A,n),poly_inv(A,g,n),poly_mul(A,f,A,n,n);
for(int i=0,iv2=math::inv[2];i<n;++i)g[i]=1ll*(g[i]+A[i])*iv2%mod;
}
void jif(int*g,int*f,int n){
for(int i=1;i<=n;++i)g[i]=1ll*f[i-1]*math::inv[i]%mod;g[0]=0;
}
void dao(int*g,int*f,int n){
for(int i=0;i<n-1;++i)g[i]=1ll*(i+1)*f[i+1]%mod;g[n-1]=0;
}
void poly_ln(int*g,int*f,int n){
static int A[M],B[M];
dao(A,f,n),clr(B,n),poly_inv(B,f,n),poly_mul(A,B,A,n,n),jif(g,A,n);
}
void poly_exp(int*g,int*f,int n){
static int A[M];
if(n==1)return g[0]=1,void();
poly_exp(g,f,(n+1)>>1);
clr(A,n),poly_ln(A,g,n);
for(int i=0;i<n;++i)fmod(A[i]=(!i)+f[i]-A[i]+mod);
poly_mul(g,A,g,n,n),clr(g+n,lim-n);
}
void poly_qpow(int*g,int *f,int k,int n){
static int A[M];
clr(A,n),poly_ln(A,f,n);
for(int i=0;i<n;++i)A[i]=1ll*A[i]*k%mod;
poly_exp(g,A,n);
}
void poly_div(int*Q,int*R,int*F,int*G,int n,int m){
static int A[M];
reverse(F,F+n),reverse(G,G+m);
poly_inv(A,G,n-m+1),poly_mul(F,A,A,n-m+1,n-m+1);
for(int i=0;i<n-m+1;++i)Q[i]=A[n-m-i];
reverse(F,F+n),reverse(G,G+m),poly_mul(G,Q,A,m,n-m+1);
for(int i=0;i<m-1;++i)fmod(R[i]=F[i]+mod-A[i]);
}
}
-
主函数第一行一定要
initmath
!!! -
千万记得清空数组,太恶心了。
-
如何检查NTT是否正确:手造一组卷积。
-
如何检查多项式求逆是否正确:先求一次逆,再把得到的多项式求逆,看看是否是原来的。
-
如何检查多项式开根是否正确:先开根,然后把这个数组和自己卷一下,看看是否是原来的。
-
卡常卡傻了:NTT里面
i+j+k
改成i|j|k
慢0.2s,fmod
不加inline
慢 0.2s。。。 -
如果没有特殊说明,函数传参都是左闭右开的。
一些简单的技巧
都是平时做题遇上的
差卷积
求出 (C_k=sum_{i=k}^{n}A_iB_{i-k})
这时候把 (B) 反转为 (B') ,即 (B'_{n-i}=B_i) 。
这时候跑 (A*B') 就会发现 (C_k=sum_{i=0}^{k}A_iB'_{k-i}=sum_{i=0}^{k}A_iB_{n+i-k}) ,所以 (C_{i}=(A*B')_{n+i})
Chirp-Z Transform
一个非常牛逼的构造:
可以发现我们求的答案
差卷积即可。
如果每次快速幂计算是 (O(nlog mod)) 的,特别慢,最好直接两次前缀积(次数的二维前缀和)计算。
void CZT(int*g,int*f,int m,int n,int c){//给定n次多项式f,求f(c^{0->m-1})
static int pw[N<<1],ipw[N<<1],ivc,A[M],B[M];
pw[0]=pw[1]=ipw[0]=ipw[1]=1,ivc=qpow(c,mod-2);
for(int i=2;i<n+m;++i)pw[i]=1ll*pw[i-1]*c%mod,ipw[i]=1ll*ipw[i-1]*ivc%mod;
for(int i=2;i<n+m;++i)pw[i]=1ll*pw[i-1]*pw[i]%mod,ipw[i]=1ll*ipw[i-1]*ipw[i]%mod;
for(int i=0;i<n;++i)A[i]=1ll*f[i]*ipw[i]%mod;
for(int i=0;i<n+m;++i)B[i]=pw[i];
reverse(A,A+n),poly_mul(A,B,A,n,n+m);
for(int i=0;i<m;++i)g[i]=1ll*ipw[i]*A[i+n-1]%mod;
}
多项式平移
用 (F(x)) 得到 (F(x+c)) 的系数。
设 (F(x)=sumlimits_{i=0}^{n}a_ix^i) 。
差卷积即可。
参考资料
各个板题的洛谷题解,不一一挂连接了
广告:多项式笔记(二)