多项式各种操作
多项式各种操作
这里的操作还不是特别完备，可能以后会更新？？gugugu
多项式求逆
给定一个多项式\(A(x)\)，求多项式\(F(x)\)满足：
\[F(x)A(x)\equiv 1\pmod{x^n}
\]
当然如果是良心出题人会让你模\(998244353\)之类的模数。
首先要把\(n\)补成\(2\)的次幂，好理解也好写一些。
假设我们现在已经求出了\(F_0(x)\)，满足：
\[F_0(x)A(x)\equiv 1\pmod{ x^{n/2}}
\]
现在要求\(F(x)\)，满足：
\[F(x)A(x)\equiv 1\pmod{x^n}
\]
显然也满足：
\[F(x)A(x)\equiv 1\pmod{x^{n/2}}
\]
那么两式相减可得：
\[F(x)-F_0(x)\equiv0\pmod{x^{n/2}}
\]
两边平方：
\[F^2(x)-2F(x)F_0(x)+F_0^2(x)\equiv0\pmod{x^{n}}
\]
注意这里模数也平方了，因为左边\(<n/2\)的每一项都是\(0\)，平方之后，对于第\(i\)项，\(i\in [n/2,n)\)，可得：\(a_i=\sum_{j=0}^i a_ja_{i-j}\)，注意这里\(j\)和\(i-j\)必然有一个\(<n/2\)，也就是说这一项也为\(0\)，所以是正确的。
那么两边都乘上一个\(A(x)\)：
\[F(x)\equiv2F_0(x)-A(x)F_0^2(x)\pmod{x^n}
\]
那么我们只要递归做就好了。
复杂度\(T(n)=T(n/2)+O(n\log n)\)，所以\(T(n)=O(n\log n)\)。常数可想而知
代码大概长这个样子：
void get_inv(int *r,int *t,int m) {
if(m==1) return t[0]=qpow(r[0],mod-2),void();
get_inv(r,t,(m+1)>>1);
N=1,bit=0;while(N<(m<<1)) N<<=1,bit++;
for(int i=0;i<N;i++) pos[i]=pos[i>>1]>>1|((i&1)<<(bit-1));
for(int i=0;i<m;i++) tmp[i]=r[i];
for(int i=m;i<N;i++) tmp[i]=0;
ntt(tmp,1),ntt(t,1);
for(int i=0;i<N;i++) t[i]=1ll*t[i]*(2%mod-1ll*tmp[i]*t[i]%mod+mod)%mod;
ntt(t,-1);for(int i=m;i<N;i++) t[i]=0;
}
完整代码会贴在下面。
多项式求ln
给定一个多项式\(A(x)\)，求\(\ln A(x)\)，保证\([0]A(x)=1\)。
多项式的对数...大概可以如下定义：
\[\ln(1-F(x))=\sum_{n=1}^\infty -\frac{F^n(x)}{n}
\]
其实就是把这玩意泰勒展开了一下。
那么计算其实很简单了，令：
\[F(x)=\ln A(x)
\]
两边对\(x\)求导：
\[F'(x)=\frac{A'(x)}{A(x)}
\]
然后注意到后面是可以算出来的，直接多项式求逆就好了。
然后两边积分：
\[F(x)=\int \frac{A'(x)}{A(x)} {\rm{d}}x
\]
注意到求导和积分都是可以\(O(n)\)算的，所以总复杂度\(O(n\log n)\)。
牛顿迭代
给定一个函数\(G(x)\)，求多项式\(F(x)\)，使得\(G(F(x))=0\)。
假定我们现在已经求出了\(F_0(x)\)，满足：
\[G(F_0(x))\equiv0\pmod{x^{n/2}}
\]
考虑对\(G(F(x))\)在\(F_0(x)\)的泰勒展开：
\[G(F(x))\equiv\sum_{n=0}^\infty G^{(n)}(F_0(x))(F(x)-F_0(x))^n\pmod{x^n}
\]
这里\(G^{(n)}\)表示\(n\)阶导，注意这里是对\(F(x)\)求导，也就是：
\[\frac{{\rm{d}}^nG(F(x))}{{\rm{d}}F(x)^n}
\]
注意到\(F(x)\)和\(F_0(x)\)的前\(n/2\)项相同，那么可以发现：
\[(F(x)-F_0(x))^2\equiv0\pmod{x^n}
\]
所以后面的项都消掉了，剩下前两项：
\[G(F(x))=G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))=0
\]
然后可以把\(F(x)\)解出来：
\[F(x)=F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}
\]
然后就做完了。
下面举两个栗子看看这玩意有啥用。
多项式求exp
给定一个多项式\(A(x)\)，求\(\exp(A(x))\)，保证\([0]A(x)=0\)。
同样的，用麦克劳林级数定义如下：
\[e^{A(x)}=\sum_{n=0}^\infty \frac{F^n(x)}{n!}
\]
现在是要求\(F(x)\)满足：
\[e^{A(x)}=F(x)
\]
两边取\(\ln\)：
\[A(x)=\ln F(x)\\
\ln F(x)-A(x)=0
\]
构造函数\(G(x)\)如下：
\[G(F(x))=\ln F(x)-A(x)
\]
那么：
\[G'(F(x))=\frac{1}{F(x)}
\]
带入先前的牛顿迭代式子：
\[F(x)=F_0(x)-(\ln F_0(x)-A(x))F_0(x)\\
=F_0(x)(1-\ln F_0(x)+A(x))
\]
带进去算就好了，复杂度\(O(n\log n)\)，这里面还套了个多项式求逆，常数可想而知
多项式开根
给定一个多项式\(A(x)\)，求\(F(x)\)，满足\(F^2(x)\equiv A(x)\pmod{x^n}\)。
同样构造函数\(G(x)\)：
\[G(F(x))=F^2(x)-A(x)=0
\]
那么带进去：
\[\begin{align}
F(x)&=F_0(x)-\frac{F_0^2(x)-A(x)}{2F_0(x)}=\frac{F_0^2(x)+A(x)}{2F_0(x)}\\
&=\frac{1}{2}(F_0(x)+\frac{A(x)}{F_0(x)})
\end{align}
\]
复杂度同样\(O(n\log n)\)。
多项式快速幂
给定一个多项式\(A(x)\)，求\(F(x)\)，满足\(F(x)\equiv A^k(x)\pmod{x^n}\)。
两边取\(\ln\)：
\[\ln F(x)=k\ln A(x)
\]
然后\(\exp\)回来：
\[F(x)=\exp(k\ln A(x))
\]
复杂度\(O(n\log n)\)。
这个也可以实现多项式开根，只是跑的慢一些而已。
模板
多项式\(\exp\)：题目来自洛谷模板：
#include<bits/stdc++.h>
using namespace std;
void read(int &x) {
x=0;int f=1;char ch=getchar();
for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
void print(int x) {
if(x<0) putchar('-'),x=-x;
if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}
const int maxn = 1e6+10;
const int mod = 998244353;
int s[maxn],pos[maxn],N,bit,n,s2[maxn];
int qpow(int a,int x) {
int res=1;
for(;x;x>>=1,a=1ll*a*a%mod) if(x&1) res=1ll*res*a%mod;
return res;
}
void ntt(int *r,int op) {
for(int i=0;i<N;i++) if(pos[i]>i) swap(r[i],r[pos[i]]);
for(int i=1;i<N;i<<=1) {
int wn=qpow(op==1?3:qpow(3,mod-2),(mod-1)/(i<<1));
for(int j=0;j<N;j+=(i<<1))
for(int k=0,w=1;k<i;k++,w=1ll*w*wn%mod) {
int x=r[j+k],y=1ll*w*r[i+j+k]%mod;
r[j+k]=(x+y)%mod,r[i+j+k]=(x-y+mod)%mod;
}
}
if(op==-1) {
int inv=qpow(N,mod-2);
for(int i=0;i<N;i++) r[i]=1ll*r[i]*inv%mod;
}
}
int tmp[maxn],tmp2[maxn],tmp3[maxn];
void get_inv(int *r,int *t,int m) {
if(m==1) return t[0]=qpow(r[0],mod-2),void();
get_inv(r,t,(m+1)>>1);
N=1,bit=0;while(N<(m<<1)) N<<=1,bit++;
for(int i=0;i<N;i++) pos[i]=pos[i>>1]>>1|((i&1)<<(bit-1));
for(int i=0;i<m;i++) tmp[i]=r[i];
for(int i=m;i<N;i++) tmp[i]=0;
ntt(tmp,1),ntt(t,1);
for(int i=0;i<N;i++) t[i]=1ll*t[i]*(2%mod-1ll*tmp[i]*t[i]%mod+mod)%mod;
ntt(t,-1);for(int i=m;i<N;i++) t[i]=0;
}
void integral(int *r,int *t,int m) {
N=1;while(N<(m<<1)) N<<=1;t[0]=0;
for(int i=0;i<m;i++) t[i+1]=1ll*r[i]*qpow(i+1,mod-2)%mod;
for(int i=m+1;i<N;i++) t[i]=0;
}
void derivative(int *r,int *t,int m) {
N=1;while(N<(m<<1)) N<<=1;
for(int i=1;i<m;i++) t[i-1]=1ll*r[i]*i%mod;
for(int i=m-1;i<N;i++) t[i]=0;
}
void get_ln(int *r,int *t,int m) {
derivative(r,tmp3,m);
get_inv(r,tmp2,m);
ntt(tmp2,1),ntt(tmp3,1);
for(int i=0;i<N;i++) tmp2[i]=1ll*tmp2[i]*tmp3[i]%mod;
ntt(tmp2,-1);integral(tmp2,t,m);
for(int i=0;i<=m<<1;i++) tmp2[i]=0;
}
int t2[maxn];
void get_exp(int *r,int *t,int m) {
if(m==1) return t[0]=1,void();
get_exp(r,t,(m+1)>>1);
N=1,bit=0;while(N<(m<<1)) N<<=1,bit++;
for(int i=m;i<N;i++) t[i]=0;
get_ln(t,t2,m);
for(int i=m;i<N;i++) t2[i]=0;
for(int i=0;i<m;i++) t2[i]=(-t2[i]+r[i]+mod)%mod;
t2[0]=(t2[0]+1)%mod;
for(int i=0;i<N;i++) pos[i]=pos[i>>1]>>1|((i&1)<<(bit-1));
ntt(t2,1),ntt(t,1);
for(int i=0;i<N;i++) t[i]=1ll*t[i]*t2[i]%mod;
ntt(t,-1);for(int i=m;i<N;i++) t[i]=0;
}
int main() {
read(n);
for(int i=0;i<n;i++) read(s[i]);
get_exp(s,s2,n);
for(int i=0;i<n;i++) printf("%d ",(s2[i]+mod)%mod);puts("");
return 0;
}
多项式开根：题目来自洛谷模板：
#include<bits/stdc++.h>
using namespace std;
void read(int &x) {
x=0;int f=1;char ch=getchar();
for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
void print(int x) {
if(x<0) putchar('-'),x=-x;
if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}
#define lf double
#define ll long long
const int maxn = 6e5+10;
const int inf = 1e9;
const lf eps = 1e-8;
const int mod = 998244353;
const int inv2 = 499122177;
int f[maxn],g[maxn],n,m,mxn,bit,N,w[maxn],rw[maxn],s[maxn],t[maxn],pos[maxn];
int qpow(int a,int x) {
int res=1;
for(;x;x>>=1,a=1ll*a*a%mod) if(x&1) res=1ll*res*a%mod;
return res;
}
void prepare() {
w[0]=1;w[1]=qpow(3,(mod-1)/mxn);
for(int i=2;i<=mxn;i++) w[i]=1ll*w[i-1]*w[1]%mod;
rw[0]=1,rw[1]=qpow(qpow(3,mod-2),(mod-1)/mxn);
for(int i=2;i<=mxn;i++) rw[i]=1ll*rw[i-1]*rw[1]%mod;
}
void ntt(int *r,int op) {
for(int i=1;i<N;i++) if(pos[i]>i) swap(r[i],r[pos[i]]);
for(int i=1,d=mxn>>1;i<N;i<<=1,d>>=1)
for(int j=0;j<N;j+=i<<1)
for(int k=0;k<i;k++) {
int x=r[j+k],y=1ll*r[i+j+k]*(op==1?w:rw)[k*d]%mod;
r[j+k]=(x+y)%mod,r[i+j+k]=(x-y+mod)%mod;
}
if(op==-1) {
int inv=qpow(N,mod-2);
for(int i=0;i<N;i++) r[i]=1ll*r[i]*inv%mod;
}
}
int tmp1[maxn],tmp2[maxn],tmp3[maxn];
void get_pos(int len) {
for(bit=0,N=1;N<len;N<<=1,bit++);
for(int i=1;i<N;i++) pos[i]=pos[i>>1]>>1|((i&1)<<(bit-1));
}
void poly_inv(int *r,int *b,int len) {
if(len==1) return b[0]=qpow(r[0],mod-2),void();
poly_inv(r,b,len>>1);
for(int i=0;i<len;i++) tmp1[i]=b[i],tmp2[i]=r[i];
get_pos(len<<1);
ntt(tmp1,1),ntt(tmp2,1);
for(int i=0;i<N;i++) b[i]=((2ll*tmp1[i]%mod-1ll*tmp2[i]*tmp1[i]%mod*tmp1[i]%mod)%mod+mod)%mod;
ntt(b,-1);
for(int i=len;i<N;i++) b[i]=0;
for(int i=0;i<len<<1;i++) tmp1[i]=tmp2[i]=0;
}
void poly_sqrt(int *r,int *b,int len) {
if(len==1) return b[0]=r[0],void();
poly_sqrt(r,b,len>>1);
poly_inv(b,tmp3,len);
get_pos(len<<1);
for(int i=0;i<len;i++) tmp2[i]=r[i];
ntt(tmp2,1),ntt(tmp3,1);
for(int i=0;i<N;i++) tmp3[i]=1ll*tmp3[i]*tmp2[i]%mod;
ntt(tmp3,-1);
for(int i=0;i<len;i++) b[i]=1ll*inv2*(b[i]+tmp3[i])%mod;
for(int i=0;i<len<<1;i++) tmp3[i]=tmp2[i]=0;
}
int main() {
read(n);
for(int i=0;i<n;i++) read(g[i]);
for(mxn=1;mxn<=n<<1;mxn<<=1);
prepare();
poly_sqrt(g,t,mxn>>1);
for(int i=0;i<n;i++) printf("%d ",t[i]);puts("");
return 0;
}
代码是不同时期写的，可能相同的函数有点出入。
