声明：因本人水平有限，本文将引用大量链接来进行算法的解释~学算法的同时要学会读论文

这是一种利用C++模板递归实现的FFT，和常见的迭代法FFT相比有可读性好，性能高的优势。

先给出洛谷的提交结果，目前为P1919第一：

目前运行时间为洛谷高精度乘法模板题的第一，内存占用也是C++提交里最少的（目前的第二使用的是AVX2优化的NTT，而我没有手动使用SIMD函数，纯C++编程）。

再贴一下这个FFT的源代码地址：github.com/With-Sky/HintFFT

和使用了这个FFT的高精度整数库的地址：github.com/With-Sky/HyperInt-mini

几种FFT的性能测试：github.com/With-Sky/FFT-Benchmark

给个star呗~

这个高精度库的乘法在计算两个十进制16000000位以内的数相乘（这个范围内使用的是FFT）已经可以基本赶上大名鼎鼎的GMP库了，例如计算一百万个9和一百万个7相乘：

上面是使用 VisualStudio 编译的，实际上这个库是跨平台的（clang，gcc，mingw64均测试过），只要使用 C++14 以上标准进行编译即可~

下面是算法的实现方法。

核心知识是模板函数递归，见上一篇专栏。

这里再用一个打印斐波那契数列第n项的小例子回顾一下：

这是数字信号处理（DSP）中一个很重要的知识点，阿B有很多视频讲解，我就不班门弄斧了，这里推荐一些视频，如果感兴趣请务必看完，想必能感受到数学和算法的美妙之处~

偏理论的三个视频：

3Blue1Brown讲解傅里叶变换的视频

Veritasium真理元素讲解FFT的视频

鹤翔万里的互动视频

偏实际算法的视频：

两个来自油管大佬的视频（上集，下集）

SherlockGn的介绍FFT加速高精度乘法的视频

如果想深入学习，推荐直接读DSP的书，和FFT算法相关论文（需要一定的检索能力）

这里简单概述一下，要实现高精度乘法就需要先实现多项式乘法，而多项式乘法相当于多项式系数序列的卷积。将多项式系数视作离散的信号序列，根据傅里叶变换的卷积定理，信号在时域的卷积相当于在频域直接相乘，在频域相乘是很快的，时间复杂度为 ，n为序列长度。问题就在如何快速地做傅里叶变换和反变换，对于离散傅里叶变换（DFT），我们有与之等价的快速算法，即快速傅里叶变换（FFT），可以实现时间复杂度的DFT。因此高精度乘法总体的时间复杂度也是。

这里给出一个递归的FFT算法的C++代码：

下面是一个简单的迭代FFT算法的C++代码：

这里使用了标准库的复数，有的人对标准库避之不及，认为性能很差，其实复数是一个很简单的库，在开启优化后和手写是没有区别的，除非要使用SIMD加速复数运算，否则不要浪费时间去再实现一遍啦~

上面的FFT代码是最常见的迭代FFT代码，注意到我的命名 radix2 ，代表基2，实际上FFT根据变换的数组长度（或者说分治方法）分为基2，基4，基8等等，它代表数组长度是，，（n为大于等于0的整数），此外还有混合基和分裂基FFT。这些FFT的蝶形图有所不同，一般基4的性能高于基2，基8性能高于基4（但灵活度下降了，因为长度只能是4或8的n次方），基8开始性能提升就不明显了，而且代码也会变得复杂。

DIT与DIF，时间抽取和频率抽取，上面的代码就是DIT的FFT，因为这个算法是将时域的下标按奇数偶数分开，DIF就是将频域下标按奇数偶数分开，具体看DSP相关的书籍，我看的是西安电子科技大学出版社的数字信号处理。这里给出基2频率抽取的代码：

基2是按下标和抽取的，而基4则是，，，抽取，基8同理。

很多讲FFT算法的文章都会忽略对蝶形图的讲解，而着重讲递归的思想，而DSP相关书籍则是由公式直接推出蝶形图，而没有讲FFT的递归，实际上迭代的FFT和其蝶形图所描述的是一种东西。

这里定义一个函数，代表旋转因子：，这是一个复数，实部等于，虚部等于

下面是8点FFT，基2时间抽取的蝶形图（流图）：

频率抽取的蝶形图：

所谓的蝶形单元，就是迭代FFT中每次对两个元素所做的变换：

我们可以把它写成函数，参数是旋转因子，数组指针和两个元素之间的间隔（rank）：

注意时间抽取和频率抽取有不同的蝶形运算单元，是因为它们的蝶形图不同（其实运算顺序是反过来的）

分裂基和基2，基4等抽取算法不同，它是先将下标为偶数的抽取为一组，奇数为一组，同时奇数组当中再按下标分为奇数偶数组，举个例子，对于8点的分裂基，第一组下标为；第二组下标为；第3组下标为。

对这三部分进行FFT再合并即可得到当前长度的FFT。

