零知识证明 - 理解FFT的蝶形运算

好几个星期没写文章，主要原因是和小伙伴们一起优化bellman的性能。新的，有趣的，让人兴奋的优化想法一个个蹦出来，性能一点点的提升，让原先的不可能变成了一个个的可能。性能优化，需要沉淀，需要专注分析，需要实践，需要集思广益。更奇妙的是，一步步的优化，之前感觉没有必要优化的小点，都变的重要起来。 Trapdoor Tech的小伙伴们相信：持续，专注，才有好的品质。

这个周末有空，讲讲我对FFT的理解。方便学习零知识证明的小伙伴理解这部分内容。从头讲起，零知识证明Groth16算法基于QAP问题：

零知识证明 - Groth16算法介绍

利用Groth16计算证明之前，需要计算出H。目前，普遍采用的是FFT算法。

01 FFT计算原理

学习FFT计算原理，还是要看算法导论（Introduction-to-Algorithms）的第30章（Polynomials and the FFT)。算法导论，讲东西，来龙去脉讲的比较清楚，逻辑性比较强，也比较细致。该章节从多项式的表示开始，多项式可以用两种表示方法：系数表示和点值表示。

系数表示，对多项式加法友好，时间复杂度是O(n)，但是多项式乘法不友好，采用多项式一项项的相乘，时间复杂度为 O(n**2)。点值表示，对多项式乘法友好，时间复杂度是 O(n)，但是多项式加法不友好。

DFT，采用特殊的点，计算值。这些特殊的点，就是单位根的幂次。采用这些特殊点的情况下，系数表示和点值表示之间可以通过FFT（快速傅立叶变换）计算完成。时间复杂度是O(nlogn)。

算法导论上这个图，可以很好的解释这些计算之间的关系：

如果是系数表示，多项式乘法的计算复杂度为O(n**2)。如果先转化为点值表示，点值表示相乘，再转化回系数表示，计算复杂度为O(nlogn)。

n次单位根

存在n个n次的单位根。这些单位根满足w**n = 1。

这些单位根，存在一些性质或者定理：

FFT的基本原理

DFTn(a)的计算采用FFT的计算，时间复杂度为O (nlogn)。原理是，将多项式分解成奇偶两部分（假设n为2的幂次，n不为2的幂次也有相应的算法，后话）：

也就是说，计算A(x)的对应的值，可以分解成两个子多项式（奇偶），每个多项式的点是之前的平方。递归的角度可以看出，这些子多项式，可以进一步划分分更小的多项式。整个计算的复杂度：

如果把输入点也分成两部分，每部分是n/2个，k的范围是0~n/2，则:

这个就是传说中的“蝶形计算”。当两个子多项式的值计算完成后，原来多项式的两个相差n/2的点的值，是两个子多项式的值的一加一减。

整个计算算法的伪代码如下：

其中的作用在第二个子多项式结果上的w_n^k，称为旋转因子。整个蝶形运算可以用可视化的图形表示如下：

8阶多项式示例

算法导论清晰的给出了8阶多项式的FFT的计算过程：

在stage s=1的情况下，也就是最“底层”的两个1阶多项式合并。旋转因子是w_2^0。注意这一层的输出结果，会作为下一层stage s=2的输入，这一层的旋转因子是w_4^{0,1}。这一层要处理的是“两个”2阶多项式的合并。在奇偶划分的情况下，第一个/第二个元素和第三个/四个元素（跨度为2）做蝶形运算。依次类推，直到stage s=3 做完。

整个过程，可以通过Merkle树形象表示：

注意最底层的元素的位置，经过递归的划分已经不是顺序的。算法导论也总结了相应的计算方法，感兴趣的小伙伴，可以自行查看。

可以返回头，体会体会这个公式：

因为特殊的取点的原因（w^2相当于阶降低了一半），子多项式的值也是结果多项式的一半。因为在取点分一半（恰好是子多项式的阶），又因为旋转因子正好是正负，所以存在蝶形关系。每一层的子多项式的元素翻倍。

02 FFT一般性计算扩展

在理解了2分的FFT的计算逻辑后，可以扩展到一般性计算。推荐另外一个网站（注意多项式符号和上文有点区别）：

http://www.bealto.com/gpu-fft2_fft.html

多项式不再拆解成“奇偶”两部分，而是分成P份，每份为Q个元素：

其中，x_u为系数。u的范围是0-P-1，先定义一下，在这种划分下的DFT计算：

DFT_Q代表了有Q个元素，x[u]代表了Q个元素序列，每个元素序列跨度为P，P*Q=n。

如果把输入也切割成P份（就像二分情况下，将输入分成两份一样），k可以表达成k=v+Qj，v=0..Q-1，j=0..P-1。其中，DFT_Q (x[u])记为z[u]。

也就是说，P*Q分组的情况下，蝶形运算如下图：

进一步，可以将y的每个元素的计算看成一个新的DFT_P：

在理解这些的前提下，看看网站给出的按照P*Q分组计算的示意图就比较清晰：

整个FFT的计算由三部分组成：

1/ P个DFT_Q计算（z[0], z[1], ...z[P-1]) 2/ 乘上Twiddle 3/ Q个DFT_P计算

在Q个DFT_P计算之后，需要将P个计算数据“分散”到以Q为偏移的具体的位置。

总结：

零知识证明Groth16的计算中为了计算H，采用了FFT计算。FFT的蝶形运算是理解FFT计算的基础。二分情况下的计算，相对清晰，算法导论给出了详细的推导过程。P*Q的一般性计算推导，相对复杂一些。