快速傅里叶变换

# 前言

这篇文章咕的很久，三角函数似乎没啥用。

# 三角函数

# 前前言

三角函数似乎没啥用。三角函数似乎没啥用。三角函数似乎没啥用。

# 单位圆

一个以原点为中心且半径为 $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;

}