浅谈 FFT （快速傅里叶变换）

直接开始吧。

前置芝士

简单数学（平面直角坐标系，多项式乘法）
码力（好吧其实并没有那么多代码要写）
简单分治（或者，用过线段树也可以）

数论的话，学 NTT 之前再看也可以。

如果你都会就直接开始了。

引入：多项式乘法

七年级课本数学，应该不用怎么讲吧。

例如：

$(x^2+5x+6)(3x^3+x^2-2x-4)$

拆括号：

$3x^5+x^4-2x^3-4x^2+15x^4+5x^3-10x^2-20x+18x^3+6x^2-12x-24$

合并同类项：

$3x^5+16x^4+21x^3-8x^2-6x-24$

上面的方法就相当于：

若 $h(x)=f(x)\times g(x)$ ，则

$h_k=\sum_{i+j=k}f_ig_j$

这样做复杂度是 $O(n^2)$ 的。

暴力代码：（模板题， $55pts$ ）

 #include<bits/stdc++.h>
using namespace std;
int n,m;
int f[1000005],g[1000005],h[1000005];
int main(){
    scanf("%d%d",&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(int i=0;i<=n;i++)
        for(int j=0;j<=m;j++)
            h[i+j]+=f[i]*g[j];
    for(int i=0;i<=n+m;i++)printf("%d ",h[i]);
    return 0;
}

点值转化

考虑优化这个暴力。

我们想，如果我们的乘法只要 $O(n)$ ，那就很妙了。

根据拉格朗日插值，我们知道，在平面直角坐标系中， $n$ 个点就可以确定一个 $n-1$ 次多项式。

用 $n$ 个点来表示一个 $n-1$ 次多项式的方法，就是点值表示。

我们就会想到，能否将一个 $n-1$ 次多项式转化为 $n$ 个点呢？那么两个多项式相乘的点值表示，不就是这 $2n$ 个点值的相乘吗？

例如，我们有两个多项式 $f(x)=x^2+5x+6$ 和 $g(x)=x^2-2x-4$ 。

那么 $f(0)=6,f(1)=12,f(2)=20$ （确定了 $f$ 这个多项式）

$g(0)=-4,g(1)=-5,g(2)=-4$ （确定了 $g$ ）

我们设 $h(x)=f(x)\times g(x)$ 。

那么很显然， $h(0)=-24,h(1)=-60,h(2)=-80$ 。

好像有点不对啊， $h$ 的次数应该是 $4$ ，需要 $5$ 个点啊。

那我们再找两个点不就好了……

好了现在我们知道 $h(0)=-24,h(1)=-60,h(2)=-80,h(3)=-30,h(4)=168$ 。

根据拉格朗日插值公式， $h(x)=\sum_{0\le i\le 4}h(i)\times \prod_{1\le j\le 4,i\ne j}\dfrac{x-j}{i-j}=x^4+3x^3-8x^2-32x-24$

嗯，就是这样。

对于两个次数为 $n$ 的多项式相乘，我们有一种策略：分别代入 $2n+1$ 个数，求出在这两个多项式上的点值，然后把这些点值相乘，最后再把点值转化为多项式。

我们把第一步——求点值叫做 DFT ，第二步——求多项式叫做 IDFT 。综合起来，就叫做 FFT 。

一个显而易见的事实是， IDFT 是 DFT 的CP逆运算。

有关复数

~~不会吧不会吧不会真的有人还没学会复数吧~~

一个复数的形式是 $a+bi$ ，这里 $i$ 是 $-1$ 的平方根之一（另一个是 $-i$ ，虽然他们两个哪个都不大于 $0$ （也不小于 $0$ （因为复数之间好像不能作大小比较来着）））。

复数的运算很简单，就跟多项式差不多（你可以认为 $i$ 是一个未知数）。

$(a+bi)\pm(c+di)=(a\pm c)+(b\pm d)i$

$(a+bi)(c+di)=(ac-bd)+(ad+bc)i$

$\dfrac{a+bi}{c+di}=\dfrac{(a+bi)(c-di)}{(c+di)(c-di)}=\dfrac{(ac+bd)i+(bc-ad)i}{c^2+d^2}$

而这些运算在几何上也有一些~~神奇的~~性质：

假设我们有一个复数 $z=a+bi$ ：

我们让它对应平面直角坐标系上的点 $(a,b)$ 。

定义一个复数的模长为 $|z|=\sqrt{a^2+b^2}$ ，也就是它对应的点到原点 $(0,0)$ 的距离。

定义它的幅角为它对应的点到原点的这条线段与 $x$ 轴的夹角（带正负的角，如果不会自行百度），用式子就是 $\arctan\dfrac{b}{a}$ 。

然后，我们看两个复数 $z_1=a+bi$ 和 $z_2=c+di$ 。

它们的乘积是 $(ac-bd)+(ad+bc)i$ 。

我们把模长代进去： $|z_1\times z_2|=\sqrt{(ac-bd)^2+(ad+bc)^2}=\sqrt{a^2c^2+b^2d^2-2abcd+a^2d^2+b^2c^2+2abcd}=\sqrt{a^2c^2+b^2d^2+a^2d^2+b^2c^2}=\sqrt{(a^2+b^2)(c^2+d^2)}=\sqrt{a^2+b^2}\times\sqrt{c^2+d^2}=|z_1|\times|z_2|$

把幅角的 $\tan$ 值代进去： $\dfrac{ad+bc}{ac-bd}=\dfrac{\dfrac{b}{a}+\dfrac{d}{c}}{1-\dfrac{b}{a}\times\dfrac{d}{c}}$

所以我们得到：复数相乘时，模长相乘，幅角相加。

单位根

学过三角函数的，一定都知道单位圆。

我们把单位圆搬到复平面上，那么在单位圆上的复数，模长一定为 $1$ 。

$Tips:$ 接下来讨论的复数，均是单位圆上的复数，也即模长为 $1$ 的复数。

我们像切蛋糕一样，从 $(1,0)$ 开始，把这个圆切成 $n$ 块，切出的点就是 $\omega_n^k(0\le k<n)$ 。

拿张图过来：

图中标出了 $\omega_3^0,\omega_3^1,\omega_3^2$ 。

显然， $n$ 次单位根是 $x^n=1$ 的复数根。

因此 $\omega_n^{k\pm n}=\omega_n^k$ 。

同时，我们由图容易知道， $\omega_{2n}^{2k}=\omega_n^k$ （或者更一般地，上下标可以同时乘除任意一个正整数~~虽然没用~~）。

DFT 过程

我们把单位根代进多项式 $F(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}$ ！

我们把 $F(x)$ 拆成奇次项和偶次项，提出一个 $x$ ： $F(x)=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2})$ ，设 $FL(x)=(a_0+a_2x+\cdots+a_{n-2}x^\frac{n-2}2),FR(x)=x(a_1+a_3x+\cdots+a_{n-1}x^\frac{n-2}2)$

