部分内容引自command_block的博客
多项式基础知识
多项式的表示法
系数表示法
$$f(x)=sumlimits_{i=0}^na_i*x^i$$
点值表示法
引理：$n+1$个不同的点可以唯一确定一个最高为$n$次的多项式
给出$n+1$维二元向量组$(vec{x},vec{y})$，满足$f(x_i)=y_i(forall i,j s.t. x_i
e x_j)$，确定一个多项式的形式
求值和插值
多项式求值：将系数表示转化为点值表示
多项式插值：将点值表示转化为系数表示
多项式插值及上述引理的解释：
多项式插值的本质是范德蒙矩阵求逆
$$egin{bmatrix}1&x_1&x_1^2&cdots&x_1^n\1&x_2&x_2^2&cdots&x_2^n\vdots&vdots&vdots&ddots&vdots\1&x_{n+1}&x^2_{n+1}&cdots&x^n_{n+1}end{bmatrix}egin{bmatrix}a_0\a_1\vdots\a_nend{bmatrix}=egin{bmatrix}y_1\y_2\vdots\y_{n+1}end{bmatrix}$$
可以证明若$x_i$各不相同则范德蒙矩阵一定可逆，则系数向量$vec{a}$有唯一解
多项式加减法
系数表示法
$$A(x)±B(x)=sumlimits_{i=0}^n(a_i±b_i)*x^i$$
点值表示法
$$(vec{x},vec{y_a})±(vec{x},vec{y_b})=(vec{x},vec{y_a}±vec{y_b})$$
多项式乘法
系数表示法
$$A(x)×B(x)=sumlimits_{i=0}^nsumlimits_{j=0}^na_i×b_j×x^{i+j}=sumlimits_{i=0}^{2n}sumlimits_{j=0}^{min(i,n)}a_j×b_{i-j}×x^i$$
点值表示法（向量$vec{x},vec{y_a},vec{y_b}$应为$2n+1$维以确定$2n$次多项式）
$$(vec{x},vec{y_a})*(vec{x},vec{y_b})=(vec{x},vec{y}),y_k=y_{a_k}y_{b_k}$$
卷积运算
$$c_k=sumlimits_{i⊕j=k}a_ib_j$$
其中⊕为某种运算
点值表示下的多项式乘法即为加法卷积
拉格朗日插值法
给定$n+1$个互不相同的点，确定$n$次多项式$f(x)$, 第$i$个点的坐标为$(x_i,y_i)$，我们需要找到该多项式在$k$点的取值，根据拉格朗日插值法有
$$f(k)=sumlimits_{i=0}^ny_iprodlimits_{i
e j}frac{k-x[j]}{x[i]-x[j]}$$
复杂度为$O(n^2)$
快速傅里叶变换$(FFT)$
基本思想：将系数表达转化为点值表达，利用点值表达相乘计算多项式乘法，再将点值表达相乘的结果转化为系数表达
将系数表达转化为点值表达$(DFT)$:
若$n$不是$2$的整数次幂，则补若干系数为$0$的项使$n$成为$2$的整数次幂
对$n$次多项式$f(x)=sumlimits_{i=0}^na_ix^i$，按$x$次数的奇偶性划分可得$f(x)=sumlimits_{i=0}^{2ile n}a_{2i}x^{2i}+sumlimits_{i=0}^{2i+1le n}a_{2i+1}x^{2i+1}$
令多项式
$f_1(x)=sumlimits_{i=0}^{2ile n}a_{2i}x^{i},f_2(x)=sumlimits_{i=0}^{2i+1le n}a_{2i+1}x^{i}$
则$f(x)=f_1(x^2)+xf_2(x^2)$
分别代入$w_{n+1}^k,w_{n+1}^{k+frac{n+1}{2}}$可得
$$f(w_{n+1}^k)=f_1(w_{frac{n+1}{2}}^k)+w_{n+1}^kf_2(w_{frac{n+1}{2}}^k)qquad(1)$$
$$f(w_{n+1}^{k+frac{n+1}{2}})=f_1(w_{frac{n+1}{2}}^k)-w_{n+1}^kf_2(w_{frac{n+1}{2}}^k)qquad(2)$$
所以如果我们知道两个多项式$f_1(x)$和$f_2(x)$分别在$w_{frac{n+1}{2}}^0...w_{frac{n+1}{2}}^{frac{n+1}{2}-1}$的点值表示
利用式$(1)$可求$f(x)$在$w_{n+1}^0...w_{n+1}^{frac{n+1}{2}-1}$的点值表示，利用式$(2)$可求$f(x)$在$w_{n+1}^{frac{n+1}{2}}...w_{n+1}^{n}$的点值表示，时间复杂度为$O(n)$
而对于$f_1(x),f_2(x)$也满足同样的性质，递归求解即可，时间复杂度$O(nlogn)$
单位根求解：$w_n^k=cosfrac{2pi k}{n}+i*sinfrac{2pi k}{n}$
将点值表达转化为系数表达$(IDFT)$
设$g=DFT(f)$即$g(k)=sumlimits_{i=0}^n(w_{n+1}^k)^if(i)$
则有$f=IDFT(g)$即$f(k)=frac{1}{n}sumlimits_{i=0}^n(w_{n+1}^{-k})^ig(i)$
挂一手代码，朴素递归写法（多项式乘法）
#include<iostream> #include<cstdio> #include<cmath> #define PAI acos(-1) using namespace std; const int maxn=4e6+10; struct complex{ double a,b; }f[maxn],g[maxn],mul[maxn]; int n,m,N; complex operator *(const complex &x,const complex &y) { return (complex){x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a}; } complex operator +(const complex &x,const complex &y) { return (complex){x.a+y.a,x.b+y.b}; } complex operator -(const complex &x,const complex &y) { return (complex){x.a-y.a,x.b-y.b}; } int up(int x) { int ret=1; for(;ret<x;ret*=2); return ret; } complex makew(int len,int k) { return (complex){cos(2*PAI*k/len),sin(2*PAI*k/len)}; } void fft(int ty,int len,complex *a) { if(len==1) return; complex f1[len/2],f2[len/2]; for(int i=0;2*i<len;i++) f1[i]=a[2*i]; for(int i=0;2*i+1<len;i++) f2[i]=a[2*i+1]; fft(ty,len/2,f1),fft(ty,len/2,f2); for(int i=0;i<len/2;i++) { complex w=makew(len,ty*i); a[i]=f1[i]+w*f2[i]; a[i+len/2]=f1[i]-w*f2[i]; } } int main() { scanf("%d%d",&n,&m); N=up(n+m+1); for(int i=0;i<=n;i++) scanf("%lf",&f[i].a); for(int i=0;i<=m;i++) scanf("%lf",&g[i].a); fft(1,N,f),fft(1,N,g); for(int i=0;i<N;i++) mul[i]=f[i]*g[i]; fft(-1,N,mul); for(int i=0;i<n+m+1;i++) printf("%d ",(int)(mul[i].a/N+0.5)); printf(" "); return 0; }
蝴蝶操作的引入与$FFT$的实现优化
考虑$FFT$的递归与回溯过程
| $Level1$ | $0$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ |
| $Level2$ | $0$ | $2$ | $4$ | $6$ | $1$ | $3$ | $5$ | $7$ |
| $Level3$ | $0$ | $4$ | $2$ | $6$ | $1$ | $5$ | $3$ | $7$ |
| $Level4$ | $0$ | $4$ | $2$ | $6$ | $1$ | $5$ | $3$ | $7$ |
我们在递归写法中对每一层存储相应的值，时空代价高昂
但可以发现$FFT$的实现具有这样的性质：
上层值传递给下层后不再被使用，下层值生成上层后不再被使用
$FFT$的实现优化：基于这种性质我们不妨将原系数数组通过调整顺序转换为递归底层的顺序从而省去实际分配传递序列值的过程，然后自底向上按照原合并过程进行合并，这一改进使得我们不再依赖递归完成$FFT$的实现，从而省去了不必要的时空开支
蝴蝶定理：$FFT$递归底层元素序号是其自身位置的二进制反序
| $Top$ | $000$ | $001$ | $010$ | $011$ | $100$ | $101$ | $110$ | $111$ |
| $Bottom$ | $000$ | $100$ | $010$ | $110$ | $001$ | $101$ | $011$ | $111$ |
二进制反序预处理
for(int i=0;i<N;i++) tr[i]=(tr[i>>1]>>1)|((i&1)*(N>>1)); 理解： ~~~~~ x 其他部分 最后一位 翻转的话就相当于 x 连接上 其他部分的反转 其他部分的反转=tr[i>>1]>>1 然后判一下最后一位,如果是1则加上N>>1
完全版$FFT$代码
#include<iostream> #include<cstdio> #include<cmath> using namespace std; const int maxn=3e6+10; const double PAI=acos(-1); struct complex{ double a,b; }f[maxn],g[maxn],mul[maxn]; int n,m,N,tr[maxn]; complex operator *(const complex &x,const complex &y) { return (complex){x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a}; } complex operator +(const complex &x,const complex &y) { return (complex){x.a+y.a,x.b+y.b}; } complex operator -(const complex &x,const complex &y) { return (complex){x.a-y.a,x.b-y.b}; } int up(int x) { int ret=1; for(;ret<x;ret*=2); return ret; } complex makew(int len,int k) { return (complex){cos(2*PAI*k/len),sin(2*PAI*k/len)}; } void fft(int ty,complex *a) { for(int i=0;i<N;i++) if(i<tr[i]) swap(a[i],a[tr[i]]); for(int len=2;len<=N;len<<=1) { complex stw=makew(len,ty); for(int pos=0;pos<N;pos+=len) { complex wn=(complex){1,0}; for(int i=0;i<len/2;i++,wn=wn*stw) { complex w=wn*a[i+pos+len/2]; a[i+pos+len/2]=a[i+pos]-w; a[i+pos]=a[i+pos]+w; } } } } int main() { scanf("%d%d",&n,&m); N=up(n+m+1); for(int i=0;i<=n;i++) scanf("%lf",&f[i].a); for(int i=0;i<=m;i++) scanf("%lf",&g[i].a); for(int i=0;i<N;i++) tr[i]=(tr[i>>1]>>1)|((i&1)*(N>>1)); fft(1,f),fft(1,g); for(int i=0;i<N;i++) mul[i]=f[i]*g[i]; fft(-1,mul); for(int i=0;i<n+m+1;i++) printf("%d ",(int)(mul[i].a/N+0.5)); printf(" "); return 0; }
数论变换$(NTT)$
$FFT$的不足：必须做复数而且是浮点数的运算，因此计算量会比较大，而且浮点数运算产生的误差会比较大
关于数论意义下$FFT$算法优化的可能性探讨：
对于质数$p=qk+1(k=2^m)$，设其原根为$g$，我们用$G_n$来代替$w_n$，其中$G_n=g^{frac{p-1}{n}}$，由原根的性质（见数论内容简要整理）可得$G_n$满足如下性质：
$1,$任取$0leqslant i,j<n$有$G_n^i
eq G_n^j$
$2,G_{2n}^{2k}=G_n^k$
$3,G_n^{k+frac{n}{2}}=-G_n^k$
因此在数论意义下采用原根代替单位根是可行的
局限性：仅在模$p$意义下有效
常用质数及其原根：$998244353,1004535809,469762049$原根均为$3$
代码
#include<iostream> #include<cstdio> #define int long long using namespace std; const int maxn=3e6+10,gs=3,p=998244353; int n,m,N,tr[maxn],f[maxn],g[maxn],mul[maxn]; int up(int x) { int ret=1; for(;ret<x;ret*=2); return ret; } int qpow(int x,int y) { int ret=1; for(;y;y>>=1,(x*=x)%=p) if(y&1) (ret*=x)%=p; return ret; } inline void mod(int &x) { x=(x%p+p)%p; } void ntt(bool ty,int *a) { for(int i=0;i<N;i++) if(i<tr[i]) swap(a[i],a[tr[i]]); for(int len=2;len<=N;len<<=1) { int stg=qpow(gs,(p-1)/len); if(ty) stg=qpow(stg,p-2); for(int pos=0;pos<N;pos+=len) for(int gn=1,i=0;i<len/2;i++,(gn*=stg)%=p) { int G=gn*a[i+pos+len/2]%p; a[i+pos+len/2]=a[i+pos]-G; a[i+pos]=a[i+pos]+G; mod(a[i+pos+len/2]),mod(a[i+pos]); } } } signed main() { scanf("%lld%lld",&n,&m); N=up(n+m+1); for(int i=0;i<=n;i++) scanf("%lld",&f[i]); for(int i=0;i<=m;i++) scanf("%lld",&g[i]); for(int i=0;i<N;i++) tr[i]=(tr[i>>1]>>1)|((i&1)*(N>>1)); ntt(0,f),ntt(0,g); for(int i=0;i<N;i++) mul[i]=f[i]*g[i]%p; ntt(1,mul); for(int i=0;i<n+m+1;i++) printf("%lld ",mul[i]*qpow(N,p-2)%p); printf(" "); return 0; }
多项式求逆
问题描述：给定多项式$f(x)$求$f^{-1}(x)$
假设我们已经求出$f(x)$在模$x^{lceilfrac{n}{2}
ceil}$意义下的逆元$f_0^{-1}(x)$，则有
$$f(x)f^{-1}_0(x)equiv 1(mod x^{lceilfrac{n}{2}
ceil})$$
$$f(x)f^{-1}(x)equiv 1(mod x^{lceilfrac{n}{2}
ceil})$$
$$f^{-1}(x)-f^{-1}_0(x)equiv 0(mod x^{lceilfrac{n}{2}
ceil})$$
两边平方可得
$$f^{-2}(x)-2f^{-1}(x)f^{-1}_0(x)+f^{-2}_0(x)equiv 0(mod x^n)$$
两边同乘$f(x)$并移项可得
$$f^{-1}(x)equiv 2f^{-1}_0(x)-f(x)f^{-2}_0(x)(mod x^n)$$
递归计算即可
后记：感觉这个算法还是挺好理解的，但是实现过程中有很多需要注意的地方，有助于加深对多项式的理解
1，直接将$f(x)$转换为点值表达分别对各点求逆（取倒数），再转化为系数表达，由于多项式的逆是无穷级数，在模意义下高位累加到低位上，而点值表达求逆只能保证求出的多项式与原多项式在选取的点上满足互为倒数
2，在点值表示下模拟上述过程，最后转为系数表示，假设我们已经求得在$mod x^{lceilfrac{n}{2}
ceil}$意义下的$f_0^{-1}(x)$，我们取$mod x^{lceilfrac{n}{2}
ceil}$意义下的$f(x)$，由上式可知，直接使用点值表示相乘的$f^{-1}(x)$是一个高于$n$次的多项式，用这一多项式作$f^{-1}_0(x) mod x^n$继续运算自然也是错误的
3，由$2$，直接使用点值表示相乘的$f^{-1}(x)$是一个高于$n$次的多项式，取点个数应设为$2N$，只取$N$个点无法正确转化为系数表达
最后放上代码
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define int long long using namespace std; const int maxn=5e5+10,gs=3,p=998244353; int n,N,tr[maxn],f[maxn],f0[maxn],g[maxn],invN; int up(int x) { int ret=1; for(;ret<x;ret<<=1); return ret; } int qpow(int x,int y) { int ret=1; for(;y;y>>=1,(x*=x)%=p) if(y&1) (ret*=x)%=p; return ret; } inline void mod(int &x) { x=(x%p+p)%p; } void ntt(bool ty,int *a) { for(int i=0;i<N;i++) if(i<tr[i]) swap(a[i],a[tr[i]]); for(int len=2;len<=N;len<<=1) { int stg=qpow(gs,(p-1)/len); if(ty) stg=qpow(stg,p-2); for(int pos=0;pos<N;pos+=len) for(int gn=1,i=0;i<len/2;i++,(gn*=stg)%=p) { int G=gn*a[i+pos+len/2]%p; a[i+pos+len/2]=a[i+pos]-G; a[i+pos]=a[i+pos]+G; mod(a[i+pos+len/2]),mod(a[i+pos]); } } } int fx(int x) { int ret=0; for(int i=0,xp=1;i<N;i++,(xp*=x)%=p) (ret+=f[i]*xp)%=p; return ret; } void inv(int pw) { if(pw==1) { memset(g,0,sizeof(g)); g[0]=qpow(f[0],p-2); return; } inv((pw+1)/2); N=up(pw)<<1,invN=qpow(N,p-2); for(int i=0;i<N;i++) tr[i]=(tr[i>>1]>>1)|((i&1)*(N>>1)); for(int i=0;i<pw;i++) f0[i]=f[i]; for(int i=pw;i<N;i++) f0[i]=0; ntt(0,g),ntt(0,f0); for(int i=0;i<N;i++) g[i]=2*g[i]-g[i]*g[i]%p*f0[i]%p,mod(g[i]); ntt(1,g); for(int i=0;i<pw;i++) (g[i]*=invN)%=p; for(int i=pw;i<N;i++) g[i]=0; } signed main() { scanf("%lld",&n); for(int i=0;i<n;i++) scanf("%lld",&f[i]); inv(n); for(int i=0;i<n;i++) printf("%lld ",g[i]); printf(" "); return 0; }
分治$FFT$
多项式求逆解法
设$g_0=0$，则有$f_i=sumlimits_{j+k=i}f_jg_k$
不妨设$F(x)=sumlimits_{i=0}^{n-1}f_ix^i,G(x)=sumlimits_{i=0}^{n-1}g_ix^i(mod x^n)$
则有$F(x)G(x)equivsumlimits_{i=0}^{n-1}x^isumlimits_{j+k=i}f_jg_kequiv F(x)-1(mod x^n)$
即$F(x)equiv (1-G(x))^{-1}(mod x^n)$，那么就是一个多项式求逆的模板了
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define int long long using namespace std; const int maxn=5e5+10,gs=3,p=998244353; int n,N,tr[maxn],f[maxn],f0[maxn],g[maxn],invN; int up(int x) { int ret=1; for(;ret<x;ret<<=1); return ret; } int qpow(int x,int y) { int ret=1; for(;y;y>>=1,(x*=x)%=p) if(y&1) (ret*=x)%=p; return ret; } inline void mod(int &x) { x=(x%p+p)%p; } void ntt(bool ty,int *a) { for(int i=0;i<N;i++) if(i<tr[i]) swap(a[i],a[tr[i]]); for(int len=2;len<=N;len<<=1) { int stg=qpow(gs,(p-1)/len); if(ty) stg=qpow(stg,p-2); for(int pos=0;pos<N;pos+=len) for(int gn=1,i=0;i<len/2;i++,(gn*=stg)%=p) { int G=gn*a[i+pos+len/2]%p; a[i+pos+len/2]=a[i+pos]-G; a[i+pos]=a[i+pos]+G; mod(a[i+pos+len/2]),mod(a[i+pos]); } } } int fx(int x) { int ret=0; for(int i=0,xp=1;i<N;i++,(xp*=x)%=p) (ret+=f[i]*xp)%=p; return ret; } void inv(int pw) { if(pw==1) { memset(g,0,sizeof(g)); g[0]=qpow(f[0],p-2); return; } inv((pw+1)/2); N=up(pw)<<1,invN=qpow(N,p-2); for(int i=0;i<N;i++) tr[i]=(tr[i>>1]>>1)|((i&1)*(N>>1)); for(int i=0;i<pw;i++) f0[i]=f[i]; for(int i=pw;i<N;i++) f0[i]=0; ntt(0,g),ntt(0,f0); for(int i=0;i<N;i++) g[i]=2*g[i]-g[i]*g[i]%p*f0[i]%p,mod(g[i]); ntt(1,g); for(int i=0;i<pw;i++) (g[i]*=invN)%=p; for(int i=pw;i<N;i++) g[i]=0; } signed main() { scanf("%lld",&n); f[0]=1; for(int i=1;i<n;i++) scanf("%lld",&f[i]),f[i]*=-1; inv(n); for(int i=0;i<n;i++) printf("%lld ",g[i]); printf(" "); return 0; }