原文：https://rje.li/24-07-03-alg-fft.html
作者：Eli Riggs
译者：Kurt Pan
FFT（“快速傅立叶变换”）是有史以来最伟大的算法之一，但大多数人只见过它的经典形式，即将时域中的离散样本转换为频域中的复数表示。它的同类算法 NTT（“数论变换”）对模数上的多项式插值，虽然更少为人所知，但仍在广泛使用，主要用于密码学。
然而，为了更特殊的目的、更陌生的领域和更不熟悉的性质，人们发明了更加奇异的“FFT”。在本文中，我将说明任何一种 FFT 背后的那一个单一算法，然后快速介绍这些奇异的 FFT，并展示如何通过在我们的共用算法中插入一小组参数来实现它们。
最终我们会看到，FFT 让我们可以在 的时间内更改一个在 个点上定义的函数的基。朴素的方法是通过应用线性变换矩阵来改基，而矩阵乘法需要 ~ 时间，因此这是一个巨大的优势。
在继续之前，先说几点：
我将一组允许我们应用 FFT 算法 的参数称为一个FFT 。这起初有点令人困惑，但很快就会变得不那么含糊也更方便。为清楚起见，我仅将复平面上众所周知的 FFT 称为 经典 FFT。 始终表示 ，通常是在指数变得太小的时候。 我使用了很多下标，因此向量索引由上标 表示。下标 表示标签或参数。
本文是用 Sage 笔记本写的。我写了一些实用函数，它们的含义应该从名称中就能看出来。唯一奇怪的是bind——它让我可以在定义之后向类添加方法，这样我就可以在文章中间创作了。
%display latex
from collections import namedtuple, defaultdict
from itertools import islice, chain
take = lambda n, xs: list(islice(xs, n))
interleave = lambda *xss: list(chain.from_iterable(zip(*xss)))
def iterate(x, f): yield x; yield from iterate(f(x), f)
powers = lambda base, start=None: iterate(start or base.parent().one(), lambda x: x*base)
dot = lambda *xss: sum(product(xs) for xs in zip(*xss))
rands = lambda f, n: [f.random_element() for _ in range(n)]
bind = lambda cls: lambda f: (setattr(cls, f.__name__, f), f)[1]
一个FFT的定义
对于给定域 和大小 ，一个对 和 的 FFT 由其包含 个点的初始定义域 ，和一个 层的序列：
完全确定。
Layer = namedtuple('Layer', 'π t')
Fft = namedtuple('Fft', 'domain layers')
其中，要么 ，意味着 只有一个点， 为空，要么：
是一个函数，将 的点映射到一个大小为原来一半的新定义域 上，而且它是二对一的，意味着对于任意输出点 ，恰好有两个输入点 映射到它：。 是一个从 到 的函数，它对映射到 中同一点的 中的任意两点进行区分，或者输出一个不同的值。 也是一个 FFT，其中 ，或者用python写为 layers[1:]。
根据该定义，我们得到一个由 个由映射 和除 之外的每个定义域上的区分器（通常称为 twiddles）生成的定义域 组成的序列，每个定义域的大小减半，
显式地构造出逆映射 将会很有用，该映射将 中的每个输出点映射到其两个输入 上（又名其 纤维），以及 和 的值。让我们为此向 Fft 添加一个方法，它还将检查我们的要求是否得到满足。返回的字典具有以下形式 ，其所有键的集合恰好是 。
@bind(Fft)
def build_inverse_map(self: Fft):
π, t = self.layers[0]
inv_map = defaultdict(lambda: [])
for x in self.domain:
inv_map[π(x)].append((x, t(x)))
assert all(
len(fiber) == 2 and fiber[0][1] != fiber[1][1]
for fiber in inv_map.values()
)
return inv_map
如果满足这些要求，我们可以将 FFT 算法 应用于任意函数 来得到一个由 个系数组成的向量 ：
给定任意一个 FFT，我们可以导出其 基 ，它是 的 函数 向量。然后我们就可以理解系数 的含义了：它们是基函数的 权重，而基的加权和将重构出我们的原始函数 。
嗯……如果你的 FFT 构成的基是具有足够的维度来重构原始函数的，那么它应该如此。我们假设这是真的，因为我们只对应用 FFT 算法感兴趣，而不是结果到底有多有用。
我将解释如何导出这个基，因为它是由 FFT 的结构形成的，但需要注意的是，基纯粹是概念性的：它告诉你如何解释得到的系数，但如果你只是想运行算法，它并不必要。
FFT 算法及其基
给定 （ 上的一个函数），FFT 算法有三个步骤：
使用当前层 将 分解为两个较小的定义域 上的函数 和 ，满足：
在 和 上递归调用 以得到它们的系数。 以“合理的方式”组合 和 的系数，并返回这些组合系数。
在这些步骤中，分解是最容易解释的，但却最难说清楚。现在假设我们会编写一个函数decompose，该函数输入逆映射和函数在 上的求值，返回 和 在 上的求值。
那么除decompose之外，该算法可以完整描述为：
@bind(Fft)
def fft(self: Fft, f: dict):
if not self.layers:
return [f[list(self.domain)[0]]]
inv_map = self.build_inverse_map()
f0, f1 = decompose(inv_map, f)
next_fft = Fft(inv_map.keys(), self.layers[1:])
return interleave(next_fft.fft(f0), next_fft.fft(f1))
我知道这还没什么意义，我只是想把它写出来，因为它非常短。
我们会从头开始，使用上面列出的内容，尝试得出相同的解决方案。
我们的目标是取一个函数 并返回一个系数向量 作为基函数 的权重，使得加权和可以重构 。（在此过程中，我们还需要定义 ）。所以我们需要定义这个等式的左边：
在 的基础情况下，定义域中只有一个点 。 只是将 映射到一个单个值，因此可以用一个始终返回该值 的常数函数来完整描述之。将我们的单个基函数定义为 ，可以返回 作为该函数的系数或权重，从而满足我们的标准：。
否则，继续对 进行操作。首先，将 decompose 函数应用于 ，将其拆分为 和 ，或者“依赖于 的部分”和“不依赖于 的部分”。
因此，我们从每个 处的 值开始，现在有了每个 处的 和 值，或者等效地说，每个点 处的 和 值。
现在开始递归，但我们需要小心。我们是将对 和 应用下一个（更小的）FFT，而不是 和 。我的意思是，例如对 ，更小的 FFT 将返回由较小定义域的函数组成较小的基 的系数 。为了将其代回我们的方程，必须移动内部的 以将定义域投影到更小的定义域，对于 也是如此：
再次地 - 我们将更小的 FFT 应用于 和 ，为了使结果准确地重构 ，我们需要将更小的基函数应用于 而不是 。
现在我们试着推导 和 的表达式，以使得等式始终成立。我们会试着将右侧表示为单个求和，因此首先将 移入内部：
但是……在一般情况下，我们能达到的也就只有这些了。我们有两个不同的内积式的求和，而构造 和 以使求和始终相等的唯一方法就是直接连接，设置 ，对于 也类似，所以我们就是这么做的。
差不多吧。实际上，我们选择交错 和 的元素，但这种选择或任何其他选择都是任意的。我们通常选择交错的原因是 通常是 1 次函数，而 通常是 2 次函数。因此，当我们交错或将 中的项放入偶数索引并将 中的项放入奇数索引时：
于是索引将对应于基函数的次数，假设对 也是如此。换句话说，假设 的次数为 ，则 的次数将为 ，而 的次数将为 ，与我们的外部索引 完全对应。
掌握了以上所有数学知识后，你现在应该可以理解我们算法的最后两行了：
next_fft = Fft(inv_map.keys(), self.layers[1:])
return interleave(next_fft.fft(f0), next_fft.fft(f1))
我们还可以添加一种方法，将 FFT 的每个基函数应用于一个点，并将结果作为向量返回。这很方便，因为我们现在可以使用dot(coeffs, fft.basis(pt))在任何点对系数进行求值。
@bind(Fft)
def basis(self: Fft, pt):
if not self.layers:
return [1]
domain, ((π, t), *next_layers) = self
next_fft = Fft({ π(x) for x in domain }, next_layers)
next_basis = next_fft.basis(π(pt))
return interleave(
next_basis,
[t(pt)*b for b in next_basis],
)
此时定义逆 FFT 也非常简单，它将系数转换为求值：
@bind(Fft)
def ifft(self: Fft, coeffs: list) -> dict:
if not self.layers:
return { x: coeffs[0] for x in self.domain }
domain, ((π, t), *next_layers) = self
next_fft = Fft({ π(x) for x in domain }, next_layers)
f0 = next_fft.ifft(coeffs[0::2])
f1 = next_fft.ifft(coeffs[1::2])
return {
x: f0[π(x)] + t(x) * f1[π(x)]
for x in domain
}
分解函数
在 FFT 的每一层，我们将函数 （定义在 上）分解为两个函数 和 ，每个函数定义在 上：
为此我们需要利用 和 的特殊性质。考虑 中映射到 中的同一点的两个点 和 ，
我们知道 和 的值，因为 （函数或其所有求值的向量）是作为算法的输入给出的，并且我们知道 和 ，因为 是我们 FFT 的一个参数。所以这里有 4 个未知数，这是无解的……但是等等！ 和 映射到同一点，这意味着 ，我们将其简单地写为 。重写为：
现在我们有 2 个方程和 2 个未知数，可以求解了。如果你仔细看看，你可能会认出这是 ，其中 和 ，所以你可能已经知道了如何手动求解之。在我看来，将其视为矩阵乘法会更容易，我们使用 2x2 求逆公式：
这为我们提供了 和 的显式公式，仅和可以访问的 和 相关。你可以看到为什么我们需要 来区分 的纤维，从而保证 ，否则我们的矩阵将不可逆，从而导致出现除以零。
或者我们也可以让 Sage 为我们完成这件事：
f_x0, f_x1, t_x0, t_x1 = SR.var('f_x0 f_x1 t_x0 t_x1')
(matrix([[1, t_x0],
[1, t_x1]]).inverse() * vector([f_x0, f_x1])).simplify_rational()
注意到 的表达式对应于计算升降的“斜率” ，而 的表达式对应于“y轴截距”，如果你将其重新带入“”，则会得到该值。
如果我们预先计算定义域内所有点处标量的逆 ，就可以通过对的2 次乘法 + 2 *（1 次减法和 1 次乘以逆标量）= 4 次乘法和 2 次减法。
实践中，你会希望选择一个巧妙的twiddle 和巧妙的定义域来减少这种计算。例如，如果我们知道域中的所有点都有 ，就像在乘法的情况下一样，那么它的开销会低得多：
_.substitute(t_x1 = -t_x0).simplify_rational()
这可以通过 1 次加法、1 次减法和 2 次乘法来计算。如果将缩放 延迟到最后，则可以将其减少为 1 次乘法，再加上一次最终的缩放，最终开销更低（ 次总缩放乘法，而不是 ）。
因此，这里是decompose函数：
def decompose(inv_map, f):
f0 = {
π_x: (f[x0] * t1 - f[x1] * t1) / (t1 - t0)
for π_x, ((x0, t0), (x1, t1)) in inv_map.items()
}
f1 = {
π_x: (f[x1] - f[x0]) / (t1 - t0)
for π_x, ((x0, t0), (x1, t1)) in inv_map.items()
}
return f0, f1
再来一遍，但写在一起
如果将我们的三个函数合并为一个，它看起来会是这样的。我发现对这段代码片段沉思会让人平静。将其打印出来并塞到枕头下面以求得好运。
@bind(Fft)
def fft(self: Fft, f: dict):
if not self.layers:
return [f[list(self.domain)[0]]]
domain, ((π, t), *next_layers) = self
inv_map = defaultdict(lambda: [])
for x in domain:
inv_map[π(x)].append(x)
f0 = {
π_x: (f[x0] * t(x1) - f[x1] * t(x0)) / (t(x1) - t(x0))
for π_x, (x0, x1) in inv_map.items()
}
f1 = {
π_x: (f[x1] - f[x0]) / (t(x1) - t(x0))
for π_x, (x0, x1) in inv_map.items()
}
next_fft = Fft(inv_map.keys(), next_layers)
return interleave(next_fft.fft(f0), next_fft.fft(f1))
在每一层，分解都是线性于，然后我们对一半大小的问题进行两次递归。层数为 ，因此总运行时间为 。
矩阵表示
由于应用 FFT 是在基变换，或函数的线性变换，因此可用应用于函数求值的矩阵来表示。矩阵乘法是 （量级）的，因此性能要差得多，但我们可以将这种表示用于其他目的，例如计算当用作一个线性码时 FFT 的最小距离。
我们矩阵的列将来自把 FFT 应用于“独热”向量 ，可以将其视为将求值的拉格朗日基映射到 FFT 基。
one_hot = lambda n, i: vector([1 if j == i else 0 for j in range(n)])
@bind(Fft)
def fft_matrix(self: Fft):
domain = list(self.domain)
N = len(domain)
return matrix.column([
self.fft(dict(zip(domain, one_hot(N, i))))
for i in range(N)
])
例子
每个这些例子都值得写一篇自己的文章，但目前我只想给出一个参考，并展示它们对于其 、 和 的选择是有效的。
def demo_fft(field, sizes, fft_factory, symbolic_point):
for n in sizes:
print(f'n={n}:')
fft = fft_factory(n)
print(f' basis:', fft.basis(symbolic_point))
set_random_seed(0)
evals = dict(zip(fft.domain, rands(field, 2^n)))
# FFT，然后在每个原始点进行求值，得出原始值
coeffs = fft.fft(evals)
evals2 = { x: dot(coeffs, fft.basis(x)) for x in fft.domain }
assert evals == evals2
# IFFT 正常工作
evals3 = fft.ifft(coeffs)
assert evals == evals3
# FFT矩阵正确计算系数
M = fft.fft_matrix()
assert M * vector([evals[x] for x in fft.domain]) == vector(coeffs)
print(' tests succeeded!')
IM = M.inverse().T[:2^(n-1)]
# 这里计算非常慢
dist = LinearCode(IM).minimum_distance()
print(' code minimum distance:', dist, f'(best possible is {2^(n-1)+1})')
乘法 FFT
经典 FFT 在 的乘法子群中起作用，其中 是第 个单位根。由于精度问题，这实际上有点难展示。如果我们用浮点数， 映射就不是精确的二到一的，我们的结果也不会准确。因此，我将展示 NTT，它是完全相同的算法，但应用于有限域的乘法子群。
注意这里的基正是你所期望的那些。
SquaringLayer = Layer(
π = lambda x: x^2,
t = lambda x: x,
)
GF17 = GF(2^4 + 1)
def make_ntt(n):
gen = GF17.multiplicative_generator() ^ (2^(4-n))
return Fft(
domain = take(2^n, powers(gen)),
layers = [SquaringLayer] * n,
)
demo_fft(GF17, [2, 3], make_ntt, polygen(GF17, 'X'))
n=2:
basis: [1, X, X^2, X^3]
tests succeeded!
code minimum distance: 3 (best possible is 3)
n=3:
basis: [1, X, X^2, X^3, X^4, X^5, X^6, X^7]
tests succeeded!
code minimum distance: 5 (best possible is 5)
圆 FFT
来自 Circle STARKs
https://eprint.iacr.org/2024/278
CircleYLayer = Layer(
π = lambda xy: xy[0],
t = lambda xy: xy[1],
)
CircleXLayer = Layer(
π = lambda x: 2*x^2 - 1,
t = lambda x: x,
)
GF127 = GF(2^7 - 1)
C127.<i> = GF127.extension(polygen(GF127)^2 + 1)
circle_gen = C127.multiplicative_generator()^(GF127.order()-1)
def make_cfft(n):
gen = circle_gen^(2^(7-n-1))
return Fft(
domain = take(2^n, powers(gen^2, start=gen)),
layers = [CircleYLayer] + [CircleXLayer] * (n-1),
)
demo_fft(GF127, [2, 3], make_cfft, polygens(GF127, 'X,Y'))
n=2:
basis: [1, Y, X, X*Y]
tests succeeded!
code minimum distance: 2 (best possible is 3)
n=3:
basis: [1, Y, X, X*Y, 2*X^2 - 1, 2*X^2*Y - Y, 2*X^3 - X, 2*X^3*Y - X*Y]
tests succeeded!
code minimum distance: 4 (best possible is 5)
加法 NTT
一开始出自 Novel Polynomial Basis and Its Application to Reed-Solomon Erasure Codes，这里我使用 FRI-Binius中的参数。
https://arxiv.org/abs/1404.3458 https://eprint.iacr.org/2024/504
GF256 = GF(2^8, repr='int')
subspace = lambda n: [GF256.from_integer(i) for i in range(2^n)]
beta = lambda i: GF256.from_integer(2^i)
W = lambda i, x: product(x - u for u in subspace(i))
q = lambda i, x: (W(i, beta(i))^2 / W(i+1, beta(i+1))) * x * (x + 1)
def make_additive(n):
return Fft(
domain = subspace(n),
layers = [
Layer(
# i=i to work around lambda capture quirk
π = lambda x, i=i: q(i, x),
t = lambda x: x,
)
for i in list(range(n))
]
)
demo_fft(GF256, [2, 3], make_additive, polygen(GF256, 'X'))
n=2:
basis: [1, X, 122*X^2 + 122*X, 122*X^3 + 122*X^2]
tests succeeded!
code minimum distance: 3 (best possible is 3)
n=3:
basis: [1, X, 122*X^2 + 122*X, 122*X^3 + 122*X^2, 251*X^4 + 219*X^2 + 32*X, 251*X^5 + 219*X^3 + 32*X^2, 81*X^6 + 81*X^5 + 170*X^4 + 81*X^3 + 251*X^2, 81*X^7 + 81*X^6 + 170*X^5 + 81*X^4 + 251*X^3]
tests succeeded!
code minimum distance: 5 (best possible is 5)
相关文章：
