FFT抄袭笔记
你看我都不好意思说是学习笔记了，毕竟\(FFT\)我怎么可能学得会
那就写一篇抄袭笔记吧ctrl+c真舒服
先从多项式说起吧
1.多项式
我们定义一个多项式
这就是一个\(n-1\)次的多项式了
比如说\(F(x)=x^3+2x^2+x+1\)就是一个三次的多项式了
我们还可以把多项式理解成函数，比如说上面那个多项式\(F(2)=2^3+2\times2^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\)轴的夹角，就是\(\theta =cot\frac{b}{a}\)
于是\(a+bi=cos\theta l+sin\theta l\times i\)
根据一些我背不过的三角函数运算法则，虚数相乘的几何意义就是模长相乘，幅角相加
之后是神奇的单位根
首先对于半径为一的圆，我们把它的圆心定为坐标原点这样的圆叫做单位圆
神奇的单位根来源于这样一个方程
注意这里的\(a\)可是一个虚数
\(1\)自然也可以被看成一个虚数，模长为\(1\)，幅角为\(2\pi\)的虚数
所有说\(a\)的模长\(l\)，幅角\(\theta\)肯定会满足如下的性质
所以
对于这样的虚数我们就叫做\(n\)次单位根，对于上面那一个幅角为\(\frac{2\pi k}{n}\)的记为\(\omega^k_n\)
由于模长都是\(1\)，所以单位根就分部在单位圆上，而且还将单位圆\(n\)等分
单位根有一些非常好的性质这些都是从慎老师那里抄来的
- \(\omega_{2n}^{2k}=\omega_n^k\)
这个非常显然，因为模长都是\(1\)，而\(\frac{2\pi\times 2k}{2n}=\frac{2\pi k}{n}\)
- 如果\(n\)是偶数，\(\omega_{n}^k=-\omega_{n}^{k+n/2}\)
其实还有更为重要的性质\(\omega_{n}^k \times \omega_{n}^j=\omega_{n}^{k+j}\)
这里的话\(\omega_n^{k+n/2}=\omega_n^k\times \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}\times \omega_{n}^1\)
所以我们知道了\(\omega_n^1\)就可以推出来所有的\(n\)次单位根了
而\(\omega_n^1\)的幅角是\(\frac{2\pi}{n}\)，所以\(\omega_n^1=cos\frac{2\pi}{n}+sin\frac{2\pi}{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{2\pi\times 2k}{n}=\frac{2\pi 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\times B=C\)，我们现在已经求出了\(C\)这个矩阵，需要求出\(B\)这个矩阵
显然\(B=C\times 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\times T=S\)
那么
当\(i=j\)时
所以
如果\(i\neq 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\times 17\times 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}\times 2k}=g^{qk}\)，因为\(p=\frac{q}{2}\times 2n+1\)
非常美妙的是\(\omega_n^k=g^{qk}\)，所以这样的性质原根也有
- \(\omega_n^k\times \omega_n^j=\omega_{n}^{k+j}\)
这一条应该这非常显然
\(\omega_n^k=g^{qk},\omega_n^j=g^{qj}\)，于是\(\omega_n^k\times \omega_n^j=g^{qk}\times 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}}\times \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;
}
