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

本文为日报特供版QAQ

考虑到前面已经有人写过FFT了，所以部分FFT的基础内容不会出现在本篇，如基本的公式，NTT（它仅涉及单位复根循环的性质）等。

本文主要是讲解一些竞赛中会用到的围绕着FFT算法的一些内容，如DFT过程的几何直观、卷积和*CZT等，另外，最后也会介绍一些冷门却好用的tricks。

之前日报的FFT可以看这里

如果是FFT入门建议看下我之前的文章（有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 所述的优化）：

//常见卷积形式代码 
//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 对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]

我们可以表示为

 $Q[k]$:  于是很显然的有： $A[k]=\frac {P[k]+Q[k]} 2,B[k]=i\frac {Q[k]-P[k]} 2

.

变换后：

 $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也两两打包为一次。

.

在这篇论文中其实有用式子推导DFT的打包的方法。但这种方法没有一个明显的几何直观，这里就不展开了。

I.III 代码示例

为了简洁，这里仅给出打包的核心代码：

//该程序计算了(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 ）发现了一种非常初等的针对计算 (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就可以加减计算出答案了（而上文介绍的方法只需四次，两者间仅差一次）。

具体实现可以参考原文（是一篇MTT的题解）

• II.预处理单位根+指针

这是一个很好用的技巧：预处理单位根，然后将代码中访问数组下标的操作替换为指针，可以节省不少时间。

II.I 核心代码

这里直接给出对于FFT示范代码（NTT可以照这样子写）：

//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单位根预处理是这样子的：

for(int i =0; i < N; ++i) W[i] =cp(cos(2*pi*i/N), sin(2*pi*i/N));//直接对于最小的步长求单位复根
//这里也可以递推求，但精度会下降很多

做IFFT时单位根，可以直接对原单位根序列取共轭复数即可。像这样：

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）。

代码如下：

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]项不考虑）。

代码如下：

//此程序可以快速计算任意长度的 循环卷积 
//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后仍是整数。

所以我们可以拆成均是相邻整数相乘的式子：

+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 的优化）

//此程序可以快速计算任意长度的 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}