快速傅里叶变换
前言
这篇文章咕的很久，三角函数似乎没啥用。
三角函数
前前言
三角函数似乎没啥用。三角函数似乎没啥用。三角函数似乎没啥用。
单位圆
一个以原点为中心且半径为 \(1\) 的圆。
公式
\[\cos(\alpha-\beta)=\cos\alpha\cos\beta+\sin\alpha\sin\beta
\]
\[\cos(\alpha+\beta)=\cos\alpha\cos\beta-\sin\alpha\sin\beta
\]
\[\sin(\alpha+\beta)=\sin\alpha\cos\beta+\cos\alpha\sin\beta
\]
\[\sin(\alpha-\beta)=\sin\alpha\cos\beta-\cos\alpha\sin\beta
\]
作以 \(x\) 正半轴为始边，逆时针旋转 \(\alpha\),\(\beta\),\(\alpha-\beta\) 的终边与单位圆分别交于 \(A(\cos(\alpha),\sin(\alpha))\),\(B(\cos(\beta),\sin(\beta))\),\(C(\cos(\alpha-\beta),\sin(\alpha-\beta))\), \(x\) 正半轴与单位圆交于 \(D(1,0)\),有 \(\lvert AB \rvert = \lvert CD \rvert\)。
\[(\cos\alpha-\cos\beta)^2+(\sin\alpha-\sin\beta)^2=(1-\cos(\alpha-\beta))^2+\sin^2(\alpha-\beta) \]\[\begin{aligned} (\cos\alpha-\cos\beta)^2+(\sin\alpha-\sin\beta)^2&=\cos^2\alpha -2\cos\alpha\cos\beta+\cos^2\beta +\sin^2\alpha-2\sin\alpha\sin\beta+\sin^2\beta\\ &=2-2\cos\alpha\cos\beta-2\sin\alpha\sin\beta\\ (1-\cos(\alpha-\beta))^2+\sin^2(\alpha-\beta)&=1-2\cos(\alpha-\beta)+\cos^2(\alpha-\beta)+\sin^2(\alpha-\beta)\\ &=2-2\cos(\alpha-\beta) \end{aligned}\]\[2-2\cos\alpha\cos\beta-2\sin\alpha\sin\beta=2-2\cos(\alpha-\beta) \]\[\cos\alpha\cos\beta+\sin\alpha\sin\beta=\cos(\alpha-\beta) \]\[\cos(-\beta)=\cos\beta \]\[\sin(-\beta)=-\sin\beta \]\[\cos(\alpha+\beta)=\cos\alpha\cos\beta-\sin\alpha\sin\beta \]\[\begin{aligned} \sin(\alpha+\beta)&=\cos(\dfrac{\pi}{2}-\alpha-\beta)=\cos((\dfrac{\pi}{2}-\alpha)-\beta)\\ &=\cos(\dfrac{\pi}{2}-\alpha)\cos\beta+\sin(\dfrac{\pi}{2}-\alpha)\sin\beta\\ &=\sin\alpha\cos\beta+\cos\alpha\sin\beta \end{aligned}\]\[\begin{aligned} \sin(\alpha-\beta)&=\sin\alpha\cos(-\beta)-\cos\alpha\sin(-\beta)\\ &=\sin\alpha\cos\beta-\cos\alpha\sin\beta \end{aligned}\]
虚数
定义
\[i=\sqrt{-1}
\]
\[a+bi
\]
公式
\[(a+bi)+(c+di)=(a+c)+(b+d)i
\]
\[(a+bi)-(c+di)=(a-c)+(b-d)i
\]
\[(a+bi)\times(c+di)=ac+adi+bci+bd\times(\sqrt{-1})^2=(ac-bd)+(ad+bc)i
\]
单位复根
坐标轴
\(x\) 代表实数 \(y\) 代表虚数。
定义
\(n\) 次单位复根是满足 \(\omega_n^n=1\) 的复数 \(\omega_n\)，一共有 \(n\) 个。
如何理解？画一个单位圆，将其 \(n\) 等分，圆上 \(n\) 个点就是 \(n\) 次单位复根。
因此有：
\[\omega_n=\cos(\dfrac{2\pi}{n})+\sin(\dfrac{2\pi}{n})i
\]
\[\omega_n^k=\cos(\dfrac{2\pi}{n}k)+\sin(\dfrac{2\pi}{n}k)i
\]
消去定理
\[\omega_{nm}^{km}=\omega_n^k
\]
\[\omega_{nm}^{km}=\cos(\dfrac{2\pi}{nm}km)+\sin(\dfrac{2\pi}{nm}km)i=\cos(\dfrac{2\pi}{n}k)+\sin(\dfrac{2\pi}{n}k)i=\omega_n^k
\]
消去定理的推论
\[\omega_n^{n/2}=\omega_2=-1
\]
折半定理
若 'n&1=0' 且 \(n>0\)，有：
\[(\omega_n^{k+n/2})^2=(\omega_n^{k})^2
\]
\[(\omega_n^{k+n/2})^2=(\omega_n^{k}\omega_n^{n/2})^2=(\omega_n{^k}\times(-1))^2=(\omega_n^{k})^2
\]
求和定理
对于任意整数 \(1\le n\) 和不能被 \(n\) 整除的非负整数 \(k\)，有
\[\sum_{j=0}^{n-1}(\omega_{n}^k)^j=0
\]
多项式的表达
点值表达
一个次数界为 \(n\) 的多项式 \(A(x)\) 的点值表达是一个由 \(n\) 个点对组成的集合:
\[\left\{(x_0,y_0),(x_1,y_1),\cdots,(x_{n-2}, y_{n-2}),(x_{n-1},y_{n-1})\right\}
\]
其中 \(x_i\) 互不相同。
一个傻逼的性质， \(A(x) = B(x) + C(x)\) 时有：
\[A(x) = \left\{(x_0,By_0+Cy_0),(x_1,By_1+Cy_1),\cdots,(x_{n-2},B y_{n-2}+C y_{n-2}),(x_{n-1},By_{n-1}+Cy_{n-1})\right\}
\]
\(A(x) = B(x) \times C(x)\) 时有：
\[A(x) = \left\{(x_0,By_0\times Cy_0),(x_1,By_1\times Cy_1),\cdots,(x_{n-2},By_{n-2}\times C y_{n-2}),(x_{n-1},By_{n-1}\times Cy_{n-1})\right\}
\]
系数表达
一个次数界为 \(n\) 的多项式 \(A(x)\) 的系数表达是一个由 \(n\) 个系数组成的多项式:
\[a_0x^0 + a_1x^1 + a_2x^2 + \cdots + a_{n-2}x^{n-2}+a_{n-1}x^{n-1}
\]
系数表达与点值表达的关系
对于任意一个由 \(n\) 个点对组成的点值表达都有唯一一个次数界为 \(n\) 的多项式与其对应。
有 \(y_i = A(x_i)\)。
把其写成矩阵形式：
\( \begin{bmatrix} 1 \quad x_0 \quad x_0^2 \quad &\cdots \quad x_0^{n-1}\\1 \quad x_1 \quad x_1^2 \quad &\cdots \quad x_1^{n-1}\\\vdots\quad\vdots\quad\vdots\quad&\ddots\quad\vdots\\1 \quad x_{n-1} \quad x_{n-1}^2 \quad &\cdots \quad x_{n-1}^{n-1}\\\end{bmatrix}\) \(\begin{bmatrix}a_0\\a_1\\\vdots\\a_{n-1}\end{bmatrix}\) \(=\) \(\begin{bmatrix}y_0\\y_1\\\vdots\\y_{n-1}\end{bmatrix}\)
范德蒙德矩阵在 \(x_i\) 皆不同的情况下有逆矩阵，因此有点值表达，能够确定系数 \(a_i\)。
快速傅里叶变换&快速傅里叶逆变换
规定 \(n\) 是 \(2\) 的 \(q\) 次幂。
已知
\[A(x)=a_0x^0 + a_1x^1 + a_2x^2 + \cdots + a_{n-2}x^{n-2}+a_{n-1}x^{n-1}
\]
求
\[\left\{(\omega_n,y_0),(\omega_n^1,y_1),\cdots,(\omega_n^{n-2}, y_{n-2}),(\omega_n^{n-1},y_{n-1})\right\}
\]
分治！！！
\[A_0(x)=a_0x^0 + a_2x^1 + \cdots + a_{n-2}x^{n/2-1}
\]
\[A_1(x)=a_1x^0 + a_3x^1 + \cdots +a_{n-1}x^{n/2-1}
\]
\[A_2(x)=x
\]
\[A(x)=A_0(x^2)+A_2(x)A_1(x^2)
\]
\[A_0=\left\{(\omega_{n/2},y_0),(\omega_{n/2}^1,y_1),\cdots,(\omega_{n/2}^{n/2-2}, y_{n/2-2}),(\omega_{n/2}^{n/2-1},y_{n/2-1})\right\}
\]
\[A_1=\left\{(\omega_{n/2},y_0),(\omega_{n/2}^1,y_1),\cdots,(\omega_{n/2}^{n/2-2}, y_{n/2-2}),(\omega_{n/2}^{n/2-1},y_{n/2-1})\right\}
\]
\[A_2=\left\{(\omega_{n},y_0),(\omega_{n}^1,y_1),\cdots,(\omega_{n}^{n-2}, y_{n-2}),(\omega_{n}^{n-1},y_{n-1})\right\}
\]
\[A(\omega_{n}^{i})=A_0(\omega_{n}^{2i})+\omega_{n}^{i}A_1(\omega_{n}^{2i})
\]
\[A(\omega_{n}^{i})=A_0(\omega_{n/2}^{i})+\omega_{n}^{i}A_1(\omega_{n/2}^{i})
\]
注意，\(\omega_{n/2}^{i}\) 在 \(A_0\) 和 \(A_1\) 的点值表达中出现了！！！
巨丑的伪代码：
FFT(a):
n=a.length
if n==1
return a
w=e^{2 * pi * i / n}
W=1
A0={a0,a2,...,a(n-2)}
A1={a1,a3,...,a(n-1)}
y0=FFT(A0)
y1=FFT(A1)
for k in range(n/2): // 0 至 n/2-1
y[k]=y0[k]+W*y1[k]
y[k]=y0[k]-W*y1[k]
W*=w
return y
这样写为什么是对的？？？
\[\begin{aligned}
y0_k&=A0(\omega_{n/2}^{k})\\
y1_k&=A1(\omega_{n/2}^{k})\\
y_k&=A(\omega_{n}^{k})\\
A(\omega_{n}^{k})&=y0_k+\omega_{n}^{k}y1_k\\
&=A0(\omega_{n/2}^{k})+\omega_{n}^{k}A1(\omega_{n/2}^{k})\\
A(\omega_{n}^{k+n/2})&=A0((\omega_{n}^{k+n/2})^2)+\omega_{n}^{k+n/2}A1((\omega_{n}^{k+n/2})^2)\\
&=A0(\omega_{n}^{2k})+\omega_{n}^{n/2}\omega_{n}^{k}A1(\omega_{n}^{2k})\\
&=A0(\omega_{n/2}^{k})+(-1)\omega_{n}^{k}A1(\omega_{n/2}^{k})\\
&=A0(\omega_{n/2}^{k})-\omega_{n}^{k}A1(\omega_{n/2}^{k})\\
\end{aligned}
\]
已知
\[\left\{(\omega_n,y_0),(\omega_n^1,y_1),\cdots,(\omega_n^{n-2}, y_{n-2}),(\omega_n^{n-1},y_{n-1})\right\}
\]
求
\[A(x)=a_0x^0 + a_1x^1 + a_2x^2 + \cdots + a_{n-2}x^{n-2}+a_{n-1}x^{n-1}
\]
有
\[
\begin{bmatrix}
1 &1 &1 &1 &1 &\cdots &1\\
1 &\omega_{n}^{1} &\omega_{n}^{2} &\omega_{n}^{3} &\omega_{n}^{4} &\cdots &\omega_{n}^{n-1}\\
1 &\omega_{n}^{2} &\omega_{n}^{4} &\omega_{n}^{6} &\omega_{n}^{8} &\cdots &\omega_{n}^{2(n-1)}\\
1 &\omega_{n}^{3} &\omega_{n}^{6} &\omega_{n}^{9} &\omega_{n}^{12} &\cdots &\omega_{n}^{3(n-1)}\\
1 &\vdots &\vdots &\vdots &\vdots &\ddots &\vdots\\
1 &\omega_{n}^{n-1}&\omega_{n}^{2(n-1)}&\omega_{n}^{3(n-1)}&\omega_{n}^{4(n-1)}&\cdots &\omega_{n}^{(n-1)(n-1)}
\end{bmatrix}
\begin{bmatrix}
a_0\\
a_1\\
a_2\\
a_3\\
\vdots\\
a_{n-1}
\end{bmatrix}
=
\begin{bmatrix}
y_0\\
y_1\\
y_2\\
y_3\\
\vdots\\
y_{n-1}
\end{bmatrix}
\]
\[V_n=
\begin{bmatrix}
1 &1 &1 &1 &1 &\cdots &1\\
1 &\omega_{n}^{1} &\omega_{n}^{2} &\omega_{n}^{3} &\omega_{n}^{4} &\cdots &\omega_{n}^{n-1}\\
1 &\omega_{n}^{2} &\omega_{n}^{4} &\omega_{n}^{6} &\omega_{n}^{8} &\cdots &\omega_{n}^{2(n-1)}\\
1 &\omega_{n}^{3} &\omega_{n}^{6} &\omega_{n}^{9} &\omega_{n}^{12} &\cdots &\omega_{n}^{3(n-1)}\\
1 &\vdots &\vdots &\vdots &\vdots &\ddots &\vdots\\
1 &\omega_{n}^{n-1}&\omega_{n}^{2(n-1)}&\omega_{n}^{3(n-1)}&\omega_{n}^{4(n-1)}&\cdots &\omega_{n}^{(n-1)(n-1)}
\end{bmatrix}\]
\[V_n^{-1}=
\dfrac{1}{n}
\begin{bmatrix}
1 &1 &1 &1 &1 &\cdots &1\\
1 &\omega_{n}^{-1} &\omega_{n}^{-2} &\omega_{n}^{-3} &\omega_{n}^{-4} &\cdots &\omega_{n}^{-(n-1)}\\
1 &\omega_{n}^{-2} &\omega_{n}^{-4} &\omega_{n}^{-6} &\omega_{n}^{-8} &\cdots &\omega_{n}^{-2(n-1)}\\
1 &\omega_{n}^{-3} &\omega_{n}^{-6} &\omega_{n}^{-9} &\omega_{n}^{-12} &\cdots &\omega_{n}^{-3(n-1)}\\
1 &\vdots &\vdots &\vdots &\vdots &\ddots &\vdots\\
1 &\omega_{n}^{-(n-1)}&\omega_{n}^{-2(n-1)}&\omega_{n}^{-3(n-1)}&\omega_{n}^{-4(n-1)}&\cdots &\omega_{n}^{-(n-1)(n-1)}
\end{bmatrix}\]
\[V_n^{-1}V_n=\begin{bmatrix}
1 &0 &0 &\cdots &0\\
0 &1 &0 &\cdots &0\\
0 &0 &1 &\cdots &0\\
0 &\vdots &\vdots &\ddots &\vdots\\
0 &0 &0 &\cdots &1
\end{bmatrix}\]
\[[V_n^{-1}V_n]_{i,j}=\sum_{k=0}^{n-1} (\omega_n^{-ki}/n)(\omega_n^{kj}) =\sum_{k=0}^{n-1}\omega_n^{kj-ki}/n
\]
应用求和引理。
在看原来的柿子：
\[\left\{(\omega_n,y_0),(\omega_n^1,y_1),\cdots,(\omega_n^{n-2}, y_{n-2}),(\omega_n^{n-1},y_{n-1})\right\}
\]
求
\[A(x)=a_0x^0 + a_1x^1 + a_2x^2 + \cdots + a_{n-2}x^{n-2}+a_{n-1}x^{n-1}
\]
有
\[a_j=\dfrac{1}{n}\sum_{k=0}^{n-1}y_k\omega_{n}^{-kj}
\]
稍微变换一下 \(FFT\) 的过程即可。
快速傅里叶变换&快速傅里叶逆变换的优化
改成迭代计算即可，因为在同一层的计算会用到相同的单位复根，可以减少常数。
简单的说，将序列位逆序置换，得到最低层的计算顺序，然后不断向上计算，调整位置（实际上在合并时就自动调整位置了）。
最终代码
#include <bits/stdc++.h>
using namespace std;
const int N = 2e7+11;
const double pi = acos(-1);
struct Comp{
double a, b;
Comp operator + (const Comp &oth)const{return (Comp){a+oth.a, b+oth.b};};
Comp operator - (const Comp &oth)const{return (Comp){a-oth.a, b-oth.b};};
Comp operator * (const Comp &oth)const{return (Comp){a*oth.a-b*oth.b, a*oth.b+b*oth.a};}
}f[N], g[N], w_i[N];
int n, m, rev[N];
void FFT(int len, Comp *a, bool type){
for(int i = 0; i < len; i ++){
rev[i] = (rev[i>>1]>>1) + (i%2?(len>>1):0);
if(rev[i] > i) swap(a[rev[i]], a[i]);
}
for(int d = 1; d < len; d<<=1){
Comp w_n = (Comp){cos(pi/d), sin(pi/d)};
if(type) w_n.b = -w_n.b;
w_i[0] = (Comp){1, 0};
for(int i = 1; i < d; i ++) w_i[i] = w_i[i-1]*w_n;
for(int fir = 0; fir < len; fir += d<<1){
int sec = fir+d;
for(int i = 0; i < d; i ++){
Comp A0 = a[fir+i], A1 = w_i[i] * a[sec+i];
a[fir+i] = A0 + A1, a[sec+i] = A0 - A1;
}
}
}
if(type) for(int i = 0; i < len; i ++)
a[i].a /= len, a[i].b /= len;
return;
}
int main(){
scanf("%d %d", &n, &m);
for(int i = 0; i <= n; i ++) scanf("%lf", &f[i].a);
for(int i = 0; i <= m; i ++) scanf("%lf", &g[i].a);
int len = 1; while(len <= n+m) len <<= 1; cout << len << endl;
FFT(len, f, 0), FFT(len, g, 0);
for(int i = 0; i < len; i ++) f[i] = f[i] * g[i], printf("%.2lf ", f[i].a);
puts("");
FFT(len, f, 1);
for(int i = 0; i <= n+m; i ++)
printf("%d ", (int)(f[i].a+0.5));
return 0;
}
浙公网安备 33010602011771号