分裂基对数组长度的要求和基2一样，只要为（n为大于等于0的整数）即可。

分裂基算法在中文互联网上讲解较少，仅推荐一个在线PPT：max.book118.com/html/2017/0803/125780228.shtm

外网则推荐Jörg Arndt的Matters Computational Ideas, Algorithms, Source Code 的第425页，21.5 Split-radix algorithm链接是：link.springer.com/book/10.1007/978-3-642-14764-7

以及ooura的FFT对分裂基的解释：www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmn1_24.html

分裂基适用于序列长度为 2^n 的FFT，并且是运算次数比基8更少，并且更灵活，运算次数的比较见ooura的FFT：www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmn1_26.html#mul3-3

顺便贴上ooura FFT的代码地址，这也是很快的FFT：www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmndl.html

分裂基算法的迭代版相比基2基4等写起来较为复杂，这里仅给出递归版：

这里给出分裂基的蝶形单元：

分裂基的流图：

可以发现分裂基的DIF和基2的DIF一样，要最后对数据进行二进制逆序变换。

设计一个高效的查找表

在上面我们的代码中，旋转因子omega的更新是由每次乘以unit来完成的：

理论上这样做没什么问题，但是忽略了浮点数的运算精度问题，随着omega不断乘以unit，误差会累积到无法接受，具体到加速多项式乘法和高精度乘法就是序列长度过长时得到的结果有误差，或者高精度乘法无法使用过大（如10000）进制来计算，这样运算效率会大大降低。

我们要设计一个类，初始化时传入的参数是旋转因子下标的最大值的对数，并且有一个成员函数，参数是旋转因子下标对数和其指数，为我们返回该旋转因子。比如我们最大计算1024点的FFT，则我们的旋转因子下标最大为1024，是2的10次方，则初始化时传入10,；我们要获得旋转因子下标为512，指数为100的旋转因子，则为该成员函数传入：

可以再写一个expand成员函数，控制查找表扩展到多少范围内的旋转因子，避免初始化而没有用到造成性能的浪费。

由于C++的零成本抽象特性，这样写同时兼顾了可读性和性能，而这个查找表的具体实现则参考github的代码（github.com/With-Sky/HintFFT），我平衡了查找表的查找性能和其对内存的空间占用（内存空间也影响时间，因为更大的表就要初始化更多的元素，在洛谷是吃亏的）。这里设计了两个表，一个是ComplexTable，空间占用较小（初始化快），一个是ComplexTableX，空间占用大（初始化较慢）但查找速度快，大家可以自行选择。

用正变换实现逆变换

简单来说就是序列的共轭做正变换，结果再进行一次共轭，还有归一化（除以序列长度），等价于该序列做逆变换：

DIT与DIF混合使用

时间抽取DIT需要在一开始进行下标的二进制逆序交换，而频率抽取DIF是在最后进行的，根据这一特点，我们进行高精度乘法或者多项式乘法时，只需要先进行一次DIF的正变换，在频域相乘后再做DIT的逆变换即可，而这两次的变换都不需要进行二进制逆序交换。DIF之后的频域是逆序的，但是不会影响频域的相乘，而后逆序又可以直接作为DIT逆变换的输入，最后就是卷积的结果：

只进行一次正变换和逆变换

上面的算法要进行两次正变换和一次逆变换，同时要创建两个复数组，其实有更好的方法。

我们把要进行卷积的两个数组分别放到复数组的虚部和实部，只进行一次正变换，在频域每个元素和自身相乘，再进行逆变换，最后我们的卷积结果就等于复数虚部的二分之一：

如此优化可以减少一次FFT的运行时间，以及减少一个复数组的空间占用。

回到递归算法（关键优化）

常规的的递归算法慢就慢在需要临时数组存储奇偶下标的数，而将输入数组根据下标奇偶分开和存储数组需要额外的时间空间开销，于是我们使用迭代算法，无需额外空间实现原位计算。

普通的迭代算法的问题在哪里呢，观察一下普通迭代算法，对长度为N的数组操作时，最外层的每一次循环都会遍历一次数组，共遍历次（总算法复杂度就是），看起来没什么问题，实际上当数组较长时，每次遍历都会超出CPU缓存，进而会访问更慢的一级缓存甚至内存，而且这个痛苦的遍历要重复次。

解决方法，我们回到递归，这次根据蝶形图来递归：

观察发现，8点的DIT变换依赖于数组前半部分的4点DIT变换和对数组后半部分的4点DIT变换，我们可以先处理前半部分，而前半部分同理依赖于两个两点DIT，这样就是一个递归的结构。这样实现递归，没有常规递归使用额外空间和分割数组的开销，同时，在递归到某一阶段时数据全部都在缓存里。