这里保证 $n$ 是 $2$ 的幂，不会分得不均匀。

我们发现， $FL(x)$ 和 $FR(x)$ 都是关于 $x^2$ 的多项式。

现在我们把单位根代进去。

当 $2k<n$ 时，代入 $\omega_n^k$ ：

$F(\omega_n^k)=FL(\omega_\frac{n}2^k)+\omega_n^kFR(\omega_\frac{n}2^k)$

代入 $\omega_n^{k+\frac{n}2}$ ：

$F(\omega_n^{k+\frac{n}2})=FL(\omega_\frac{n}2^{k+\frac{n}2})+\omega_n^{k+\frac{n}2}FR(\omega_\frac{n}2^{k+\frac{n}2})=FL(\omega_\frac{n}2^k)-\omega_n^kFR(\omega_\frac{n}2^k)$

哈！这两个等式不就只有一个正负号的区别吗？

辣么我们就可以用分治的思路，每次把 $n$ 减半，当然，到 $n=1$ 时，我们代入的只有 $\omega_1^0=1$ ，此时什么也不要做。

IDFT 过程

上面说过， IDFT 是 DFT 的CP逆运算，它可以把点值表示转成系数表示。

其实它们就只有一句话的区别：只要把 DFT 中的 $\omega$ 全部换成 $\omega^{-1}$ ，做完之后除以 $n$ 即可。

证明只需证（使用矩阵求逆的证法）：

显然

我们假设

那么

验证：

所以 $A\times A^{-1}$ 确实是单位矩阵。

递归实现 FFT

 #include<bits/stdc++.h>
using namespace std;
const double pi=acos(-1.0);
struct CP{
    double x,y;
    CP(double xx=0,double yy=0){
        x=xx,y=yy;
    }
    CP operator +(CP const &b)const{
        return CP(x+b.x,y+b.y);
    }
    CP operator -(CP const &b)const{
        return CP(x-b.x,y-b.y);
    }
    CP operator *(CP const &b)const{
        return CP(x*b.x-y*b.y,x*b.y+y*b.x);
    }
}f[2097160],g[2097160],sav[2097160];
int n,m;
void fft(CP *f,int len,bool flag){
    if(len==1)return;
    for(int k=0;k<len;k++)sav[k]=f[k];
    for(int k=0;k<len/2;k++)f[k]=sav[k<<1],f[k+len/2]=sav[k<<1|1];
    fft(f,len/2,flag);
    fft(f+len/2,len/2,flag);
    CP Ur(cos(2*pi/len),sin(2*pi/len));
    if(!flag)Ur.y*=-1;
    CP now(1,0);
    for(int k=0;k<len/2;k++){
        sav[k]=f[k]+now*f[k+len/2];
        sav[k+len/2]=f[k]-now*f[k+len/2];
        now=now*Ur;
    }
    for(int k=0;k<len;k++)f[k]=sav[k];
}
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);
    m+=n;
    n=1;
    while(n<=m)n<<=1;
    fft(f,n,1);fft(g,n,1);
    for(int i=0;i<n;i++)f[i]=f[i]*g[i];
    fft(f,n,0);
    for(int i=0;i<=m;i++)printf("%d ",(int)(f[i].x/n+0.499)); 
    return 0;
}

