多项式学习笔记
Post time: 2021-03-11 18:34:53
从今日起开始学习多项式所有基本操作，并加以部分练习。
UPDATE 2022.3.17 结尾放上了正在更新的多项式模板。
一、多项式乘法（快速傅里叶变换FFT）
1.1 定义
快速傅里叶变换（fast Fourier transform，FFT），是用来计算离散傅里叶变换（discrete Fourier transform，DFT）及其逆变换（IDFT）的快速算法。
1.2 直观感受
DFT相当于是把时域信号（多项式系数表示）转化为频域信号（多项式点值表示）的过程，IDFT相反。
另外不懂的一句话：多项式乘法实际上是多项式系数向量的卷积。
1.3 单位根
点值表示法对应的单位根是 \(2n\) 次单位根。所谓“\(n\) 次单位根”，指满足 \(x^n=1\) 的复数。在复数域中，\(1\) 有 \(n\) 个单位根 \(e^{\frac{2k\pi i}{n}}\)，其中 \(i^2=-1\)。
在复数平面上，设 \(n\) 次单位根为 \(\omega_n^1,\omega_n^2,\ldots,\omega_n^n\)，则有 \(\omega_{n}^{k}=\cos\ k \frac{2\pi}{n}+i\sin k\frac{2\pi}{n}\)。它们有这样几个性质：
性质1 \(\ \omega_{2n}^{2k}=\omega_n^k\)；
性质2 \(\ \omega_n^{\frac{n}{2}}=-1\)；
性质3 \(\ \omega_n^{k+\frac{n}{2}}=-\omega_n^k\)。
1.4 离散傅里叶变换(DFT)
设一个 \(n-1\) 次多项式 \(A\) 的 \(n\) 个系数为 \(a_0,a_1,\ldots,a_{n-1}\)，这里默认 \(n\) 为 \(2\) 的整数次幂。
则
设
不难发现，\(A(x)=A_0(x^2)+xA_1(x^2)\)。
在多项式中带入单位根 \(\omega_n^k\) 和 \(\omega_n^{k+\frac{n}{2}}(k<\frac{n}{2})\)，得
不难发现，因为只有常数项不同，所以只要我们求出了 \(\omega_n^k(k<\frac{n}{2})\)，就可以 \(O(1)\) 得出 \(\omega_n^{k+\frac{n}{2}}\)，从而可以将问题递归求解，复杂度 \(O(n\log n)\)。
1.5 离散傅里叶逆变换(IDFT)
这里需要把点值表示（设为 \((y_0,y_1,\ldots,y_{n-1})\)）转化为系数表示。
设一组向量为 \((c_0,c_1,\ldots,c_{n-1})\)，满足：
即 \(B(x)=y_0+y_1x+y_2x^2+\ldots+y_{n-1}x^{n-1}\) 在 \(\omega_n^0,\omega_n^{-1},\ldots,\omega_n^{-(n-1)}\) 处的点值表示，则有：
设 \(S(x)=\sum\limits_{i=0}^{n-1}x^i\)，则
显然，当 \(k\neq 0\) 时，\(S(\omega_n^k)=0\)；当 \(k=0\) 时， \(S(\omega_n^0=1)=n\)。
代回原式，可得
于是，递归求解的FFT所有理论都已经证毕。但是有些（甚至是大部分）时候，递归版FFT由于需要重新开数组，导致常数非常大，会被卡。所以需要研究一种迭代版本的FFT。
1.6 迭代版FFT
考虑如何优化之前的递归版做法。可以首先举一个简单的例子：
\(n=8\) 时的递归情况：
step 1: 0,1,2,3,4,5,6,7
step 2: 0,2,4,6; 1,3,5,7
step 3: 0,4; 2,6; 1,5; 3,7;
step 4: 0; 4; 2; 6; 1; 5; 3; 7;
乍一看好像没什么规律，那么写出它们的二进制表示：
000,100,010,110,001,101,011,111
再写出它们反转后的二进制表示：
000,001,010,011,100,101,110,111
惊奇地发现它是递增的！
所以，可以将原数组按照二进制反转后表示的数的大小进行排序，从而得到一个最终的操作序列，这样就可以直接进行操作了。
（然后就咕了很久
二、快速数论变换NTT
将FFT板子中的单位根换成原根，即可得到NTT。由于复数和浮点运算带来的不利影响，FFT比NTT慢一些。
原根的定义大概是说：若一个数 \(g\) 是 \(P\) 的原根，则需要满足 \(g^i(0\leq i<n) \mod P\) 互不相等。原根拥有所有上述FFT中需要用到的单位根的性质，所以可以直接替换。
一般 \(P\) 取 \(998244353\)，\(g\) 取 \(3\)，逆元为 \(332748118\)。
三、任意模数多项式乘法
NTT+CRT，只要使CRT最后那个式子模的那个 \(\prod p_i>\) 答案最大值即可。一般就选三个模数 \(998244353,469762049,1004535809\)，它们的原根都是 \(3\)。
点击查看代码
#include<iostream>
#include<cstdio>
using namespace std;
inline int qpow(int a,int k,int p){int s=1;for(;k;k>>=1,a=1ll*a*a%p)if(k&1)s=1ll*s*a%p;return s;}
const int N=4e5+13,g=3;
const int p[4]={0,469762049,998244353,1004535809};
const int gi[4]={0,qpow(g,p[1]-2,p[1]),qpow(g,p[2]-2,p[2]),qpow(g,p[3]-2,p[3])};
int n,m,P,a[N],b[N],aa[N],bb[N],r[N],limit=1,ans[N][4];
__int128 Inv[4][4],MUL;
inline void init(){
for(int i=1;i<=3;++i)
for(int j=1;j<=3;++j)
if(i!=j) Inv[i][j]=qpow(p[i]%p[j],p[j]-2,p[j]);
}
inline void fft(int *f,int type,int op){
for(int i=0;i<limit;++i)
if(i<r[i]) swap(f[i],f[r[i]]);
for(int mid=1;mid<limit;mid<<=1){
int Wn=qpow(type==1?g:gi[op],(p[op]-1)/(mid<<1),p[op]);
for(int j=0;j<limit;j+=(mid<<1)){
int w=1;
for(int k=0;k<mid;++k,w=1ll*w*Wn%p[op]){
int x=f[j+k],y=1ll*w*f[j+k+mid]%p[op];
f[j+k]=(x+y)%p[op];f[j+k+mid]=(x-y+p[op])%p[op];
}
}
}
}
inline void get(int x){
__int128 res=0,a=ans[x][1],b=ans[x][2],c=ans[x][3],A=p[1],B=p[2],C=p[3];
__int128 M=A*B*C;
res+=a*B%M*C%M*Inv[2][1]%M*Inv[3][1]%M;res%=M;
res+=b*A%M*C%M*Inv[1][2]%M*Inv[3][2]%M;res%=M;
res+=c*A%M*B%M*Inv[1][3]%M*Inv[2][3]%M;res%=M;
printf("%lld ",(long long)(res%P));
}
int main(){
init();
scanf("%d%d%d",&n,&m,&P);
for(int i=0;i<=n;++i) scanf("%d",&aa[i]);
for(int i=0;i<=m;++i) scanf("%d",&bb[i]);
int l=0;
while(limit<=n+m) limit<<=1,++l;
for(int i=0;i<limit;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
for(int t=1;t<=3;++t){
for(int i=0;i<limit;++i) a[i]=aa[i];
for(int i=0;i<limit;++i) b[i]=bb[i];
fft(a,1,t);fft(b,1,t);
for(int i=0;i<limit;++i) a[i]=1ll*a[i]*b[i]%p[t];
fft(a,-1,t);int inv=qpow(limit,p[t]-2,p[t]);
for(int i=0;i<=n+m;++i) ans[i][t]=1ll*a[i]*inv%p[t];
}
for(int i=0;i<=n+m;++i) get(i);
return 0;
}
四、多项式乘法逆
给定 \(F(x)\)，求使 \(F(x)\times G(x)\equiv 1\pmod{x^n}\) 成立的 \(G(x)\) 的各项系数。
首先考虑如果已经求得了 \(\mod{x^{\lceil\frac n2\rceil}}\) 的答案 \(G'(x)\)，如何推出 \(G(x)\)？
可得
于是可以利用这个性质倍增求解，复杂度 \(O(n\log n)\)。
五、任意模数多项式乘法逆
仍然像任意模数NTT一样，考虑 三模NTT+CRT 来解决。注意到任意模数多项式乘法需要满足的是最大可能的值小于几个模数的乘积（因为这样才能保证CRT最后那个式子有唯一解）。
考虑这个乘法逆的式子：
显然，直接计算 \(FG'^2\) 会超上限。所以可以先算一次 \(FG'\)，然后再算 \(G\) 即可。
虽然总的时间复杂度为 \(O(n\log n)\)，但是由于每一层计算 \(G\) 需要 \(18\) 次 NTT，常数可想而知。
点击查看代码
#include<iostream>
#include<cstdio>
using namespace std;
typedef __int128 iint;
const int N=4e5+13,p1=998244353,p2=469762049,p3=1004535809,P=1e9+7,G=3;
inline int qpow(int a,int k,int p){int s=1;for(;k;k>>=1,a=1ll*a*a%p)if(k&1)s=1ll*s*a%p;return s;}
inline int inv(int a,int p){return qpow(a,p-2,p);}
const int g1=inv(G,p1),g2=inv(G,p2),g3=inv(G,p3);
int INV[4][4];
int n,a[N],b[N],r[N],c[N],aa[N],b1[N],b2[N],b3[N],d[N];
inline void init(){
INV[1][2]=inv(p1%p2,p2);INV[1][3]=inv(p1%p3,p3);
INV[2][1]=inv(p2%p1,p1);INV[2][3]=inv(p2%p3,p3);
INV[3][1]=inv(p3%p1,p1);INV[3][2]=inv(p3%p2,p2);
}
inline void fft(int *f,int limit,int type,int p,int g){
for(int i=0;i<limit;++i)
if(i<r[i]) swap(f[i],f[r[i]]);
for(int i=1;i<limit;i<<=1){
int Wn=qpow(type==1?G:g,(p-1)/(i<<1),p);
for(int j=0;j<limit;j+=(i<<1)){
int w=1;
for(int k=0;k<i;++k,w=1ll*w*Wn%p){
int x=f[j+k],y=1ll*w*f[j+k+i]%p;
f[j+k]=(x+y)%p,f[j+k+i]=(x-y+p)%p;
}
}
}
if(type==1) return;
int iinv=inv(limit,p);
for(int i=0;i<limit;++i) f[i]=1ll*f[i]*iinv%p;
}
inline int crt(iint a,iint b,iint c){
iint res=0,A=p1,B=p2,C=p3,M=A*B*C;
res+=a*B%M*C%M*INV[2][1]%M*INV[3][1]%M;res%=M;
res+=b*A%M*C%M*INV[1][2]%M*INV[3][2]%M;res%=M;
res+=c*A%M*B%M*INV[1][3]%M*INV[2][3]%M;res%=M;
return (int)(res%P);
}
inline void mul(int *a,int *b,int limit,int p,int g){
fft(a,limit,1,p,g),fft(b,limit,1,p,g);
for(int i=0;i<limit;++i) a[i]=1ll*a[i]*b[i]%p;
fft(a,limit,-1,p,g);
}
inline void Mul(int *a,int *b,int *A,int limit){
for(int i=0;i<limit;++i) aa[i]=a[i],b1[i]=b[i];
mul(b1,aa,limit,p1,g1);
for(int i=0;i<limit;++i) aa[i]=a[i],b2[i]=b[i];
mul(b2,aa,limit,p2,g2);
for(int i=0;i<limit;++i) aa[i]=a[i],b3[i]=b[i];
mul(b3,aa,limit,p3,g3);
for(int i=0;i<limit;++i) A[i]=crt(b1[i],b2[i],b3[i]);
}
inline void Inv(int n,int *a,int *b){
b[0]=inv(a[0],P);
int l=1;
while(l<(n<<1)){
for(int i=0;i<(l<<1);++i) r[i]=(r[i>>1]>>1)|((i&1)?l:0);
for(int i=0;i<l;++i) c[i]=a[i];
for(int i=l;i<(l<<1);++i) c[i]=0;
Mul(c,b,d,l<<1);
for(int i=1;i<(l<<1);++i) c[i]=(-d[i]+P)%P;c[0]=1;
Mul(c,b,d,l<<1);
for(int i=0;i<l;++i) b[i]=d[i];
for(int i=l;i<(l<<1);++i) b[i]=0;
l<<=1;
}
}
int main(){
init();
scanf("%d",&n);
for(int i=0;i<n;++i) scanf("%d",&a[i]);
Inv(n,a,b);
for(int i=0;i<n;++i) printf("%d ",b[i]);
return 0;
}
六、多项式对数函数（多项式 ln）
求一个多项式 \(B(x)\) 满足 \(B(x)\equiv \ln A(x)\pmod{x^n}\)。
设 \(f(x)=\ln x\)，则
求导，得
这样可以很显然地得到，要从 \(A\) 求 \(B\)，即用 \(A\) 的导数乘 \(A\) 的逆元再积分回去即可。
附上求导求积分代码：
点击查看代码
inline void getdet(int n,int *f,int *g){
for(int i=1;i<n;++i) g[i-1]=1ll*i*f[i]%P;
}
inline void getinvdet(int n,int *f,int *g){
for(int i=0;i<n-1;++i) g[i+1]=1ll*inv(i+1)*f[i]%P;
g[0]=0;
}
七、多项式指数函数（多项式exp）
前置：泰勒展开
如果要求一个函数 \(f(x)\) 某一点上的值，我们可以构造一个函数 \(g(x)\)，那么对于 \(x\) 点的值，有 \(f(x)\approx g(x)=g(x_0)+\frac{f^1(x_0)}{1!}(x-x_0)+\frac{f^2(x_0)}{2!}(x-x_0)^2+……+\frac{f^n(x_0)}{n!}(x-x_0)^n\)
其中 \(f^n(x_0)\) 表示对原函数图像上 \(x_0\) 这个点进行 \(n\) 阶求导。
令 \(G'(x)\equiv e^{F(x)}\pmod {x^{\lceil\frac{n}{2}\rceil}}\)，\(G(x)\equiv e^{F(x)}\pmod {x^n}\)。对于 \(\ln G(x)\) 在 \(x_0=G'(x)\) 处泰勒展开，可得
神奇的是，由于 \(G(x)-G'(x)\equiv 0\pmod {x^{\lceil\frac{n}{2}\rceil}}\)，从 \((G(x)-G'(x))^2\) 开始，由于 \((G(x)-G'(x))^i\) 指数 \(\geq 2\) 所以 \((G(x)-G'(x))^i\equiv 0\pmod{x^n}\)，也就是后面的项全没了！
化一化式子就可以得到，\(G(x)\equiv G'(x)+(F(x)-\ln G'(x))G'(x)\pmod {x^n}\)。复杂度 \(O(n\log n)\)。
八、多项式除法
给定 \(n\) 次多项式 \(A(x)\) 和 \(m\) 次多项式 \(B(x)\)，求满足 \(A(x)=B(x)P(x)+Q(x)\) 的 \(n-m\) 次多项式 \(P(x)\) 以及 \(m-1\) 次多项式 \(Q(x)\)。相当于除法中的商和余数。
一个很 nb 的转化：考虑把等号两边所有系数全反过来，变成 \(A'(x)=B'(x)P'(x)+Q'(x)x^{n-m+1}\)。然后在 \(\bmod x^{n-m+1}\) 意义下，后面这一项直接没了，也就是说可以直接多项式求逆求出 \(P'(x)\)，转为 \(P(x)\) 之后，再用 \(Q(x)=A(x)-B(x)P(x)\) 求出 \(Q(x)\) 即可。注意到 \(P(x)\) 是 \(n-m\) 次多项式，所以在 \(\bmod x^{n-m+1}\) 意义下不会丢失任何信息。复杂度 \(O(n\log n)\)。
九、多项式多点求值
求 \(F(x)\) 在 \(x_1,x_2,\ldots,x_k\) 处的点值。
注意到 \(F(x)\bmod (x-x_1)(x-x_2)\ldots(x-x_k)\) 之后的值在 \(x_1,x_2,\ldots,x_k\) 处的点值是不变的（\(F(x)=P(x)R(x)+Q(x)\)，\(x=x_i\) 的时候 \(R(x)=0\)，所以 \(F(x)=Q(x)\)），所以考虑分治求解，对一段区间求点值先对整个区间的 \(\prod (x-x_i)\) 取模，这样可以把次数限制在区间长度内，于是得到了复杂度 \(O(n\log^2 n)\) 的做法。
十、多项式快速插值
考虑拉格朗日插值式子
对于每个 \(i\)，用极限的方法可以这样表示：
然后这一步是我不太懂的洛必达：
计算 \(\sum d_i\prod_{j\not= i} (x-x_j)\)，直接分治 ntt，每次记录两个值 \(A(x)=\sum d_i\prod (x-x_i),B(x)=\prod(x-x_i)\)，然后合并区间的时候 \(newA=A_lB_r+A_rB_l\) 即可。复杂度是常数很大的 \(O(n\log^2 n)\)，一般 \(50000\) 都要跑 \(2\sim 3s\)。
My Poly
点击查看代码
namespace poly{
const int G=3,Gi=qpow(G,mod-2);
inline void poly_det(int n,int *a,int *b){for(int i=1;i<n;++i) b[i-1]=(ll)i*a[i]%mod;b[n-1]=0;}
inline void poly_int(int n,int *a,int *b){for(int i=0;i<n-1;++i) b[i+1]=(ll)ni[i+1]*a[i]%mod;b[0]=0;}
int r[N];
inline void fft(int *f,int limit,int type){
for(int i=0;i<limit;++i)
if(i<r[i]) swap(f[i],f[r[i]]);
for(int i=1;i<limit;i<<=1){
int Wn=qpow(type==1?G:Gi,(mod-1)/(i<<1));
for(int j=0;j<limit;j+=(i<<1)){
int w=1;
for(int k=0;k<i;++k,w=(ll)w*Wn%mod){
int x=f[j+k],y=(ll)w*f[j+k+i]%mod;
f[j+k]=md(x+y),f[j+k+i]=md(x-y+mod);
}
}
}
if(type==1) return;
int Inv=qpow(limit,mod-2);
for(int i=0;i<limit;++i) f[i]=(ll)f[i]*Inv%mod;
}
inline void poly_mul(int n,int *a,int *b){
static int mul_a[N],mul_b[N];
#define x mul_a
#define y mul_b
int limit=1,l=0;
while(limit<n) limit<<=1,++l;
for(int i=0;i<limit;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
for(int i=0;i<limit;++i) x[i]=a[i],y[i]=b[i];
fft(x,limit,1),fft(y,limit,1);
for(int i=0;i<limit;++i) x[i]=(ll)x[i]*y[i]%mod;
fft(x,limit,-1);
for(int i=0;i<n;++i) a[i]=x[i];
#undef x
#undef y
}
inline void poly_inv(int n,int *a,int *b){
static int inv_a[N];
#define c inv_a
for(int i=0,lim=min(N-1,n<<2);i<lim;++i) b[i]=0;
b[0]=qpow(a[0],mod-2);int l=2;
while(l<(n<<1)){
for(int i=0;i<(l<<1);++i) r[i]=(r[i>>1]>>1)|((i&1)?l:0);
for(int i=0;i<l;++i) c[i]=a[i];
for(int i=l;i<(l<<1);++i) c[i]=0;
fft(c,l<<1,1),fft(b,l<<1,1);
for(int i=0;i<(l<<1);++i) b[i]=(ll)md(2-(ll)b[i]*c[i]%mod+mod)*b[i]%mod;
fft(b,l<<1,-1);
for(int i=l;i<(l<<1);++i) b[i]=0;
l<<=1;
}
#undef c
}
inline void poly_ln(int n,int *a,int *b){
assert(a[0]==1);
static int ln_a[N],ln_b[N];
#define c ln_a
#define d ln_b
poly_det(n,a,c);
poly_inv(n,a,d);
poly_mul(n<<1,c,d);
poly_int(n,c,b);
#undef c
#undef d
}
inline void poly_exp(int n,int *a,int *b){
assert(a[0]==0);
static int exp_a[N],exp_b[N];
#define c exp_a
#define d exp_b
for(int i=0,lim=min(N-1,n<<2);i<lim;++i) b[i]=0;
b[0]=1;int l=2;
while(l<(n<<1)){
poly_ln(l,b,d);
for(int i=0;i<l;++i) c[i]=md(a[i]-d[i]+mod);modadd(c[0],1);
for(int i=l;i<(l<<1);++i) c[i]=0;
poly_mul(l<<1,b,c);
for(int i=l;i<(l<<1);++i) b[i]=0;
l<<=1;
}
#undef c
#undef d
}
inline void poly_div(int n,int m,int *a,int *b,int *p,int *q){
static int div_a[N],div_b[N],div_c[N];
#define c div_a
#define d div_b
#define e div_c
memset(c,0,sizeof c);memset(d,0,sizeof d);memset(e,0,sizeof e);
for(int i=0;i<n;++i) c[i]=a[i];
for(int i=0;i<m;++i) d[i]=b[i];
std::reverse(c,c+n),std::reverse(d,d+m);
poly_inv(n+m,d,e);
poly_mul((n+m)<<1,e,c);
for(int i=0;i<n-m+1;++i) p[i]=e[i];
std::reverse(p,p+n-m+1);std::reverse(c,c+n),std::reverse(d,d+m);
poly_mul(n,d,p);
for(int i=0;i<m;++i) q[i]=md(c[i]-d[i]+mod);
#undef c
#undef d
#undef e
}
}
