现开个坑，发现自己不会多项式，学习一下。
分治 FFT/NTT
计算
\[f_i=\sum^i_{j=1}f_{i-j}g_j
\]
其中 \(f_0=1\)，\(g\)已给出。
考虑CDQ分治，令分治函数 solve(l,r) 表示计算 \(f_l\) 至 \(f_r\) 的结果。
假设我们已经计算了 \(f_l\) 至 \(f_{mid}\) 的结果，考虑计算它们对 \(f_{mid+1}\) 至 \(f_r\) 的贡献 ，容易发现形式是一个卷积，用 NTT 计算即可。
时间复杂度 \(O(n \log^2 n)\) ，记得清空部分数组。
模板:
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int MX=3e6+7;
const int MOD=998244353,g=3,gi=332748118;
int rev[MX];
ll f[MX],G[MX];
ll ksm(ll a,ll b){
ll ret=1;
while(b){
if(b&1)ret=ret*a%MOD;
a=a*a%MOD;
b>>=1;
}
return ret;
}
struct poly{
ll d[MX];
};
void ntt(poly &z,int len,int inv){
for(int i=0;i<len;i++)
if(i<rev[i])swap(z.d[i],z.d[rev[i]]);
for(int i=1;i<len;i<<=1){
ll gn=ksm(inv?gi:g,(MOD-1)/(i<<1));
for(int j=0;j<len;j+=(i<<1)){
ll g0=1;
for(ll k=0;k<i;k++,g0=g0*gn%MOD){
int x=z.d[j+k],y=g0*z.d[i+j+k]%MOD;
z.d[j+k]=(x+y)%MOD,z.d[i+j+k]=(x-y+MOD)%MOD;
}
}
}
if(inv){
int rev=ksm(len,MOD-2);
for(int i=0;i<len;i++)z.d[i]=z.d[i]*rev%MOD;
}
}
poly A,B;
void calc(int L,int R){
int d=R-L+1,k=ceil(log2(d));
int l=1<<max(k,1),mid=L+R>>1;
for(int i=0;i<l;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(max((int)ceil(log2(d)),1)-1));
ntt(A,l,0);
ntt(B,l,0);
for(int i=0;i<=l;i++)A.d[i]=A.d[i]*B.d[i]%MOD;
ntt(A,l,1);
for(int i=mid+1;i<=R;i++)f[i]=(f[i]+A.d[i-L-1])%MOD;
for(int i=0;i<l;i++)A.d[i]=B.d[i]=0;
}
void solve(int l,int r){
if(l==r)return;
int mid=l+r>>1;
solve(l,mid);
int len=r-l+1;
for(int i=l;i<=mid;i++)A.d[i-l]=f[i];
for(int i=1;i<=r-l;i++)B.d[i-1]=G[i];
calc(l,r);
solve(mid+1,r);
}
int main(){
int n;
cin>>n;
for(int i=1;i<n;i++)cin>>G[i];
f[0]=1;
solve(0,n-1);
for(int i=0;i<n;i++)cout<<f[i]<<' ';
cout<<endl;
}
用处(我会的)：优化dp
例题：#575. 「LibreOJ NOI Round #2」不等关系
题目大意：给定一个由 < 与 > 组成的字符串，长度为 \(n\) ，求有多少个长度为 \(n+1\) 的排列满足填入字符串后不等式成立，答案对 \(998244353\) 取模。
考虑先不计算所有 > 的贡献，容易发现 < 号组成了非常多的连通块。设 \(sz_i\) 表示第 \(i\) 个连通块大小， 总方案数为
\[\frac{n!}{\prod sz_i}
\]
考虑容斥计算 > 的贡献。设 \(f_i\) 表示结束点为 \(i\) 的答案，考虑从之前的联通块转移，记录 \(cnt_i\) 为第 \(i\) 个位置之前有多少个 > ，可以得到转移方程：
\[f_i=\sum_{j=0}^{i-1}f_j\times(-1)^{cnt_i-cnt_j}\times\frac{1}{(i-j)!}
\]
其中 \((-1)^{cnt_i-cnt_j}\) 为加入 \(cnt_i-cnt_j\) 个 > 的贡献的容斥系数。
显然，这个式子是可以用分治 NTT 优化转移的，考虑 \(F(x)=f_x\times(-1)^{cnt_x},G(x)=\frac{1}{x!}\)，则
\[F(x)=-1^{cnt_{i-1}+cnt_i} \times \sum_{j=0}^{i-1}F(j)G(i-j)
\]
用分治 NTT 优化即可，时间复杂度 \(O(n\log^2n)\)。
自己想的板子题：
计算
\[f_n=\sum_{i+j=n,c_n>c_i}f_ig_j
\]
考虑分治 FFT/NTT ，在算内部贡献的时候还要在套一个 cdq ，时间复杂度 \(O(n\log^3n)\)
多项式乘法逆
我怎么还会过这种东西...
给出多项式 \(f(x)\) ，求多项式 \(f^{-1}(x)\)，满足
\[f^{-1}(x)f(x)\equiv 1(\mod x^n)
\]
考虑分治，假设已经求出 \(f(x)\) 在模 \(x^{\frac{n}{2}}\) 下的乘法逆 $f^{-1}_{0}(x) $
则
\[f^{-1}(x) \equiv f_{0}^{-1}(x) (2-f(x)-f_{0}^{-1}(x))
\]
code:
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int MX=3e5+7;
const int MOD=998244353,gx=3,gi=332748118;
ll rev[MX],a[MX],b[MX],c[MX];
ll ksm(ll a,ll b){
ll ret=1;
while(b){
if(b&1)ret=ret*a%MOD;
a=a*a%MOD;
b>>=1;
}
return ret;
}
ll inv(ll x){return ksm(x,MOD-2);}
void ntt(ll *d,int len,int inv){
for(int i=0;i<len;i++)
if(i<rev[i])swap(d[i],d[rev[i]]);
for(int i=1;i<len;i<<=1){
ll gn=ksm(inv?gi:gx,(MOD-1)/(i<<1));
for(int j=0;j<len;j+=(i<<1)){
ll g0=1;
for(ll k=0;k<i;k++,g0=g0*gn%MOD){
int x=d[j+k],y=g0*d[i+j+k]%MOD;
d[j+k]=(x+y)%MOD,d[i+j+k]=(x-y+MOD)%MOD;
}
}
}
if(inv){
int rev=ksm(len,MOD-2);
for(int i=0;i<len;i++)d[i]=d[i]*rev%MOD;
}
}
void inv(int len,int n,ll *f,ll *g){
if(len==1)g[0]=inv(f[0]);
else{
int l=(len+1)>>1;
inv(l,n,f,g);
int k=1<<(int)(ceil(log2(max(len+n,1))));
for(int i=0;i<k;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(max((int)ceil(log2(len+n)),1)-1));
for(int i=0;i<k;i++){
a[i]=f[i];
if(i<l)b[i]=g[i];
else b[i]=0;
}
ntt(a,k,0),ntt(b,k,0);
for(int i=0;i<k;i++)a[i]=b[i]*((2-a[i]*b[i]%MOD+MOD)%MOD)%MOD;
ntt(a,k,1);
for(int i=0;i<len;i++)g[i]=a[i];
}
}
ll f[MX],g[MX];
int main(){
ios::sync_with_stdio(false);
int n;
cin>>n;
for(int i=0;i<n;i++)cin>>f[i];
inv(n,n,f,g);
for(int i=0;i<n;i++)cout<<g[i]<<' ';
cout<<endl;
return 0;
}
多项式除法/取模
给出多项式 \(f(x),g(x)\)，求出 \(Q(x) \times g(x) + R(x) = F(x)\)
直接写式子，设 \(f^R(x)\) 表示将 \(f(x)\) 系数翻转后的多项式，则
\[f^R(x) \equiv Q^R(x)\times g^R(x)(\mod x^{n-m+1})
\]
将 \(Q(x)\) 计算后带入即可计算出 \(R(x)\)
code:
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int MX=3e5+7;
const int MOD=998244353,gx=3,gi=332748118;
ll rev[MX],a[MX],b[MX],tg[MX],tf[MX],vg[MX];
ll ksm(ll a,ll b){
ll ret=1;
while(b){
if(b&1)ret=ret*a%MOD;
a=a*a%MOD;
b>>=1;
}
return ret;
}
ll inv(ll x){return ksm(x,MOD-2);}
void ntt(ll *d,int len,int inv){
for(int i=0;i<len;i++)
if(i<rev[i])swap(d[i],d[rev[i]]);
for(int i=1;i<len;i<<=1){
ll gn=ksm(inv?gi:gx,(MOD-1)/(i<<1));
for(int j=0;j<len;j+=(i<<1)){
ll g0=1;
for(ll k=0;k<i;k++,g0=g0*gn%MOD){
int x=d[j+k],y=g0*d[i+j+k]%MOD;
d[j+k]=(x+y)%MOD,d[i+j+k]=(x-y+MOD)%MOD;
}
}
}
if(inv){
int rev=ksm(len,MOD-2);
for(int i=0;i<len;i++)d[i]=d[i]*rev%MOD;
}
}
void mul(ll *f,ll *g,ll *h,int n,int m){
int k=1<<(int)(ceil(log2(max(n+m,1))));
for(int i=0;i<k;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(max((int)ceil(log2(n+m)),1)-1));
fill(a,a+k,0);
fill(b,b+k,0);
for(int i=0;i<n;i++)a[i]=f[i];
for(int i=0;i<m;i++)b[i]=g[i];
fill(h,h+k,0);
ntt(a,k,0);
ntt(b,k,0);
for(int i=0;i<k;i++)h[i]=a[i]*b[i]%MOD;
ntt(h,k,1);
}
void inv(int len,int n,ll *f,ll *g){
if(len==1)g[0]=inv(f[0]);
else{
int l=(len+1)>>1;
inv(l,n,f,g);
int k=1<<(int)(ceil(log2(max(len+n,1))));
for(int i=0;i<k;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(max((int)ceil(log2(len+n)),1)-1));
for(int i=0;i<k;i++){
a[i]=f[i];
if(i<l)b[i]=g[i];
else b[i]=0;
}
ntt(a,k,0),ntt(b,k,0);
for(int i=0;i<k;i++)a[i]=b[i]*((2-a[i]*b[i]%MOD+MOD)%MOD)%MOD;
ntt(a,k,1);
for(int i=0;i<len;i++)g[i]=a[i];
}
}
void reverse(ll *f,int l,int r){
int cnt=l-1;
for(int i=r;i>=l;i--)b[++cnt]=f[i];
for(int i=l;i<=r;i++)f[i]=b[i];
}
void div(ll *f,ll *g,ll *q,ll *r,int n,int m){
for(int i=0;i<n;i++)tf[i]=f[i];
for(int i=0;i<m;i++)tg[i]=g[i];
reverse(tf,0,n-1);
reverse(tg,0,m-1);
inv(n,n,tg,vg);
mul(tf,vg,q,n,n);
reverse(q,0,n-m);
fill(q+n-m+1,q+n+m,0);
mul(q,g,tf,n-m+1,m);
for(int i=0;i<m;i++)r[i]=((f[i]-tf[i])%MOD+MOD)%MOD;
}
ll f[MX],g[MX],q[MX],r[MX];
int main(){
ios::sync_with_stdio(false);
int n,m;
cin>>n>>m;
n++;
m++;
for(int i=0;i<n;i++)cin>>f[i];
for(int i=0;i<m;i++)cin>>g[i];
div(f,g,q,r,n,m);
for(int i=0;i<n-m+1;i++)cout<<q[i]<<' ';
cout<<endl;
for(int i=0;i<m-1;i++)cout<<r[i]<<' ';
cout<<endl;
return 0;
}
码风太答辩了
多项式 ln
给出多项式 \(F(x)\) ， 求多项式 \(\ln(F(x))\)
设 \(G(x) = \ln(F(x))\) ，两边同时求导，得
\[G'(x)=\ln'(F(x))F'(x)=\frac{F'(x)}{F(x)}
\]
对 \(F(x)\) 求导与逆即可求出 \(G'(x)\)，积分一下就可以算出 \(G(x)\)
code:
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int MX=3e5+7;
const int MOD=998244353,gg=3,gi=332748118;
int rev[MX];
ll ln_t[MX],a[MX],b[MX];
ll f[MX],g[MX];
ll ksm(ll a,ll b){
ll ret=1;
while(b){
if(b&1)ret=ret*a%MOD;
a=a*a%MOD;
b>>=1;
}
return ret;
}
ll inv(ll x){return ksm(x,MOD-2);}
void ntt(ll *d,int len,int inv){
for(int i=0;i<len;i++)
if(i<rev[i])swap(d[i],d[rev[i]]);
for(int i=1;i<len;i<<=1){
ll gn=ksm(inv?gi:gg,(MOD-1)/(i<<1));
for(int j=0;j<len;j+=(i<<1)){
ll g0=1;
for(ll k=0;k<i;k++,g0=g0*gn%MOD){
int x=d[j+k],y=g0*d[i+j+k]%MOD;
d[j+k]=(x+y)%MOD,d[i+j+k]=(x-y+MOD)%MOD;
}
}
}
if(inv){
int rev=ksm(len,MOD-2);
for(int i=0;i<len;i++)d[i]=d[i]*rev%MOD;
}
}
void polyinv(int len,int n,ll *f,ll *g){
if(len==1)g[0]=inv(f[0]);
else{
int l=(len+1)>>1;
polyinv(l,n,f,g);
int k=1<<(int)(ceil(log2(max(len+n,1))));
for(int i=0;i<k;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(max((int)ceil(log2(len+n)),1)-1));
for(int i=0;i<k;i++){
a[i]=f[i];
if(i<l)b[i]=g[i];
else b[i]=0;
}
ntt(a,k,0),ntt(b,k,0);
for(int i=0;i<k;i++)a[i]=b[i]*((2-a[i]*b[i]%MOD+MOD)%MOD)%MOD;
ntt(a,k,1);
for(int i=0;i<len;i++)g[i]=a[i];
}
}
void derivative(ll *f,int n,ll *g){
for(int i=1;i<n;i++)
g[i-1]=f[i]*i%MOD;
g[n-1]=0;
}
void integrate(ll *f,int n,ll *g){
for(int i=1;i<n;i++)
g[i]=f[i-1]*inv(i)%MOD;
g[0]=0;
}
void polyln(ll *f,int n,ll *g){
if(f[0]!=1){
cout<<"Error"<<endl;
return ;
}
int t=n<<1,k=ceil(log2(t));
fill(ln_t,ln_t+t,0);
int l=1<<max(k,1);
derivative(f,n,ln_t);
polyinv(n,n,f,g);
for(int i=0;i<l;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(max((int)ceil(log2(t)),1)-1));
ntt(ln_t,l,0);
ntt(g,l,0);
for(int i=0;i<l;i++)ln_t[i]=ln_t[i]*g[i]%MOD;
ntt(ln_t,l,1);
integrate(ln_t,n,g);
}
int main(){
int n;
cin>>n;
for(int i=0;i<n;i++)cin>>f[i];
polyln(f,n,g);
for(int i=0;i<n;i++)cout<<g[i]<<' ';
cout<<endl;
return 0;
}
多项式exp
给定多项式 \(h(x)\) ，求多项式 \(f(x)=\exp(h(x))\)
使用牛顿迭代法，构造函数 \(G(f(x))=\ln(f(x))- h(x)\) ，显然 \(f(x)\) 为 \(G(x)\) 的一个零点。
运用牛顿迭代法可知
\[f(x) = f_0(x)(1-\ln(f_0(x))+h(x))
\]
code:
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int MX=3e5+7;
const int MOD=998244353,gg=3,gi=332748118;
int rev[MX];
ll ln_t[MX],exp_t[MX],a[MX],b[MX];
ll f[MX],g[MX];
ll ksm(ll a,ll b){
ll ret=1;
while(b){
if(b&1)ret=ret*a%MOD;
a=a*a%MOD;
b>>=1;
}
return ret;
}
ll inv(ll x){return ksm(x,MOD-2);}
void init(int len){
int k=1<<(int)(ceil(log2(max(len,1))));
for(int i=0;i<k;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(max((int)ceil(log2(len)),1)-1));
}
void ntt(ll *d,int len,int inv){
for(int i=0;i<len;i++)
if(i<rev[i])swap(d[i],d[rev[i]]);
for(int i=1;i<len;i<<=1){
ll gn=ksm(inv?gi:gg,(MOD-1)/(i<<1));
for(int j=0;j<len;j+=(i<<1)){
ll g0=1;
for(ll k=0;k<i;k++,g0=g0*gn%MOD){
int x=d[j+k],y=g0*d[i+j+k]%MOD;
d[j+k]=(x+y)%MOD,d[i+j+k]=(x-y+MOD)%MOD;
}
}
}
if(inv){
int rev=ksm(len,MOD-2);
for(int i=0;i<len;i++)d[i]=d[i]*rev%MOD;
}
}
void polyinv(int len,int n,ll *f,ll *g){
if(len==1)g[0]=inv(f[0]);
else{
int l=(len+1)>>1;
polyinv(l,n,f,g);
int k=1<<(int)(ceil(log2(max(len+n,1))));
init(len+n);
for(int i=0;i<k;i++){
a[i]=f[i];
if(i<l)b[i]=g[i];
else b[i]=0;
}
ntt(a,k,0),ntt(b,k,0);
for(int i=0;i<k;i++)a[i]=b[i]*((2-a[i]*b[i]%MOD+MOD)%MOD)%MOD;
ntt(a,k,1);
for(int i=0;i<len;i++)g[i]=a[i];
}
}
void derivative(ll *f,int n,ll *g){
for(int i=1;i<n;i++)
g[i-1]=f[i]*i%MOD;
g[n-1]=0;
}
void integrate(ll *f,int n,ll *g){
for(int i=1;i<n;i++)
g[i]=f[i-1]*inv(i)%MOD;
g[0]=0;
}
void polyln(ll *f,int n,ll *g){
if(f[0]!=1){
cout<<"Error"<<endl;
return ;
}
int t=n<<1,k=ceil(log2(t));
fill(ln_t,ln_t+t,0);
int l=1<<max(k,1);
derivative(f,n,ln_t);
polyinv(n,n,f,g);
init(t);
ntt(ln_t,l,0);
ntt(g,l,0);
for(int i=0;i<l;i++)ln_t[i]=ln_t[i]*g[i]%MOD;
ntt(ln_t,l,1);
integrate(ln_t,n,g);
}
void polyexp(ll *h,int n,ll *f){
if(h[0]!=0)return;
f[0]=1;
n=1<<(int)(ceil(log2(max(n,1))));
//cout<<n<<endl;
for(int t=2;t<=n;t<<=1){
int t2=t<<1,k=1<<(int)(ceil(log2(max(t2,1))));
polyln(f,t,exp_t);
exp_t[0]=(h[0]+1+(MOD-exp_t[0]))%MOD;
for(int i=1;i<t;i++)exp_t[i]=(h[i]+MOD-exp_t[i])%MOD;
fill(exp_t+t,exp_t+t2,0);
init(t2);
ntt(f,k,0);
ntt(exp_t,k,0);
for(int i=0;i<k;i++)f[i]=f[i]*exp_t[i]%MOD;
ntt(f,k,1);
fill(f+t,f+t2,0);
}
}
int main(){
int n;
cin>>n;
for(int i=0;i<n;i++)cin>>f[i];
polyexp(f,n,g);
for(int i=0;i<n;i++)cout<<g[i]<<' ';
cout<<endl;
return 0;
}
例题：
本题有比 \(O(n\log n)\) 跑得更快的 \(O(k \sqrt k)\) 做法，但是单 \(\log\) 的做法还是值得思考的。
注：两题均需生成函数相关知识。
注2：多项式exp有分治NTT的 \(O(n\log^2n)\) 的做法，常数比牛顿迭代法要小非常多
应用：多项式快速幂
给出 \(n\) 次多项式 \(f(x)\) ，求多项式 \(f^k(x) \mod(x^n)\)
当 \(f(x)\) 常数项为1时
\[f^k(x)\equiv\exp(k\ln(f(x)))\mod(x^n)
\]
否则找到最低不为0的项 \(x^i\)
\[f^k(x)\equiv f_i^kx^{ik}\exp(k\ln(\frac{f(x)}{f_ix^i}))\mod(x^n)
\]
todo:
多项式优化递推