你看我都不好意思说是学习笔记了，毕竟(FFT)我怎么可能学得会
那就写一篇抄袭笔记吧ctrl+c真舒服
先从多项式说起吧
1.多项式
我们定义一个多项式
这就是一个(n-1)次的多项式了
比如说(F(x)=x^3+2x^2+x+1)就是一个三次的多项式了
我们还可以把多项式理解成函数，比如说上面那个多项式(F(2)=2^3+2 imes2^2+2+1=19)
很休闲吧，我会的也就这么多了
之后多项式还有两种表达形式
-
系数表示法，就是把系数存下来啊，上面那个多项式就是({1,2,1,1})
-
点值表示法，任何一个(n-1)的多项式都可以被(n)个点完全表示出来，这个好像是拉格朗日插值里的内容了吧，比如上面那个多项式我们可以用((0,1),(1,5),(2,19),(3,49))来表示
这两种表示方法当然是可以互相转化的，系数转点值我们直接暴力找出(n)个点带进去，点值转系数我们直接列出方程来高斯消元就好了
但是一个是(O(n^2))的，一个是(O(n^3))的
还有一个东西是多项式乘法，比如说两个多项式(G(x),H(x))
有
那么
这个东西直接来做也是(O(n^2))的，但是如果给出的是两个点值表示的多项式，那么在(O(n))的时间内就可以做了，就是横坐标相等的点纵坐标直接乘起来
而(FFT)就是为了解决上面的问题的
2.虚数和单位根
再来扯一扯虚数
就是数轴上根本不存在的数了
一个虚数通常长这个样子(a+bi)，其中(i=sqrt{-1})
这个东西长得这么奇怪肯定不在数轴上了，它飘了起来
比如说(a+bi)就对应平面上的((a,b))这个点
是不是很像向量啊，我所以虚数的运算跟向量差不多
我们可以定义结构体来进行虚数的操作
struct complex
{
double r,c;
complex (double r=0,double c=0):r(r),c(c) {};
};
complex operator +(complex a,complex b) {return complex(a.r+b.r,a.c+b.c);}
complex operator -(complex a,complex b) {return complex(a.r-b.r,a.c-b.c);}
complex operator *(complex a,complex b) {return complex(a.r*b.r-a.c*b.c,a.r*b.c+a.c*b.r);}
complex operator /(complex a,complex b) {return complex((a.r*b.r+a.c*b.c)/(b.r*b.r+b.c*b.c),(a.c*b.r-a.r*b.c)/(b.r*b.r+b.c*b.c));}
好像上面最长的虚数除法在(FFT)中并不会用到
除了(a+bi)这种表示方法，虚数还可以通过模长和幅角的表示方法
模长就是(a+bi)在平面上对应的点到原点的距离，就是(l=sqrt{a^2+b^2})
幅角就是和原点连线与(x)轴的夹角，就是( heta =cotfrac{b}{a})
于是(a+bi=cos heta l+sin heta l imes i)
根据一些我背不过的三角函数运算法则，虚数相乘的几何意义就是模长相乘，幅角相加
之后是神奇的单位根
首先对于半径为一的圆，我们把它的圆心定为坐标原点这样的圆叫做单位圆
神奇的单位根来源于这样一个方程
注意这里的(a)可是一个虚数
(1)自然也可以被看成一个虚数，模长为(1)，幅角为(2pi)的虚数
所有说(a)的模长(l)，幅角( heta)肯定会满足如下的性质
所以
对于这样的虚数我们就叫做(n)次单位根，对于上面那一个幅角为(frac{2pi k}{n})的记为(omega^k_n)
由于模长都是(1)，所以单位根就分部在单位圆上，而且还将单位圆(n)等分
单位根有一些非常好的性质这些都是从慎老师那里抄来的
- (omega_{2n}^{2k}=omega_n^k)
这个非常显然，因为模长都是(1)，而(frac{2pi imes 2k}{2n}=frac{2pi k}{n})
- 如果(n)是偶数，(omega_{n}^k=-omega_{n}^{k+n/2})
其实还有更为重要的性质(omega_{n}^k imes omega_{n}^j=omega_{n}^{k+j})
这里的话(omega_n^{k+n/2}=omega_n^k imes omega_{n}^{n/2})
(omega_{n}^{n/2})并不是一个虚数，而是(-1)，所以得证
- ((omega_n^k)^j=omega_{n}^{kj})
好像从上面来看这个东西也不难理解啊
- (omega_n^k=omega_n^{k\%n})
这个不就是多转了一圈吗
之后就是单位根的求解方法，显然(omega_{n}^x=omega_{n}^{x-1+1}=omega_{n}^{x-1} imes omega_{n}^1)
所以我们知道了(omega_n^1)就可以推出来所有的(n)次单位根了
而(omega_n^1)的幅角是(frac{2pi}{n})，所以(omega_n^1=cosfrac{2pi}{n}+sinfrac{2pi}{n}i)
3.DFT
现在才是主角登场了
我们要将两个多项式乘起来显然靠系数表示并不能成功，因为根本没有办法优化
但是点值表示的话却可以在(O(n))时间内完成乘法，所以我们得先试图转化成点值表示法
但是这又是(O(n^2))的了
但是我们带进去的数如果是单位根呢，这样会不会有一些奇妙的性质呢
(DFT)就是在(O(nlogn)) 时间内把系数变成点值的
对于一个多项式(F(x))，我们先将其进行奇偶分类
(G(x))来存奇数项，(H(x))用来存偶数项
使得(F(x)=xG(x^2)+H(x^2))
比如说(F(x)=3x^3+2x^2+x+4)
那么(G(x)=3x+1)，(H(x)=2x+4)
如果我们对于一个(n-1)次多项式代入了(n)次单位根
那么
我们发现((omega_{n}^k)^2=omega_{n}^{2k})的幅角是(frac{2pi imes 2k}{n}=frac{2pi k}{frac{n}{2}})
也就是说((omega_{n}^k)^2=omega_{n/2}^k)
那么
但是如果(k>n/2)的话我们还得取个(\%)
所以
看到反复(/2)了吗，这样算下去复杂度不就是(O(nlogn))了吗
突然发现和慎老师的柿子不一样了
这是慎老师的柿子
好像也非常有道理的样子呢
懒得写了，抄一个慎老师的板子吧
void DFT (complex *f,int len)
{
if(!len) return ;
complex g[len+1],h[len+1];
for (R i=0;i<=len;++i)
g[i]=f[i<<1|1],h[i]=f[i<<1];
DFT(g,len>>1);
DFT(h,len>>1);
complex og1,og;
len<<=1;
og=complex(1,0);
og1=complex(cos(Pi*2/len),sin(Pi*2/len));
for (R k=0;k<len/2;++k)
{
f[k]=og*g[k]+h[k];
f[k+len/2]=h[k]-og*g[k];
og=og1*og;
}
}
4.IDFT
我们已经把系数变成点值了，之后我们就可以把点值大力乘起来了，这样就是两个多项式相乘之后的点值表示了
现在的问题变成了如何快速的把一个点值再转化成系数形式
我们把这个写成矩阵的形式
设上面那三个矩阵分别为(A,B,C)
也就有(A imes B=C)，我们现在已经求出了(C)这个矩阵，需要求出(B)这个矩阵
显然(B=C imes A^{-1})
我们甚至需要对那个矩阵求一个逆
那么恭喜成功优化到了n^3级别
但是(A^{-1})并不需要求逆，因为我们提前知道它是这个样子
先不管为什么，我们如果把这个矩阵和刚才得到的点值做一次矩阵乘法，就能得到系数式了
其实就是把得到的点值当做系数再来做一遍系数转点值就好了
发现这里和(DFT)唯一的不同就是我们的单位根不同了，一个是(omega_n^1)一个是(omega_n^{-1})，其余部分都是一模一样的，所以甚至我们可以把(DFT)和(IDFT)写在一起
void FFT (complex *f,int len,int v)
{
if(!len) return ;
complex g[len+1],h[len+1];
for (R i=0;i<=len;++i)
g[i]=f[i<<1|1],h[i]=f[i<<1];
FFT(g,len>>1,v);
FFT(h,len>>1,v);
complex og1,og;
len<<=1;
og=complex(1,0);
og1=complex(cos(Pi*2/len),v*sin(Pi*2/len));
for (R k=0;k<len/2;++k)
{
f[k]=og*g[k]+h[k];
f[k+len/2]=h[k]-og*g[k];
og=og1*og;
}
}
如果(v=1)，那么就是(DFT)，(v=-1)就是(IDFT)
现在再来看看那个矩阵为什么就是(A^{-1})了
设上面那个矩阵叫(T)，设(A imes T=S)
那么
当(i=j)时
所以
如果(i eq j)
设(x=i-j)
这不一等比数列求和吗，公比为(omega_{n}^x)，首项为(1)
于是
显然(omega_{n}^{nx}=omega_{n}^{nx\% n}=omega_{n}^0=1)
所以
所以我们惊奇的发现(S)是一个单位矩阵，所以(T=A^{-1})证毕
5.蝴蝶变换
我先咕了，先把板子放上
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define maxn 2000005
#define re register
#define eps 1e-6
struct complex
{
double r,c;
complex (double a=0,double b=0) {r=a;c=b;}
}f[maxn],g[maxn],og,og1,t;
complex operator +(complex a,complex b) {return complex(a.r+b.r,a.c+b.c);}
complex operator -(complex a,complex b) {return complex(a.r-b.r,a.c-b.c);}
complex operator *(complex a,complex b) {return complex(a.r*b.r-a.c*b.c,a.r*b.c+a.c*b.r);}
const double Pi=acos(-1);
int n,m,len,rev[maxn];
inline void FFT(complex *f,int v)
{
for(re int i=0;i<len;i++) if(i<rev[i]) std::swap(f[i],f[rev[i]]);
for(re int i=2;i<=len;i<<=1)
{
int ln=i>>1;
og1=complex(cos(Pi/ln),v*sin(Pi/ln));
for(re int l=0;l<len;l+=i)
{
og=complex(1,0);
for(re int x=l;x<l+ln;x++)
{
t=og*f[ln+x];
f[ln+x]=f[x]-t;
f[x]=f[x]+t;
og=og*og1;
}
}
}
}
inline int check(double x) {if(x-eps<0&&x+eps>0) return 1;return 0;}
int main()
{
scanf("%d%d",&n,&m);
for(re int i=0;i<=n;i++) scanf("%lf",&f[i].r),f[i].c=0;
for(re int i=0;i<=m;i++) scanf("%lf",&g[i].r),g[i].c=0;
len=1; while(len<n+m+1) len<<=1;
for(re int i=1;i<=len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);
FFT(f,1),FFT(g,1);
for(re int i=0;i<len;i++) f[i]=f[i]*g[i];
FFT(f,-1);
for (re int i=0;i<=m+n;++i)
{
f[i].r/=len;
if(check(f[i].r)) f[i].r=0;
}
for(re int i=0;i<=m+n;i++) printf("%.0lf ",f[i].r);
return 0;
}
蝴蝶变换就咕咕咕了，就是背一下板子好了
(FFT)尽管功能强大，但是由于做了大量虚数的运算所以经常精度堪忧
这个时候(NTT)就出场了
(NTT)就是可以取模的(FFT)了
(FFT)强大的功能取决于单位根的性质，而(NTT)在模意义下要做到这一点，就需要原根了
首先原根是个啥
定义如下
对于一个质数(p=qn+1)，(n=2^k)，存在(g)使得(g^i)((0<=i<=p-1))是不同的数，就称(g)为(p)的原根
对于最常见的(NTT)模数(998244353=7 imes 17 imes 2^{23}+1)，其原根为(3)
来看看单位根的几条性质原根是否满足
- (omega_{2n}^{2k}=omega_n^k)
我们先定义(n)次原根分别为(1,g^q,g^{2q}...g^{(n-1)q})
显然(omega_{2n}^{2k}=g^{frac{q}{2} imes 2k}=g^{qk})，因为(p=frac{q}{2} imes 2n+1)
非常美妙的是(omega_n^k=g^{qk})，所以这样的性质原根也有
- (omega_n^k imes omega_n^j=omega_{n}^{k+j})
这一条应该这非常显然
(omega_n^k=g^{qk},omega_n^j=g^{qj})，于是(omega_n^k imes omega_n^j=g^{qk} imes g^{qj}=g^{q(k+j)})
非常美妙的是(omega_n^{k+j}=g^{q(k+j)})，所以这样的性质原根也有
- (omega_n^{frac{n}{2}}equiv -1(mod p))
有了这样的性质就可以分治了
因为(omega_n^{frac{n}{2}} imes omega_n^{frac{n}{2}}=omega_{n}^n=omega_{n}^0=1)
所以(omega_n^{frac{n}{2}}equiv 1,-1)，又因为(omega_{n}^0=1)，(omega_n^{frac{n}{2}})不能和它相等，于是只能(omega_n^{frac{n}{2}}equiv -1)
之后(NTT)的板子基本和(FFT)一样真的只是把单位根换成了原根
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#define maxn 3000005
#define re register
#define LL long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read() {
char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
}
const LL g=3;
const LL mod=998244353;
const LL ig=332748118;
LL a[maxn],b[maxn];
int rev[maxn];
int n,m,len=1;
inline LL quick(LL a,LL b) {LL S=1;while(b) {if(b&1ll) S=S*a%mod;b>>=1ll;a=a*a%mod;}return S;}
inline void NTT(LL *f,int o) {
for(re int i=0;i<len;i++) if(i<rev[i]) std::swap(f[i],f[rev[i]]);
for(re int i=2;i<=len;i<<=1) {
int ln=i>>1;LL og1=quick(o==1?g:ig,(mod-1)/i);
for(re int l=0;l<len;l+=i) {
LL og=1,t;
for(re int x=l;x<l+ln;x++) {
t=og*f[ln+x]%mod;
f[ln+x]=(f[x]-t+mod)%mod;
f[x]=(f[x]+t)%mod;
og=(og*og1)%mod;
}
}
}
}
int main()
{
n=read(),m=read();
for(re int i=0;i<=n;i++) a[i]=read();
for(re int i=0;i<=m;i++) b[i]=read();
while(len<=n+m) len<<=1;
for(re int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);
NTT(a,1),NTT(b,1);
for(re int i=0;i<len;i++) a[i]=a[i]*b[i]%mod;
NTT(a,-1);
LL inv=quick(len,mod-2);
for(re int i=0;i<=n+m;i++) printf("%lld ",inv*a[i]%mod);
return 0;
}