蒟蒻数论笔记：讲给数论小白的快速傅里叶变换

题外话

我感觉自己的数论笔记从来没有写完过，从极限到微积分到积型函数到FFT

2021.11.10 为啥过了一个晚上就有15个人看了？？？

正题

1.傅里叶变换：

有关傅里叶变换，##先看视频##

2.前置芝士之泰勒展开

（这里不理解也没有什么问题，知道结论即可）

泰勒展开的本质就是用一个多项式来拟合一个函数在$x_0$处的变化，为此要做到每一次求导都与原函数相同，以此来拟合函数的变化。

$$g(x)=f(x_0)+f^{'}(x_0)(x-x_0)+\frac{f^{''}(x-x_0)^2}{2!}+\frac{f^{'''}(x-x_0)^3}{3!}+……+\frac{f^{n}(x-x_0)^n}{n!}$$

分母是幂函数求导带来的常数项约去

所以，我们可以用它来展开几个函数：$e^{ix},i\sin(x),\cos(x)$其中$i$为$\sqrt{-1}$

已知：

$$e^x=\sum_{0}^{\infty}\frac{(x-x_0)^n}{n!}$$

拓展到复数：

$$e^z=\sum_{n=0}^{\infty}{\frac{z^n}{n!}},z\in C$$

放缩：

$$e^{ix}=1+ix-\frac{x^2}{2!}-\frac{ix^3}{3!}+\frac{x^4}{4!}+\frac{ix^5}{5!}+\cdots$$

已知：

$$\sin(x)=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\cdots$$

乘上一个复数：

$$i\sin(x)=\frac{ix}{1!}-\frac{ix^3}{3!}+\frac{ix^5}{5!}-\frac{ix^7}{7!} +\cdots$$

已知：

$$\cos(x)=\frac{1}{1}-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\cdots$$

得出：

$$e^{ix}=i\sin(x)+\cos(x)$$

也就是欧拉公式。

带入$\pi$得：

$$e^{i\pi}+1=0$$

可以理解为，随着$x$的增大，$e^{ix}$也在一个圆上运动。

于是结合上面的视频可以理解傅里叶变换的公式：
这里变成了向量

$$\int_{-\infty}^{+\infty}g(t)-e^{2\pi i f t} dt$$

加负号是因为变换旋转的方向。

普通的傅里叶是$O(n^2)$的，而快速傅里叶可以做到$O(n \log n)$。

2.前置芝士之复数运算

看到前面的复数运算，有些同学可能已经有点搞不懂了，我一开始也是这样的。

这里就简单讲一下复数的运算，就是数学选修2-2的最后一章，（可以自己看）。

简单来说，复数集是对实数集的一个扩充，在定义上，我们这样来看复数：

$$C=\{a+bi|a , b\in R\}$$

我们把这个集合中的数，形如$a+bi$的叫做复数，$i$是虚数单位，$C$是复数集。

为了简洁，我们直接开始四则运算：

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

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

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

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

3.前置芝士之多项式卷积

先推荐一个回答

博主讲的很好。$QWQ$

像乘法，也可以变成多项式乘法，然后卷积。

$$h(k)=\sum_{i+j=k}f(i)g(j)$$

快乐的$O(n^2)$
在FFT中，我们就是把两个多项式变成了两组向量，然后两组向量卷积求新的向量。

4.前置芝士之点乘法

n阶的函数由n+1个点唯一确定。

传统的多项式是用系数表示法，我们可以换为点值表示法来表示。
两个与$x$有关的多项式：

$$A(x)=(x_0,y_0),(x_1,y_1) \cdots (x_n,y_n)$$

$$B(x)=(x_0,y_0^{'}),(x_1,y_1^{'}) \cdots (x_n,y_n)$$

两个相加就直接加上即可，相乘的话比较麻烦。

因为$C(x)$的项数$2n+1$，所以两个式子都要补上一堆项，然后项多项式加法一样按位运算。

$$A(x)=(x_0,y_0),(x_1,y_1) \cdots (x_{2n},y_{2n})$$

$$B(x)=(x_0,y_0^{'}),(x_1,y_1^{'}) \cdots (x_{2n},y_{2n})$$

乘法：

