没错,就是这样,用来记录修炼进度。
主要是用来存板子,有些理解了比较长时间的地方会大概讲讲(不过越写越成过程推导了)。
在形式幂级数上定义各种运算,有助于我们更快地解决一些问题,也能更好地为生成函数运算提供快速实现方式。
FFT (Fast Fourier Transform)
一切的开始。
DFT 只需要记住这两个式子就行了:
IDFT 的思想即是单位根反演。记 (G={ m DFT}(F)),即
//luogu3803
#include <bits/stdc++.h>
using namespace std;
const int N=2e6+3e5;
const double pi=acos(-1);
struct cplx
{
double x,y;
cplx(double _x=0,double _y=0) {x=_x,y=_y;}
cplx operator+(const cplx &C) {return cplx(x+C.x,y+C.y);}
cplx operator-(const cplx &C) {return cplx(x-C.x,y-C.y);}
cplx operator*(const cplx &C) {return cplx(x*C.x-y*C.y,x*C.y+y*C.x);}
}f[N],g[N];
int n,m,rv[N];
void FFT(cplx* f,int type)
{
for(int i=0;i<n;++i) if(i<rv[i]) swap(f[i],f[rv[i]]);
for(int p=2;p<=n;p<<=1)
{
int len=p>>1;
cplx w1=cplx(cos(pi/len),type*sin(pi/len));
for(int i=0;i<n;i+=p)
{
cplx wk=cplx(1,0);
for(int k=i;k<i+len;++k)
{
cplx tmp=wk*f[k+len];
f[k+len]=f[k]-tmp,f[k]=f[k]+tmp;
wk=wk*w1;
}
}
}
if(type==-1) for(int i=0;i<n;++i) f[i].x/=n;
}
int main()
{
scanf("%d%d",&n,&m);
for(int i=0;i<=n;++i) scanf("%lf",&f[i].x);
for(int i=0;i<=m;++i) scanf("%lf",&g[i].x);
for(m+=n,n=1;n<=m;n<<=1);
for(int i=1;i<n;++i)
rv[i]=(rv[i>>1]>>1)|((i&1)?n>>1:0);
FFT(f,1); FFT(g,1);
for(int i=0;i<n;++i) f[i]=f[i]*g[i];
FFT(f,-1);
for(int i=0;i<=m;++i) printf("%d ",(int)(f[i].x+0.49));
return 0;
}
拉格朗日插值(Lagrange Interpolation)
众所周知,(n+1) 个不同点可以唯一确定一个 (n) 次多项式 (F(x))。
解多项式有一个比较 naive 的做法是待定系数后高斯消元,但这么做复杂度高达 (O(n^3)) 且还有精度方面的问题。
而拉格朗日插值是一种比较优秀的方法,复杂度只有 (O(n^2))。特别地,如果给出的点横坐标连续,复杂度可以进一步降为 (O(n))。
公式很简单:
可以发现对于每一个 ((x_i,y_i)),和式中仅有 (y_i) 那一项后面的乘积为 1,其他项均为零,满足这个多项式过这 (n+1) 个点的条件。于是如果我们要算 (F(k)) 的话直接带入即可。
有时候可能会有加点操作,如果仍然 (O(n^2)) 重构不太明智,可以考虑变个形
其中
这样每次加点时只需维护 (omega_i) 即可。
#include <bits/stdc++.h>
using namespace std;
const int N=2005,P=998244353;
int n,k,x[N],y[N],ans;
int qpow(int bas,int p)
{
int res=1;
for(;p;p>>=1,bas=1LL*bas*bas%P)
if(p&1) res=1LL*res*bas%P;
return res;
}
int main()
{
n=IO::read(),k=IO::read();
for(int i=0;i<n;++i) x[i]=IO::read(),y[i]=IO::read();
for(int i=0;i<n;++i)
{
int up=1,down=1;
for(int j=0;j<n;++j)
if(i!=j) up=1LL*up*(k-x[j]+P)%P,down=1LL*down*(x[i]-x[j]+P)%P;
ans=(ans+1LL*up*qpow(down,P-2)%P*y[i]%P)%P;
}
printf("%d",ans);
return 0;
}
NTT (Number Theory Transform)
FFT 的缺点很明显:精度问题,三角函数太慢。
我们需要找到一种代替单位根的东西,它必须在形式上符合单位根的运算。
考虑原根这个东西,具体的理论不想展开说了(反正也用不上),只需记住对于形如 (p=k2^t+1) 的素数,设其模 (p) 意义下的原根为 (g),则 (g^0,g^1,dots,g^{p-1}) 互不相同。
容易验证,只需考虑 (omega^1_n=g^{frac{p-1}{n}}),则所有对于单位根的运算对原根也成立。
常见的 NTT 模数:(998244353=119 imes 2^{23}+1,2281701377=17 imes 2^{27}+1,1004535809=479 imes 2^{21}+1),它们的原根都是 (3)。
代码与 FFT 对比基本没变,只是三角函数换成了原根。
所以我为啥比 FFT 慢啊/kk
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int G=3,P=998244353,N=1350000;
ll qpow(ll b,ll p=P-2)
{
ll res=1;
for(;p;p>>=1,b=b*b%P)
if(p&1) res=res*b%P;
return res;
}
const ll invG=qpow(G);
int n,m,tr[N*2];
ll f[N*2],g[N*2],invn;
void NTT(ll f[],bool op)
{
for(int i=0;i<n;++i)
if(i<tr[i]) swap(f[i],f[tr[i]]);
for(int p=2;p<=n;p<<=1)
{
int len=p>>1,tG=qpow(op?G:invG,(P-1)/p);
for(int k=0;k<n;k+=p)
{
ll bfy=1;
for(int l=k;l<k+len;++l)
{
int tt=bfy*f[len+l]%P;
f[len+l]=f[l]-tt,f[l]=f[l]+tt;
if(f[len+l]<0) f[len+l]+=P;
if(f[l]>P) f[l]-=P;
bfy=bfy*tG%P;
}
}
}
}
int main()
{
scanf("%d%d",&n,&m); ++n,++m;
for(int i=0;i<n;++i) scanf("%d",f+i);
for(int i=0;i<m;++i) scanf("%d",g+i);
for(m+=n,n=1;n<m;n<<=1);
for(int i=0;i<n;++i)
tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
NTT(f,1); NTT(g,1);
for(int i=0;i<n;++i) f[i]=f[i]*g[i]%P;
NTT(f,0); invn=qpow(n);
for(int i=0;i<m-1;++i) printf("%lld ",f[i]*invn%P);
return 0;
}
多项式求逆(乘法逆)
求 (F(x)equiv G^{-1}(x)pmod {x^n})。
我们考虑倍增求逆,记 (F=G^{-1}pmod {x^n},F'=G^{-1}pmod {x^{n/2}}),现在我们已经求出了 (F'),考虑推出 (F)。
由主定理,复杂度为 (T(n)=T(n/2)+O(nlog n)=O(nlog n))。
由于我们需要频繁地进行一些操作,适度地封装显得非常有必要。我的多项式全是从 cmd_block 大佬那学来的,然而他的码风实在是……一言难尽。所以以下是根据我的理解封装的,大家可以自行探索合适的封装方式。我抱跳瓜的大腿了,届时会白嫖他的板子,所以代码先咕着。
多项式开方
求 (F(x)^2equiv G(x)pmod {x^n})。
继续考虑倍增,假设已求得 (F'^2equiv Gpmod {x^{n/2}}),于是有
可以看到推法和求逆很像。对于常数项,答案是模意义下的二次剩余,一般来讲做题时常数项都是 (1),如果不是可以骂死出题人。
多项式求导(求积分)
这不用说了吧,这都不会做什么多项式啊。
多项式牛顿迭代
给多项式 (G),求 (F) 使得 (G(F(x))equiv 0pmod {x^n})。
还是考虑倍增,若求得 (G(F'(x))equiv 0pmod {x^{n/2}}),直接泰勒展开:
由于 (F') 最高次项是 (x^{n/2}) 的,所以对于 (ige 2) 的情况模意义为 (0)。
于是有:
(然后你会发现这就是普通的牛顿迭代式子)
有了多项式牛迭,我们可以无脑求解很多式子。
多项式 ln
即求 (F(x)=ln G(x))。
考虑两边同时求导。但要注意两边的主元是不一样的,右边实际上是一个复合函数求导,应用链式法则,有(注意这里多项式带撇表示求导)
四部曲,求逆,求导,乘起来,求积分。
多项式 exp
即求 (F(x)=exp G(x))。
两边先同时取对数,变成 (ln F(x)equiv G(x)pmod{x^n}),移项得 (ln F(x)-G(x)equiv 0pmod {x^n})。如果设 (H(F(x))=ln F(x)-G(x)),我们发现实际要求的其实是这个函数的零点。
所以我们考虑用牛顿迭代来倍增求解。设已求得 (F'(x)equivexp G(x)pmod{x^{n/2}}),则有
注意初值,当 (F) 为常数项时值为 (1)。
多项式快速幂
即求 (F^k(x))。
考虑到 (F^k(x)=exp( kln F(x))),先 (ln) 回来,再 (exp) 回去,完事。复杂度 (O(nlog n))(但常数巨大)。
多项式带余除法
严格的描述是:给定一个 (n) 次多项式 (F(x)) 和一个 (m) 次多项式 (G(x)) ,请求出多项式 (Q(x),R(x)),满足以下条件:
- (Q(x)) 次数为 (n−m),(R(x)) 次数小于 (m);
- (F(x)equiv Q(x)G(x)+R(x)pmod {x^n})。
我们记 (F_r(x)=sum_{i=0}^n[x^{n-i}]F(x)x^i),即 (F_r) 是 (F) 系数翻转后的多项式。很容易得到一个等价的式子:
开始推式子:
式子两边同时 (mod x^{n-m+1}),余数多项式就没了:
(Q(x)) 算了出来,(R(x)) 也就可以很容易算出来了。
分治 FFT
给定序列 (g_{1,dots,n-1}),求序列 (f_{0,dots,n-1})。满足 (f_i=sum_{j=1}^if_{i-j}g_j,f_0=1)。
发现 (f) 的每一项都依赖于前面的项,考虑用 CDQ 分治解决这个过程。
假设当前已经算出了所有 (f_{lsim mid}) 的值,考虑对 (forall mid+1leq ileq r),(f_i) 所需要加的贡献是 (sum_{j=l}^{mid}f_jg_{i-j}),这刚好是一个卷积的形式。所以我们对于当前区间 ([l,r]),只需将 (f_{lsim mid}) 与 (g_{0sim r-l}) 求卷积即可。总复杂度 (O(nlog^2 n))。
考虑在模意义下这道题怎么做。对于 ({f_n}) 与 ({g_n}) 的生成函数 (F(x)=sum_{ige 0}f_ix^i,G(x)=sum_{ige 0}g_ix^i),考虑其卷积
直接多项式求逆即可,复杂度 (O(nlog n))。竟然比分治快