这样已经足以通过板子了，但众所周知递归非常慢，所以我们要卡卡常（

迭代实现 FFT

显然，递归可以直接改成迭代（手动从最后一层开始做）。可是我们没法知道哪个位置最后对应到哪个位置。

然后，我们发现：每个位置，经过原来的递归后，必定会到另一个位置~~等于没说~~。

关键在于：

$0\to 0$
$1\to 4$
$2\to 2$
$3\to 6$
$4\to 1$
$5\to 5$
$6\to 3$
$7\to 7$

哎！用二进制来表示，最后对应的位置就是原来位置的二进制反过来！

然后我们就可以递推每个数反过来是什么样子，然后就能做了。

 #include<bits/stdc++.h>
using namespace std;
const double pi=acos(-1.0);
struct CP{
    CP (double xx=0,double yy=0){x=xx,y=yy;}
    double x,y;
    CP operator +(CP const &B)const{
        return CP(x+B.x,y+B.y); 
    }
    CP operator -(CP const &B)const{
        return CP(x-B.x,y-B.y); 
    }
    CP operator *(CP const &B)const{
        return CP(x*B.x-y*B.y,x*B.y+y*B.x); 
    }
}f[2097160],g[2097160];
int n,m;
int tr[2097160];
void fft(CP *f,int n,int flag){
    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;
        CP Ur(cos(2*pi/p),sin(2*pi/p));
        if(!flag)Ur.y*=-1;
        for(int k=0;k<n;k+=p){
            CP now(1,0);
            for(int i=k;i<=k+len-1;i++){
                CP num=now*f[i+len];
                f[i+len]=f[i]-num;
                f[i]=f[i]+num;
                now=now*Ur;
            }
        }
    }
}
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);
    m+=n;
    n=1;
    while(n<=m)n<<=1;
    for(int i=0;i<n;i++)tr[i]=(tr[i>>1]>>1)|((i&1)?(n>>1):0);
    fft(f,n,1);fft(g,n,1);
    for(int i=0;i<n;i++)f[i]=f[i]*g[i];
    fft(f,n,0);
    for(int i=0;i<=m;i++)printf("%d ",(int)(f[i].x/n+0.499)); 
    return 0;
}

浅谈 FFT （快速傅里叶变换）

浅谈 FFT （快速傅里叶变换）

前置芝士

引入：多项式乘法

点值转化

有关复数

单位根

DFT 过程

IDFT 过程

递归实现 FFT

迭代实现 FFT

差分约束系统略解

umbrella_leaf

	#include<bits/stdc++.h>
	using namespace std;
	int n,m;
	int f[1000005],g[1000005],h[1000005];
	int main(){
	scanf("%d%d",&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(int i=0;i<=n;i++)
	for(int j=0;j<=m;j++)
	h[i+j]+=f[i]*g[j];
	for(int i=0;i<=n+m;i++)printf("%d ",h[i]);
	return 0;
	}

	#include<bits/stdc++.h>
	using namespace std;
	const double pi=acos(-1.0);
	struct CP{
	double x,y;
	CP(double xx=0,double yy=0){
	x=xx,y=yy;
	}
	CP operator +(CP const &b)const{
	return CP(x+b.x,y+b.y);
	}
	CP operator -(CP const &b)const{
	return CP(x-b.x,y-b.y);
	}
	CP operator *(CP const &b)const{
	return CP(xb.x-yb.y,xb.y+yb.x);
	}
	}f[2097160],g[2097160],sav[2097160];
	int n,m;
	void fft(CP *f,int len,bool flag){
	if(len==1)return;
	for(int k=0;k<len;k++)sav[k]=f[k];
	for(int k=0;k<len/2;k++)f[k]=sav[k<<1],f[k+len/2]=sav[k<<1\|1];
	fft(f,len/2,flag);
	fft(f+len/2,len/2,flag);
	CP Ur(cos(2pi/len),sin(2pi/len));
	if(!flag)Ur.y*=-1;
	CP now(1,0);
	for(int k=0;k<len/2;k++){
	sav[k]=f[k]+now*f[k+len/2];
	sav[k+len/2]=f[k]-now*f[k+len/2];
	now=now*Ur;
	}
	for(int k=0;k<len;k++)f[k]=sav[k];
	}
	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);
	m+=n;
	n=1;
	while(n<=m)n<<=1;
	fft(f,n,1);fft(g,n,1);
	for(int i=0;i<n;i++)f[i]=f[i]*g[i];
	fft(f,n,0);
	for(int i=0;i<=m;i++)printf("%d ",(int)(f[i].x/n+0.499));
	return 0;
	}