$$A(x)B(x)=(x_0,y_0 y_0^{'}),(x_1,y_1 y_1^{'}),(x_{2n},y_{2n} y_{2n}^{'})$$

现在我们见证了奇迹的诞生！！$O(n)$的时间复杂度！！！

所以，我们到现在就知道了FFT的流程：

1.得到两个多项式A，B；

2.对于每个多项式，选取n+1个点，得出点值表达$(O(n^2))$

3.点乘法

4.将$C(x)$的点值表达式转化为系数表达式。$O(n^2)$

5.前置芝士之拉格朗日插值

虽然没有用，但是还是要讲一下。

插值主要用在将点值表达式转换为系数表达式上，也就是一个$n$阶的函数由$n+1$个不同的点唯一确定。

常见的插值法是拉格朗日插值法。

由于有n个点，求n-1阶函数。我们先枚举每一个点$x_i$，建立一个函数$f_i(x)$，使得$f_i(x_i)=y_i$,而其他的$x_j,j!=i$带入，$f_i(x_j)=0$。

然后，一整个函数就是由这$n$个函数累加起来的和，在没一个点$x_i$上的值都对应其$y_i$。

然后构造出公式：

$$f(x)=\sum_{i=0}^{n-1} y_i \times \prod_{j!=i} \frac{x-x_j}{x_i-x_j}$$

这就可以用来完成FFT的第三步，但是时间复杂度是$O(n^2)$。FFT继续优化。

读到这里结合前面的芝士，你已经知道了DFT以及IDFT了，接下来就是FFT的诞生之路

FFT优化之路

6.前置芝士之单位复根

n次的单位复根是满足$\omega^{n}=1$的复数$\omega$,因为$\omega^{n}$在复数平面上的模都是一，所以相乘之后还是一，那么所有的$\omega_{i}$都会均匀度分布在单位圆上，$n=8$时；

我们再考虑欧拉公式:

$$e^{ix}=\cos x+i\sin x$$

取$x=2\pi$

$$e^{2\pi i}=1=\omega^{n}$$

从而：

$$\omega=e^{\frac{2\pi i}{n}}$$

把此时的单位根为主次单位根：

$$\omega_n=e^{\frac{2\pi i}{n}}$$

对于其他单位根，

$$\omega^k_n=e^{\frac{2\pi i k}{n}},0\le k <n$$

第一步的优化：快速傅里叶变换之$O(n \log n)$求点值

对一个多项式求点值，我们带入2n个点求点值，并强选$2^t$个单位根，保证n为2的整数次幂，那么把$2^t$带为n。

得：

$$\omega^k_n=e^{\frac{2\pi i k}{n}},0\le k <n$$

有：

$$\omega^{(j+k)\%n}_n=\omega_n^j \times \omega_n^k$$

所以多项式$A(x)$的点值为：

$$\begin{aligned} A(\omega_n^0) &=a_0+a_1+a_2+a_3+\cdots+a_{n-1}\\ A(\omega_n^1) &=a_0+a_1\omega_{n}^1+a_2\omega_n^2+a_3\omega_n^4+\cdots+a_{n-1}\omega_{n}^{n-1}\\ A(\omega_n^2)&=a_0+a_1\omega_n^2+a_2\omega_n^4+a_3\omega_n^6+\cdots+a_{n-1}\omega_n^{2n-2}\\ \cdots\cdots&\\ \end{aligned}$$

考虑把每个点值按照奇偶分开：

$$A(\omega_n^i)=(a_0+a_2\omega_n^{2i}+a_4\omega_n^{4i}+a_6\omega_n^{6i}+\cdots+a_{n-2}\omega_{n}^{i*(n-2)})+(a_1\omega_n^{i}+a_3\omega_n^{3i}+a_5\omega_n^{5i}+\cdots+a_{n-1}\omega_{n-1}^{(n-1)*i})$$

左边提一个$\omega_n^i$

$$A(\omega_n^i)=(a_0+a_2\omega_n^{2i}+a_4\omega_n^{4i}+a_6\omega_n^{6i}+\cdots+a_{n-2}\omega_{n}^{i*(n-2)})+\omega_n^i(a_1+a_3\omega_n^{2i}+a_5\omega_n^{4i}+\cdots+a_{n-1}\omega_{n-1}^{(n-2)*i})$$

然后，把左边变成一个系数为$a_0,a_2,a_4,a_6\cdots$的多项式$FL(x)$,带入$\omega_{n}^{2i}$，也就是$\omega_{n/2}^{i}$（看上面的芝士即可懂了）。

把右边也变成一个系数为$a_1,a_3,a_5,a_7\cdots$的多项式$FR(x)$,带入$\omega_{n/2}^{i}$。
得：

$$A(\omega_{n}^i)=FL(\omega_{n/2}^i)+\omega_{n}^iFR(\omega_{n/2}^i)$$

当$n=8$时：

$$A(\omega_{n}^0)=FL(\omega_{n/2}^0)+\omega_{n}^0FR(\omega_{n/2}^0)$$

$$A(\omega_{n}^1)=FL(\omega_{n/2}^1)+\omega_{n}^1FR(\omega_{n/2}^1)$$

$$A(\omega_{n}^2)=FL(\omega_{n/2}^2)+\omega_{n}^2FR(\omega_{n/2}^2)$$

$$A(\omega_{n}^3)=FL(\omega_{n/2}^3)+\omega_{n}^3FR(\omega_{n/2}^3)$$

$$A(\omega_{n}^4)=FL(\omega_{n/2}^0)+\omega_{n}^4FR(\omega_{n/2}^0)$$

$$A(\omega_{n}^5)=FL(\omega_{n/2}^1)+\omega_{n}^5FR(\omega_{n/2}^1)$$

$$A(\omega_{n}^6)=FL(\omega_{n/2}^2)+\omega_{n}^6FR(\omega_{n/2}^2)$$

$$A(\omega_{n}^7)=FL(\omega_{n/2}^3)+\omega_{n}^7FR(\omega_{n/2}^3)$$

这时我们发现前后的值是一样的，就可以折半计算了：
$T(N)=2T(\frac{N}{2})+O(n)=O(n \log n)$
啊啊啊啊！！冯巨讲的时候好迷茫现在看来真NB!!

第三步的优化：快速傅里叶逆变换之$O(n \log n)$求系数

在前置芝士之拉格朗日插值中，我们在求了$C(x)$的点值表达式后，运用拉格朗日插值法，$O(n^2)$求出系数表达式。

在这里，我们可以用快速傅里叶变换的优秀性质来进行逆变换。
首先，多项式的点值

$$A(\omega_n^{x})=\sum_{i=1}^{n-1}a_i\omega_{n}^{x_i}$$

我们用$b_i$来表示上面的式子进行快速傅里叶变换后的$A(\omega_{n}^i)$。

我们对这些$b_i$又进行一次快速傅里叶变换:

$$C(\omega_n^{x})=\sum_{i=1}^{n-1}b_i\omega_{n}^{x_i}$$

恒等变换：

$$C(\omega_n^{x})=\sum_{i=1}^{n-1}\omega_{n}^{x_i}\sum_{j=1}^{n-1}a_j\omega_{n}{ij}$$

基本操作交换sum:

$$C(\omega_n^{x})=\sum_{j=0}^{n-1}\sum_{i=0}^{n-1}\omega_{n}^{i\times (j+x)}$$

由单位复根的性质得，当$(j+x)=n$时，$j+x$才有值，后面的$\sum$才有意义。
所以：

$$C(\omega_n^{0})=na_0$$

$$C(\omega_n^{1})=na_{n-1}$$

$$C(\omega_n^{2})=na_{n-2}$$

$$C(\omega_n^{3})=na_{n-3}$$

$$\cdots$$

$$C(\omega_n^{n})=na_1$$

所以这就得到了系数了，IFFT万岁！！！

后置甜点之代码

首先是递归版本，常数有点大（板子都过不了）

转载大佬的伪代码：

点击查看代码

int Lim = 1, N, M ;
function FFT(int lenth, complex *A, int flag){
    IF (Lim == 1) return ;
    complex A0[lenth >> 1], A1[lenth >> 1] ;//分成两部分
    for(int j : 0 to lenth by_grow 2) A0[j >> 1] = A[j], A1[j >> 1] = A[j + 1] ;
    FFT(lenth >> 1, A0, flag) ;
    FFT(lenth >> 1, A1, flag) ;
    complex Wn = unit(,) , w = (1, 0) ;
  //Wn是单位根，w用来枚举幂，即我们令主次单位根不变，由于其余单位根都是其整次幂，所以可以用一个w来记录到第几次幂
  /*此处求单位根的时候会用到我们的参数flag……嗯没错就用这一次，并且flag的值域为(-1, 1)……是的，只会有两个值*/
    for(int j : 0 to (lenth >> 1) by_grow 1 with w = w * Wn){
        A[i] = A0[i] + A1[i] * w ;//应用公式，下同 
        A[i + (lenth >> 1)] = A0[i] - A1[i] * w ; //顺便求出另一半，由折半引理可显然。 
    } 
} 
function Main{
    input(N), input(M) ;
    for(i : 0 to N by_grow 1) => input(A) ;
    for(i : 0 to M by_grow 1) => input(B) ; 
    while(Lim < N + M) Lim <<= 1 ;//Lim为结果多项式的长度（暂时），化为2的幂的便于分治（二分）
    FFT(Lim, A, 1) ;//两遍FFT表示从系数化为点值 
    FFT(Lim, B, 1) ;
    for(i : 0 to Lim by_grow 2) => A[i] *= B[i] ;//点乘法，此处需要重定义乘号，因为每一项现在表示的是一个点，有x和y两个属性qwq 
}

cpp:

点击查看代码

void FFT(int Lim,complex *A,int flag){
    if(Lim == 1) return ;
    complex A0[Lim >> 1], A1[Lim >> 1] ;
    for(int i = 0; i <= Lim ; i += 2)
        A0[i >> 1] = A[i], A1[i >> 1] = A[i+1] ;
    FFT(Lim >> 1, A0, flag) ;
    FFT(Lim >> 1, A1, flag) ;
    complex unit = (complex){cos(2.0 * Pi / Lim) , flag * sin(2.0 * Pi / Lim)}, w = complex(1, 0) ;//欧拉公式 
    for(int i = 0;i < (Lim >> 1) ; i ++, w = w * unit) {
        A[i] = A0[i] + w * A1[i] ;
        A[i + (Lim>>1)] = A0[i] - w*A1[i];
    }
}

int main(){
......................
FFT(A, 1), FFT(B, 1) ;
for(i = 0; i <= Lim; i ++) A[i] = A[i] * B[i] ;
FFT(A, -1) ;
......................
}

最后

本蒟蒻的第一篇写完了的数论笔记。
后面来更FFT的优化。
感谢冯巨，
引用:https://www.cnblogs.com/pks-t/p/9251147.html
%%%%%%%%%%%%%%%%%%%%%%%%