假设你的缓存有1KB，但是要处理长度为1024的FFT，即大小为16KB的数组（假设复数的实部和虚部都是double类型，那么一个复数占用空间为16字节）。常规迭代算法的10次迭代都是超出缓存的，递归型相当于有6次（）是在缓存内，而只有4次要超出缓存。这样效率就可以提高。同时设定一个阈值，当FFT长度小于该阈值时调用普通的迭代FFT，来减少函数的调用开销。

以基2DIT为例子，先进行前半部分和后半部分的DIT变换，最后将结果合并，算出当前长度的DIT，递归到长度小于等于阈值则调用迭代的即2DIT函数：

模板递归分裂基

一旦我们使用了递归来写FFT，分裂基就可以轻易地被实现了。

我们使用递归实现分裂基的DIT和DIF，递归的思想和上面同理，比如计算8点DIT的FFT要先计算一个4点的DIT和两个2点的DIT，也就是要调用3次递归函数。

递归函数如下：

开头讲的模板函数递归可以派上用场了，模板参数是FFT的长度，用它来代替常规的递归，编译器优化时会将函数展开内联（一般的递归函数不会进行内联）：

手写短FFT

我们手动实现一些长度比较短的FFT，将循环的代码展开，优化乘一些特殊旋转因子的情况以加快速度。

上面的递归函数在递归到这个长度时直接调用短FFT即可，DIT和DIF的短FFT分别写到了长度为32点的FFT，这些代码就比较长了，直接看github上的源码吧（github.com/With-Sky/HintFFT）。

测试所用的源码：github.com/With-Sky/FFT-Benchmark

为了使查表性能不成为瓶颈，这里使用相同的查找表ComplexTableX，并提前初始化好，即初始化表的时间不计入结果。

比较的FFT均不对数组进行二进制逆序，模拟加速多项式乘法的情形，代码如下，fft_dit和fft_dif函数替换成要被测试的FFT函数，统计整个多项式乘法的运行时间：

方便起见，后面的递归FFT均为上面讲的根据蝶形图优化的递归FFT，而非本文一开始给的递归FFT。

被测试的函数依次是迭代法基2FFT，普通函数递归基2FFT，普通函数递归分裂基FFT，模板函数递归分裂基FFT，测试的FFT长度范围为2^10到2^23：

比较FFT1和FFT2，可以发现迭代法在一开始有微弱优势，但随着数据规模增加，性能下降严重，最后不敌递归的FFT2,。

比较FFT2和FFT3可以发现，分裂基和基2相比确实有一定优势，迭代法难以实现的分裂基在递归下实现起来较为轻松，只需要记住递归的规则和分裂基的蝶形单元函数即可。

比较FFT3和FFT4，和普通递归相比，模板函数递归的性能更进一步，但是要多写一些辅助函数。

FFT4和HintFFT已经非常接近了，后者只是多实现了8点、16点和32点的展开的FFT，但是性能提升不明显，还要多写几百行。

FFT4和最初的迭代FFT相比，在数据规模较大时前者运行时间仅为后者的60%，性能提升巨大。

测试平台：

Windows 11 专业工作站版 版本22H2

处理器 AMD Ryzen 7 3700X 8-Core Processor

内存 32.0 GB DDR4 3600MHz

编译器：gcc version 12.2.0 (x86_64-posix-seh-rev2, Built by MinGW-W64 project)

编译参数：-O2 -std=c++14

我的代码还有优化空间，一是没有手工写AVX指令集加速复数运算，如果手工写的话可以更快，二是那几个固定点数展开的FFT函数可以继续写下去，64点，128点...让模板函数更早结束递归，这样速度也许会更快，第三点是提高二进制逆序交换的速度，这个函数的执行时间已经和FFT本身相当了，虽然高精度整数和多项式乘法无需调用这个函数，但是在进行常规的DFT时仍需要调用，所以需要进行优化。

至此，我使用的所有优化方法都分享完了，这些优化方法对于NTT同样适用，期待大家在以上方法的基础上继续优化，刷新记录！

为了保证可移植性和开箱即用的特点（目前甚至可以直接在Platform IO IDE为ESP32编译运行），我没有去实现AVX加速其实是我写AVX也不是很熟练。第二点优化是因为每往后写一个固定点数的函数，代码量就会比前一个固定点数的函数多一倍，就没有那么多精力去实现了。最后的二进制逆序算法的优化则是因为我目前水平有限，不知道如何实现更快的二进制逆序算法，但是，代码是开源的，欢迎大家对代码进行修改，优化！这也是开源的意义之一，用群众的力量，一步步迭代优化！同时也感谢这些开源的大佬，让我有机会学习FFT的源码算法，因此，我也要将代码开源，将知识传递下去！

上面出现的论文资料是我寒假实现高精度算法看的一部分论文，很大一部分资料都是对我用处不大的，查找资料也是很痛苦的过程，因为中文社区的资料确实很少。在两个月的学习中，我从只会基本的迭代FFT学到现在的FFT，心得就是，一个领域，入门可能是很简单的，但是想做出一点成绩，需要深挖前人的成果，需要付出心血。学无止境，与君共勉！