浅谈FFT--从DFT到*CZT,及一些技巧

~~本文为日报特供版QAQ~~ 考虑到前面已经有人写过FFT了，所以部分FFT的基础内容**不会**出现在本篇，如基本的公式，NTT（它仅涉及单位复根循环的性质）等。 本文主要是讲解一些竞赛中会用到的围绕着FFT算法的一些内容，如DFT过程的几何直观、卷积和\*CZT等，另外，最后也会介绍一些冷门却好用的tricks。 之前日报的FFT可以看[这里](https://www.luogu.com.cn/blog/command-block/fft-xue-xi-bi-ji) 如果是FFT入门建议看下我之前的[文章](https://www.luogu.com.cn/blog/105254/ru-he-kuai-su-li-xie-fftifft)（有NTT）（~~**排版警告QWQ**~~）（另外这篇文章我没太大精力校对，所以其中的小漏洞可能较多。但**对于FFT实现的部分**还是很准确的） --- ## 目录 1.DFT简介 2.DFT/IDFT的几何直观 3.卷积初步 4.\*CZT 5.卷积与多项式 6.\*卷积的意义浅谈 7.一些tricks --- ## 1. DFT简介 FT，即 **傅里叶变换**，在方式上，它实质就是·是用三角函数拟近函数。傅里叶变换在不同领域有着不同的意义，但它通常用来**便捷地计算卷积**。 不过连续的傅里叶变换显然不是计算机能够处理的。在此基础上，就有了 **离散傅里叶变换** (此处指变换前后均离散)，即 **DFT**。 实际上，DFT 就是 Z变换 的一种，或者说 **Z变换 其实就是 DFT 的拓展**。具体关于 Z变换的内容会在 第四节 中讲到。 . $F[k]=\sum\limits_{n=0}^{N-1}x[n]e^{-i\frac {2\pi} N nk}$ $x[n]=\frac 1 N \sum\limits_{k=0}^{N-1}F[k]e^{i\frac {2\pi} N nk}$ 这即是DFT的式子。 **注意其中$-i$与$i$其实是不分前后的**。这和其几何直观息息相关，我们在下节马上就会讲到。 ## 2. DFT/IDFT的几何直观 本节主要是尝试建立DFT/IDFT的一种几何直观。 . 首先，自不用说是单位根。设单位复根$\omega$满足$\omega^N=1$，且这里$N$是指满足式子的最小正整数。  由复数的运算法则（幅角相加模相乘），我们其实可以将 与$\omega$的积 视为 **绕原点旋转$\frac {2\pi} N$** ； 同时，我们可以发现由不同次数的单位根组成的点是**在单位圆上平均分布的。** 注意，这里的**旋转**和**平均**的思想在下文十分重要。 . 接下来我们分析式子的几何意义。 观察式子（已用单位复根替换）： $F[k]=\sum\limits_{n=0}^{N-1}x[n]\omega^{nk}$ $x[n]=\frac 1 N\sum\limits_{k=0}^{N-1}F[k]\omega^{-nk}$ DFT时，对于所有的$x[n]$，变换后的每个分别点**绕原点旋转了$N$次，每次都是第一个角度的倍数**。这即是DFT的过程。 我们可以画出它们：  可以发现： 1.部分点所在的半径重叠 2.对于$x[n]$项，最长路径为$n$圈，每个点旋转的角度是$w^n$的倍数，且它们是在路径上**平均分布**的。 [1]即是我们FFT优化的根本(但这和本节没关系)； [2]是一个很有趣的性质，它是由遍历单位复根次数导来的，**它和IDFT的几何直观息息相关**。 我们先设设$P[n]=\sum\limits_{k=0}^{N-1}x[n]\omega^{nk}$ 发现，若求$\sum\limits_{k=0}^{N-1}P[k]$，**其和为0！**（想像一下，每个“向量”都有一个相反的“向量”与它相互抵消） 但是，还有个唯一的特殊情况：若$\omega$指数的乘积式子中有一个数为$0$，这个点就“**不会旋转**”：它**旋转的第一个角度为$0$，增倍后仍为$0$**。 这样，$\sum\limits_{k=0}^{N-1}P[k]$的得数就**为$N$倍的**未旋转的向量了。同时，这个向量也**正好在实数轴上。** 这其实就和做IDFT主过程后不需要做转换，并要乘$\frac 1 N$有关。 . 然后我们再思考下$\omega^{-nk}$在式子中的作用。 为了便于区分，我们重新定义$\omega^{-nk}$为$\omega^{-n'k'}$。 对于每项$x[n]\omega^{nk}$，其IDFT后的值为$x[n]\omega^{nk-n'k'}$。 这个式子看起来没什么可化简的，但我们**注意到$k$和$k'$的值实际上是一样的**（可以观察下求和中$k$的变化）。 于是就有$x[n]\omega^{k(n-n')}$。这意味着对于$P[n]$ IDFT后的图像，只不过是旋转的角度等倍缩小了，每个向量的分布仍是**均匀**的； 再观察DFT的式子中$k$的变化范围，发现对于$x[n]$项旋转的的圈数为$n$，而$(n-n')$后是整数，所以圈数仍为**整数**。 因此，我们上面推出的**求和为$0$的结论**也是成立的。 . 那$\omega^{-nk}$有什么具体意义呢？ 发现特别的当$n=n'$时，对于$x[n]$所有项旋转的角度**均为零**！所以$\omega^{-nk}$的实质就是**具体控制“特殊情况”出现在哪一项。** 这样，我们就构建起了DFT/IDFT的一种几何直观（~~但变换公式本身是如何推得的还要另外讨论~~） ## 3. 卷积初步 卷积的应用十分广泛，且在不同的领域有不同的意义。本节主要是单纯地介绍卷积和其代码实现，对其意义不做探讨。 . 卷积可以分为两种：**线性卷积** 和 **循环卷积** 。下面我会分别讲解它们： （科普：$[P]$指当表达式$P$成立时其值为$1$，否则为$0$） （约定：仅带条件的$\sum$中的条件变量通常为整数） - ### I. 线性卷积 其式子为$c_k=\sum\limits_{0\leq i<n}\sum\limits_{0\leq j<m}[0\leq k-i]A[k-i]B[j](0\leq k < n+m-1)$ 我们也可以从这个式子知道线性卷积结果的规模为$n+m-1$ 可以做个图，更生动些：  （图中限于篇幅只表示了一项） 仔细观察的话可以发现，这形式就像我们平时所做的 **多项式乘法** 。 - ### II. 循环卷积 其式子为$C[k]=\sum\limits_{0\leq i<n}\sum\limits_{0\leq j<n}A[(n+k-i)\mod n]B[j](0\leq k < n)$ 或$C[k]=\sum\limits_{0\leq i<n}\sum\limits_{0\leq j<n}[(i+j)\mod n=k]A[i]B[j](0\leq k<n)$ 注意到，特别的是，循环卷积要求两个序列**长度相等**，且结果的规模是**不变**的。 我们可以将取模想象为在序列后面加一个**负数编号的“假序列”**，做一张图：  或者，也可以形象地理解为先将一个序列反序，再将两个序列**分别**按 逆时针 和 顺时针 顺序放到两个同心圆上，逐位旋转相乘，每种旋转状态下**的积的和**就是结果序列该项的值（细节实现可能还要一些规定）。 就像这样：  - ### III. DFT 与 循环卷积，及一些序列操作 扯了这么多，难道我们只能朴素地$O(n^2)$计算卷积吗？ 并不，因为 DFT 和 卷积其实是大有关系的。 #### III.I 谈相加 首先扯一点简单的结论：DFT后对序列做逐项相加后IDFT，是**等同**于 对原序列做逐项相加的。 我们可以想到，幅角相同的复数相加后，其**幅角是不变的**； 也就是说，DFT后对序列做逐项相加，其**贡献**和 对原序列做逐项相加是一样的：每一位的每一个变换映射都对应相加了（如果您真的看懂了第二节，这个应该是很显然w）。由此得证。 #### III.II 谈相乘 接下来是一个关键的结论：**DFT后对序列做逐项相乘后IDFT，是等同于 对原序列做循环卷积的。** 这是一个非常重要的性质，也正是这个性质使得我们可以用 FFT 优化序列的操作。 . 但我尝试了很多办法，都不能找到一个简洁的图像表示法证明这个结论（~~窝太太菜惹QAQ~~），这里没办法，只能给出一个利用**单位根**（这里定义循环节终点一定是$1$）**反演**（$[n|k]=\frac 1 n \sum\limits_{0\leq d<n}\omega_n^{dk}$）的证明： $C[k]=\sum\limits_{0\leq i<n}\sum\limits_{0\leq j<n}[(i+j)\mod n=k]A[i]B[j]$ $C[k]=\sum\limits_{0\leq i<n}\sum\limits_{0\leq j<n}A[i]B[j][n|(i+j-k)]$ $C[k]=\sum\limits_{0\leq i<n}\sum\limits_{0\leq j<n}A[i]B[j](\frac 1 n\sum\limits_{0\leq d<n}\omega_n^{d(i+j-k)})$ $C[k]=\frac 1 n\sum\limits_{0\leq d<n}\omega_n^{-dk}(\sum\limits_{0\leq i<n}A[i]\omega_n^{di})(\sum\limits_{0\leq j<n}b_j\omega_n^{dj})$ - ### IV. 循环卷积和线性卷积的关系 经过尝试，我们可以发现：**如果将两个序列补零拓展到$n+m-1$项，再对它们做循环卷积，其结果和对它们做线性卷积是一样的** 还是上面那张图，我们对$0$相关的元素高光处理下就知道为什么：  可以发现，除去和$0$相乘的“无效项”后，其形式其实就是个**线性卷积**！ 而通过式子我们也可以发现，对于$k-i<0$时，两项总有一项是补上去的零（即下标超出原序列范围，可显然得证），**使这组的贡献无效。** 这样，我们便可以用 循环卷积 来计算 线性卷积，从而弥补 DFT只能快速计算循环卷积 的缺陷。 - ### V. 一种常见卷积形式及代码实现 在题目中，我们常常会遇到这种卷积形式： $P[k]=\sum\limits_{i=0}^{n-1}f(k-i)g(i)$ 我们可以对$f(x)$在`[-n+1, n-1]`的整数取值$x$从小到大对应序列编号从小到大做一个长$2n-1$的序列，并同样地对$g(x)$取值`[0, n-1]`地做一个序列。 （即$f(x)$对应$f[n-1+x]$，$g(x)$对应$g[x]$） 上式即可化成： $P[k]=\sum\limits_{i=0}^{n-1}f[n-1+k-i]g[i]$ 我们可以发现这很像一个完整的**线性卷积**形式。实际上，它就是在计算$f[i]$与$g[i]$的线性卷积的`[n-1, 2n-2]`项。 我们可以形象地做张图：  另外，如果$f(x)$的代入值是$(i-k)$的话，我们可以 **在函数定义上更改(紧贴着变量加负号) 或 将序列反转** ，以化成上述形式。 #### IV.I. \*小优化 其实这里的卷积长度可以只做到`[0, 2n-2]`，或者说不做末$n-1$项。我们可以保证答案不受影响。 具体的证明可以详见 ***6. IV*** 。 #### IV.II.代码示例 这里给出关键部分的代码（且已经使用了 *3.IV.I* 所述的优化）： ```cpp //常见卷积形式代码 //n 如上文所述，N 是指第一个大于等于 2n-1 的2的幂次数（即 FFT的规模） //数组缺省值为 0 inline void Conv(){ for(int i =-n+1; i <= n-1; ++i) f[n-1+i].real =F(i)/*计算 f的函数*/; for(int i =0; i < n; ++i) g[i].real =G(i)/*计算 g的函数*/; pre()/*预处理*/, FFT(f), FFT(g); for(int i =0; i < N; ++i) Q[i] =f[i]*g[i]; IFFT(Q); for(int i =0; i < n; ++i) ans[i] =int(Q[i+n-1].real/N+0.5); } ``` ## 4. \*CZT 我们知道DFT是为了快速计算卷积，但为什么要计算卷积呢？ 原因之一是卷积其实描述了一种非常常见的情景--"信号的叠加"。DFT，就可以视为快速计算的桥梁。 但显然真实的情景肯定不止这一种，于是就有了定义更广泛的 **Z变换**。 其式子的一种形式为： $Z[k]=\sum\limits_{n=0}^{N-1}x[n]z[k]^{-n}$ 注意这里的$z[k]$是规定的变换对序列。而变换对的规定，就和**具体需要变换的序列**有关。 如做DFT实际上就可以规定$z[k]=w^k$。 当然Z变换的变换对还有很多，每个变换对也有不同的含义；同时其通用的逆变换形式也较复杂。这些和下文没有太大联系，这里就不展开了。 - ### I. 计算Z变换 显然朴素计算Z变换是$O(n^2)$的，我们需要更优的算法。可Z变换**既不一定满足（类似的）蝴蝶定理，变换的规模$N$也不一定可以表达成$2$的幂次**，这要怎么办？ 其实若对于一般的Z变换，确实是很难找到快速计算的算法的。不过如果限定一些条件，我们还是有比较显然的方法的。 - ### II. CZT 现在我们限定$z[k]$是**对复平面上一段螺线的等分角抽样。** 即$z[k]=AW^{-k}(0<k<M)$。其中$A,W$是任意复数，$M$和$N$无关。 这样就可以尝试把变换化成卷积的形式了： $Z[k]=\sum\limits_{n=0}^{N-1}x[n]z[k]^{-n}$ $Z[k]=\sum\limits_{n=0}^{N-1}x[n]A^{-n}W^{nk}$ 接下来我们要想办法分别化出**仅两个（在求和式子中）自变量与$k-n$、与$n$成线性关系的函数**。 最简单地，可以想到用**平方式**： $Z[k]=\sum\limits_{n=0}^{N-1}x[n]A^{-n}W^\frac {(k-n)^2-k^2-n^2} 2$ $Z[k]=W^\frac {-k^2} 2\sum\limits_{n=0}^{N-1}A^{-n}W^\frac {-n^2} 2x[n]\cdot W^\frac {(k-n)^2} 2$ 显然，只要设$f(k-n)=W^\frac {(k-n)^2} 2,g(n)=A^{-n}W^\frac {-n^2} 2x[n]$，这就是我们上节讲的“**常用卷积形式**”了，直接套模板就行了。 代码就不放惹。 这即是 **CZT算法**。 . 当然 CZT 本身的含义、其几何直观也是可以展开来讲的。但和 Z变换 一样，其坑太深并且和 DFT、卷积 没太大关系，这里就不展开来讲了。 读者如果有兴趣的话，可以参考 [这篇ppt](https://wenku.baidu.com/view/62110fae0c22590102029d7a.html) 对CZT的讲解。 ## 5. 卷积与多项式 我们发现多项式相乘，其系数的变化规律其实就是**做线性卷积**。 于是直接带入模板求解即可。 ## 6. \*卷积的意义浅谈 上文我们曾提到：“卷积其实描述了一种非常常见的情景--"信号的叠加"”。但具体又是怎样的呢？下面，我会给出一个简单的例子： . 我们先定义信号强度函数$h(T)(0\leq T<N)$，其中$T$指**离信号发生的时刻的距离**，负数指发生前，正数亦然；$N$指总的统计范围。 然后再定义一个信号密度函数$g(t)$，指$t$时刻**发生的信号的数量**。 定义单位时刻信号总强度为**该时刻各个信号强度和**，题中单位均给定统一。 . 现在，我们要求求所有可能的时刻$t$的信号总强度。 朴素的模拟每个信号是$O(n^2)$的，但其实我们可以从另一个角度入手。 我们设$P(k)$为$k$时刻信号总强度。若我们考虑**每个信号到达$k$时刻**时的强度，便可以列出式子： $P(k)=\sum\limits_{t=0}^{N-1}h(k-t)g(t)$ 这非常显然就是个卷积的式子。 ## 7. 一些tricks - ### I. 打包FFT 我们想到，既然变换的本质是对向量做**旋转**，我们可不可以构造斜率不同的过原点的直线作为数轴，将实数一一**映射**到它们上面（即调整斜率），从而一次打包多个多项式的FFT？ 构造任意数轴并映射所需要的计算量是巨大的；但不同的是正交的两个数轴**虚数轴**和**实数轴**有特别的性质。 因此，接下来我们考虑打包两个序列做FFT，即将其中一个序列**映射到虚数轴**上。 #### I.I 对于IDFT的打包 我们只需将两个DFT后的序列（且DFT时的旋转方向一致）的其中一个沿DFT的旋转方向**再旋转90度**（即乘$i$或$-i$）后叠加，做完IDFT后直接检索 虚数轴 和 实数轴 即可（做旋转处理的序列会在**虚数轴**上）。 #### I.II 对于DFT的打包 DFT的打包比起IDFT的打包要复杂得多。 若我们要打包序列$A[k],B[k]$。 先设$P[k]=A[k]+iB[k],Q[k]=A[k]-iB[k]$ 我们可以表示为 $P[k]$:  $Q[k]$:  于是很显然的有： $A[k]=\frac {P[k]+Q[k]} 2,B[k]=i\frac {Q[k]-P[k]} 2$ . 变换后： $F_P[k]$:  $F_Q[k]$:  于是也有 $F_A[k]=\frac {F_P[k]+F_Q[k]} 2,F_B[k]=i\frac {F_Q[k]-F_P[k]} 2$ 这个等式其实是转换前后始终满足的：因为DFT/IDFT就是对这些向量做**旋转**，而$P,Q$的变换方向我们保证是**一致的**。 . 现在我们要简化$F_P,F_Q$为一次DFT运算。 我们观察$F_P[N-1]$和$F_Q[N-1]$的值，发现它们正好是一对**共轭复数**；而我们又联想到$F_P,F_Q$的定义，便可以推测它们的值序列一定存在着某种**共轭映射**。 于是我们尝试做$F_Q[k]$的共轭复数  这不就是将$P$反向旋转么？！ 即$F_Q[k]=conj(F_P[N-1-k])$（$conj(z)$为取$z$的共轭复数）（**且对于$k=N-1$特别的有**$F_Q[N-1]=F_P[N-1]$，这就是我们上文提到的）。 . 整理一下： 1.设$P[k]=A[k]+iB[k],Q[k]=A[k]-iB[k]$ 2.有$F_A[k]=\frac {F_P[k]+F_Q[k]} 2,F_B[k]=i\frac {F_Q[k]-F_P[k]} 2$ 3.有$F_Q[k]=conj(F_P[N-1-k]) (0<k<N-1)$（$conj(z)$为取$z$的共轭复数）, $F_Q[N-1]=F_P[N-1]$ 这样我们便成功将DFT也两两打包为一次。 . _在[这篇论文](http://www.doc88.com/p-4912876232995.html)中其实有用式子推导DFT的打包的方法。但这种方法没有一个明显的几何直观，这里就不展开了。_ #### I.III 代码示例 为了简洁，这里仅给出打包的核心代码： ```cpp //该程序计算了(A+B)*(C+D)（线性卷积），为了减小精度负担而拆开计算 //cp是复数类，ll是long long，conj()是取共轭复数 //MAXN2较小，MAXN较大 //变量缺省值均为 0 int N1 /*A, B长度*/, N2 /*C, D长度*/, A[MAXN2], B[MAXN2], C[MAXN2], D[MAXN2]; int N /*FFT规模。注意预处理时 N的具体大小(为两个序列相乘长度 )*/, M /*模数*/; ll ans[MAXN]; cp P[MAXN], Q[MAXN], A2[MAXN], B2[MAXN], C2[MAXN], D2[MAXN]; inline void conv(){ /*这里P打包A, B; Q打包C, D*/ for(int i =0; i < N1; ++i) P[i].real =A[i], P[i].imag =B[i]; for(int i =0; i < N2; ++i) Q[i].real =C[i], Q[i].imag =D[i]; preFFT();/*预处理*/ FFT(P), FFT(Q); for(int i =0; i < N; ++i){/*提取，这里把除法整合了，注意前面的cp(0, 0.5)|cp(0.5, 0)*/ cp P2 =conj(P[(N-i)&(N-1)]); A2[i] =cp(0.5, 0)*(P[i]+P2), B2[i] =cp(0, 0.5)*(P2-P[i]); cp Q2 =conj(Q[(N-i)&(N-1)]); C2[i] =cp(0.5, 0)*(Q[i]+Q2), D2[i] =cp(0, 0.5)*(Q2-Q[i]); } /*这里P打包ac,ad，Q打包bc,bd*/ /*程序中有简化运算*/ for(int i =0; i < N; ++i) P[i] =A2[i]*(C2[i]+cp(0, 1)*D2[i]); for(int i =0; i < N; ++i) Q[i] =B2[i]*(C2[i]+cp(0, 1)*D2[i]); IFFT(P), IFFT(Q); for(int i =0; i < N1+N2-1; ++i){/*最后检索虚实数轴即可*/ ll ac =(ll)(P[i].real/N+0.5)%M, ad =(ll)(P[i].imag/N+0.5)%M, bc =(ll)(Q[i].real/N+0.5)%M, bd =(ll)(Q[i].imag/N+0.5)%M; ans[i] =(ac+bd+ad+bc)%M; } } ``` #### \*I.IV 对于复数轴特殊性的一种代数体现 文章完成后，我又在评论区（ *感谢* [*command_block*](https://www.luogu.com.cn/user/58705) ）发现了一种非常初等的针对计算 `(A+B)*(C+D)=AC+BD+AD+BC` 的右边四项（通常是为了防止炸精度）的优化方法： $(A+Bi)*(C+Di)=(AC-BD)+(AD+BC)i$ $(A-Bi)*(C+Di)=(AC+BD)+(AD-BC)i$ 这样我们只需计算三次DFT和两次IDFT就可以加减计算出答案了（而上文介绍的方法只需四次，两者间仅差一次）。 具体实现可以参考[原文](https://www.luogu.com.cn/blog/command-block/solution-p4245#)（是一篇MTT的题解） - ### II.预处理单位根+指针 这是一个很好用的技巧：**预处理单位根**，然后将代码中访问数组下标的操作**替换为指针**，可以节省不少时间。 #### II.I 核心代码 这里直接给出对于FFT示范代码（NTT可以照这样子写）： ```cpp //cp为复数类 inline void FFT(cp *f/*操作序列*/){ for(int i =0; i < N; ++i) if(i < S[i]) swap(f[i], f[S[i]]); /*w的步长，N为总范围，除掉操作范围2，之后每次除2即可（操作范围在增倍）*/ for(int lim =1/*操作范围的1/2*/, step =N>>1; lim < N; lim <<=1, step >>=1){ for(int s/*起始点*/ =0; s < N; s +=lim<<1){ cp *w =W, *f0 =f+s, *f1 =f+s+lim;/*f0为左区间，f1为右区间*/ for(int i =0; i < lim; ++i){ cp res =*f1 * *w; *f1 =*f0-res, *f0 =*f0+res;//这里的运算顺序按照值变化的顺序有所改变 w +=step, ++f0, ++f1; } } } } ``` #### II.II 预处理 FFT单位根预处理是这样子的： ```cpp for(int i =0; i < N; ++i) W[i] =cp(cos(2*pi*i/N), sin(2*pi*i/N));//直接对于最小的步长求单位复根 //这里也可以递推求，但精度会下降很多 ``` 做IFFT时单位根，可以直接**对原单位根序列取共轭复数**即可。像这样： ```cpp inline void op(){ for(int i =0; i < N; ++i) W[i] =cp(W[i].real, -W[i].imag); } ``` . NTT单位根的预处理也大同小异。 但做INTT的单位根时，是可以直接从 NTT的单位根 推得的。 我们可以发现这样一个式子：$\omega^i\cdot\omega^{n-i}\equiv1\pmod P$。 也就是说，$\omega^i, \omega^{n-1}$**互为乘法逆元**。即可以说$\omega^{n-1}\equiv\omega^{-i}\pmod P$ 直观的理解的话可以参考上文中**打包FFT时那两个序列**的关系（$P, Q$）。 代码如下： ```cpp inline void getWinv(){ for(int i =0; i < N; ++i) Winv[i] =W[(N-1)&(N-i)]; } ``` #### II.III 实践经验小结 做题时，选择预处理的方式也是有技巧的。若题目**精度**要求不大，递推处理即可；否则一定要每一项都根据公式算，避免精度误差。同时**打包操作也是会一定地影响精度的**。 由于打包FFT的效率极高，个人还是推荐用公式保证精度再加打包；实在不行拆系数算也是比常规算法效率要高的。 - ### III. [已编辑] 这条被我验证过其实是无效的QWQ。 . 但为什么还要占一个位子呢？ 原本我分享的是利用 `double` 的玄学精度，尽量地让操作数字转变为小数来提高精度。 但实际上 `double` 是按 **数字长度** 做精度的，和 **数字大小** 则没关系。 ~~窝太菜re呜呜呜~~ - ### IV. 缩减规模 如像上面提到的"常见卷积形式"，需要求A`[0, n-1]`与B`[0, 2n-2]`的卷积C`[0, 3n-2]`的`[n-1, 3n-3]`项，常规操作是把A, B全部扩充到`[0, 3n-3]`然后再做卷积，但实际上只需要扩充到`[0, 2n-2]`就行了。 . 我们观察下线性卷积的图解发现：使用`"0"`最多的是答案序列的最右项，且正好"用"完了所有的`"0"`（如果序列长度不一的话，上下统计结果也是不一样的。**这里取剩下最少的一组**）：  而从最右项从右往左，"用"的`"0"`的数量也会逐个减少。 假设我们减少一位序列的长度，正好"用"完`"0"`的那项就会受"干扰"。 也就是说，假如我们逐个减小序列的长度，就会逐次地从右往左地影响答案。注意这里**原有序列**是不变量，这个结论要求减少的项是扩充的`"0"` 若我们只做`[0, 2n-2]`，就相当于不做末$n-1$项，也就是减少了$n-1$项，于是从右往左地影响了`[0, n-2]`，正好保证了答案不受影响。 代码就如上文提到的一样。 . 具体的式子也是可以列出来的，但没有图像直观，这里就不多提了。 - ### V. 任意长度循环卷积和Bluestein's Algorithm 有时候，题目会要求我们做任意长度的循环卷积，但显然我们不能用朴素的算法。我们需要一个复杂度和FFT近似的算法才能通过评测。 #### V.I 利用卷积性质 我们再做一次循环卷积的示意图，但这次并不做多余的序列。  我们发现，它其实是以答案的位置为分界点，把序列分成了 包括答案位置的右边 和 不包括答案位置的左边 两段，并分别做了一个"**对称**"的特殊乘法。 接下来我们再观察为了**做线性卷积**而拓展后的答案序列的 `[x]`项 与 `[x+n]`项 ： `[x]`：  `[x+n]`：  我们发现，**`[x]`项所缺的部分正好就是`[x+n]`所有的部分!** - ##### V.I.I `[x]`项组成 首先右边序列贡献不变是可以保证的。 我们再看左边序列：发现扩充时补零的长度至少为$n-1$，这样就保证了**左边序列的贡献为零**（可以想象一下）。 - ##### V.I.II `[x+n]`项组成 首先它的左边序列一定是全由`"0"`组成的，贡献也为零。 然后我们发现，若从右往左计数，至`[x+n]`项正好有$x+1$个`"0"`，这和`[x]`项右边序列的长度相等。也就是说，**保证了`[x]`右边序列的元素不做贡献。** 而对于中间那段，因为这种特殊乘法的"**对称性**"，中间的贡献和原来`[x]`项左边序列**的贡献是等价的。**（这里也建议想象一下） 另外，特殊的对于$[x=n-1]$，不需要`[x+n]`项补充。（可自证） - ##### V.I.III `[x+kn]`项的存在性 补足的长度**至少**为$n-1$；而对于$(x=n-1)$，是不需要`[x+n]`项补充的（可观察图）。 对于`[x+kn]`$(k>1)$的情况，我们可以置之不理，**实际上`[x+kn]`的贡献也是为零**的（这时从左往右计数的"0"数量已经超过原有序列长度了）。 但要注意，**保证存在的项只有`[n, 2n-2]`，实际使用时还是要特判的。** - ##### V.I.IV 实现 我们只需做一次 线性卷积。**然后仅将`[x+n]`项与`[x]`项叠加**即可（特别的对于`[n-1]`项不考虑）。 代码如下： ```cpp //此程序可以快速计算任意长度的 循环卷积 //cp是复数类 //n是序列规模，N是FFT的规模 //MAXN2较小，MAXN较大 //变量缺省值均为 0 cp P[MAXN], Q[MAXN]; int ans[MAXN2]; inline void Conv(int *f, int *g){ for(int i =0; i < n; ++i) P[i].real =f[i], Q[i].real =g[i]; preFFT()/*预处理*/, FFT(P), FFT(Q); for(int i =0; i < N; ++i) P[i] =P[i]*Q[i]; IFFT(P); for(int i =0; i < n-1; ++i) ans[i] =int(P[i].real/N+0.5)+int(P[i+n].real/N+0.5); ans[n-1] =int(P[n-1].real/N+0.5);/*为了避免越界特判了对于 [n-1]项的计算*/ } ``` . 这种算法其实已经很优秀了。在不要求对DFT（未IDFT）的序列做操作时，它甚至**是最优的**。 但一旦题目要求做如序列（循环卷积）幂，这个算法最终会多出一个$O(log n)$的复杂度（封装），这是我们不想要的。 因此，我们需要一个**切实**能做 **任意长度DFT** 的算法。 #### V.II Bluestein's Algorithm 上文中我们曾提到 CZT算法。 $Z[k]=W^\frac {-k^2} 2\sum\limits_{n=0}^{N-1}A^{-n}W^\frac {-n^2} 2x[n]\cdot W^\frac {(k-n)^2} 2$ 若设$A=1,W=\omega$，这个式子**就变成 DFT 的形式了**： $F[k]=\omega^\frac {-k^2} 2\sum\limits_{n=0}^{N-1}\omega^\frac {-n^2} 2x[n]\cdot \omega^\frac {(k-n)^2} 2$ 而对于 IDFT，我们只需改下符号重新推导，最后加个$\frac 1 N$即可： $x[k]=\frac 1 N\cdot\omega^\frac {k^2} 2\sum\limits_{n=0}^{N-1}\omega^\frac {n^2} 2F[n]\cdot \omega^\frac {-(k-n)^2} 2$ 这即是 Bluestein's Algorithm。 ##### V.II.I 兼容NTT 但是，这种拆法在做NTT时是行不通的：我们无法保证幂数乘$\frac 1 2$后仍是整数。 所以我们可以拆成均是**相邻整数相乘**的式子： $nk=\frac {(k-n)(k-n+1)-k(k +1)-n(n-1)} 2$ （这个式子其实也可用组合数解释，具体的话就是$nk=\binom{n+k}2-\binom n2-\binom k2$） 然后直接带入即可： $F[k]=\omega^\frac {-k(k+1)} 2\sum\limits_{n=0}^{N-1}\omega^\frac {-n(n-1)} 2x[n]\cdot \omega^\frac {(k-n)(k-n+1)} 2$ $x[k]=\frac 1 N\cdot\omega^\frac {k(k+1)} 2\sum\limits_{n=0}^{N-1}\omega^\frac {n(n-1)} 2F[n]\cdot \omega^\frac {-(k-n)(k-n+1)} 2$ ##### V.II.II 代码示例 (NTT版) 设$g(x)=\omega^{op\cdot\frac {-x(x-1)} 2}A[x],f(x)=\omega^{op\cdot\frac {x(x+1)} 2}$，其中$op$是外部变量，作用是**控制转换方向**；$A[x]$是指当前转换"对应的序列"。 便有： $F[k]=\omega^\frac {-k(k+1)} 2\sum\limits_{n=0}^{N-1}f(k-n)g(n)$ $x[k]=\frac 1 N\cdot\omega^\frac {k(k+1)} 2\sum\limits_{n=0}^{N-1}f(k-n)g(n)$ （注意这份代码也使用了 **6.IV** 的优化） ```cpp //此程序可以快速计算任意长度的 DFT/IDFT；此处实例为NTT //cp是复数类，ll是long long，Pow()是快速幂 //gg是原根，M是模数，n是规模，invn是模 M意义下 n的逆元，N是FFT的规模 //MAXN2较小，MAXN较大 //变量缺省值均为 0 inline int G(int x, int op, int *f){ return Pow(gg, op*-x*(x-1)/2)*f[x]%M; } inline int F(int x, int op){ return Pow(gg, op*x*(x+1)/2); } int f[MAXN2], g[MAXN2]; cp P[MAXN], Q[MAXN]; inline void CZT(int *A, int op /*1为正变换，-1为逆变换*/){ for(int i =-n+1; i <= n-1; ++i) f[n-1+i] =F(i, op); for(int i =0; i < n; ++i) g[i] =G(i, op, A); for(int i =0; i < 2*n-1; ++i) P[i].real =f[i], Q[i].real =g[i]; preFFT()/*预处理*/, FFT(P), FFT(Q); for(int i =0; i < N; ++i) P[i] =P[i]*Q[i]; IFFT(P); for(int i =n-1; i < 2*n-1; ++i){ A[i-n+1] =(ll)(P[i].real/N+0.5)%M*Pow(gg, op*-i*(i+1)/2); if(op == -1) A[i-n+1] =A[i-n+1]*invn%M; } } ``` 代码中有**很多部分是可以预处理的**。但为了简洁，这里就用最朴素的版本了。 --- ## 后记 在审核回复我前一篇文章时，我发现除了一个通过的，只有我的文章没有明确不予通过QWQ。 所以就心血来潮又写了一篇。（ _现在看来和上一篇已经有天差地别的区别惹w_ ） 这篇和上一篇不同的是，里面的概念和技巧都是较深奥的，可能不适合FFT初学者看。但考虑到前面**已经有**入门教程，想让自己的文章上日报还是要拿出一些干货的。 本文最有价值的内容应该是最后一节的优化技巧。里面我整理了尤其是在 《再探快速傅里叶变换》--毛啸 里提到的打包技巧，和我平日里总结的优化方式。（~~有了这些一定可以让您的FFT爆踩NTT~~） 最后再扯件事。本来文章开头是想从最最基本的 FT,傅里叶变换 讲起的。但中途查了许多资料，最后也没有完全理解连续意义下的傅里叶变换，还险些推翻了之前对DFT的认识（~~要是真推翻了重写让我这个鸽子怎么办~~），于是最后这部分还是咕咕了（~~你可以试着读下本文链接，猜下原来的标题是什么QWQ~~） 希望最终可以通过审核把，毕竟我的博客太惨了QAQ（~~真·一个人都没~~） . _现在通过审核了www。_ _FT部分其实多少也快扯到了，不过当然仅仅是其冰山一角QAQ。_ _文章为了通俗易懂舍弃了部分严谨性，一些分析也没有采用更系统的办法（典型的如卷积那部分）。这点我在写文时就想到了，不过也无可奈何（有兴趣的读者可以百度了解一下信号处理的领域，那是一个很大的坑，且基本与竞赛无关）。_ _能看到这里，想必您已经读完全文了把。最后感谢您的阅读~QWQ 别忘了点赞_ $\color{#ffb464}\texttt{首次校对于2020/01/17/19:37}$ $\color{#ffcc64}\texttt{最后校对于 2020/02/11/16:54 (添加了卷积相关) }$ $\color{#ffb400}\texttt{最后修改于 2020/03/06/11:56}$