原文：https://solvable.group/posts/ecfft/
作者：William Borgeaud
译者：Kurt Pan
本文是关于 Eli Ben-Sasson、Dan Carmon、Swastik Kopparty 和 David Levit 最近发表的一篇论文的。作者在这篇论文中提出了经典 FFT 算法的适用于所有有限域的令人惊叹的新推广。本文将概述该算法及其在 Sage 中的简单实现。我强烈建议阅读该论文以了解更多详细信息和背景。
https://arxiv.org/abs/2107.08473
经典FFT算法
令 为质数, 且 是大小为 的子群 。经典FFT 算法可用于在 中对次数 的多项式 在上求值。注意, 朴素的对 在 上每一点的求值算法需要 次操作。
FFT 的工作原理是将 写作
其中 是 次多项式 的偶数和奇数项系数 。
因此, 给定 和 在 ,上的求值， 我们可以用次操作恢复出 在 上的求值。
需要注意的关键在于 的大小是 的一半，因为 。因此,如果我们用表示 FFT 的运行时间, 我们有以下的递推关系
由此我们可以推断出 。
ECFFT算法
ECFFT 算法的目标是将 FFT 算法推广到不具有大小为的乘法子群的有限域上 , 即的域 。
这个想法既聪明又简单。概述如下:
找到一条椭圆曲线 使得 。令 是大小为的子群， 是 的一个陪集 。 令 是一个次数为2的同源 ，使得 的大小是 的一半。 使用 将 分解为两个较小的多项式 ，使用椭圆曲线 和陪集 ，将 ECFFT 应用于 上。
我现在将更详细地解释这些步骤。用 Sage 给出 Secp256k1 曲线的基域的实现, 即
p = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F
F = GF(p)
寻找曲线
我们需要找到一条曲线 使得 。论文里给出了这样做的算法。然而，暴力搜索对于较小的值 来说也足够有效了。
这是一个通过暴力搜索寻找曲线的脚本
def find_curve(n):
while True:
b = F.random_element()
E = EllipticCurve(F, [1, b])
if E.order() % n == 0:
return E
n = 2^12
E = find_curve(n)
对 ，这可以在大约一个小时内找到一条曲线 。这是我找到的曲线
选择陪集
令 是大小为 的子群 。我们选择一个 的随机陪集 。
gen = E.gens()[0]
G = (gen.order()//n) * gen
R = E.random_element()
H = [R + i*G for i in range(n)]
最后, 令 是 在 上的投影 。 将是我们将在其上执行 ECFFT的大小为的集合。
L = [h.xy()[0] for h in H]
寻找同源
我们需要一个将大小减半的同源。这用 Sage 很容易找到：
def find_isogeny(E, H):
for phi in E.isogenies_prime_degree(2):
if len(set([phi(h) for h in H])) == len(H)//2:
return phi
phi = find_isogeny(E, H)
使用的-坐标， 由2次有理函数 给出。
psi = phi.x_rational_map()
u, v = psi.numerator(), psi.denominator()
由此我们可以得到一条新的椭圆曲线 , 这是 的到达域,一个新的陪集 , 以及一个新的给出的 的子集。
上分解多项式
在经典 FFT 中, 我们能根据两个较小多项式 在较小的集合上的求值写出 在 上的求值 。
我们现在想要通过替换来执行相同的操作，将 替换成 替换成 , 以及平方操作替换成 。
令 是一个次数 的多项式 。那么存在次数 的多项式 ，使得
证明见论文附录 A。证明的思想是，从 次多项式对到次数 多项式的线性映射 是单射, 同时也是双射, 因为它的定义域和到达域具有与-向量空间同样的维度。
计算多项式 不像经典 FFT 那样简单。然而, 可以很容易从 在 上的求值到 和 在 上的求值，反之亦然。
实际上, 给定 使得 , 通过上面的公式写出 在 处的求值，我们得到以下线性关系（令 ）
可以看出矩阵是可逆的, 因此我们可以很容易地从 到 ，反之亦然。这给出了我们稍后将使用的以下有效的一一对应关系：
更重要的是, 矩阵不依赖于 , 因此我们可以预先计算它，并将其重用于 ECFFT 的所有实例。
inverse_matrices = []
nn = n // 2
q = nn - 1
for j in range(nn):
s0, s1 = L[j], L[j + nn]
assert psi(s0) == psi(s1)
M = Matrix(F, [[v(s0)^q,s0*v(s0)^q],[v(s1)^q, s1*v(s1)^q]])
inverse_matrices.append(M.inverse())
EXTEND操作
最后一个难题是 EXTEND 操作。令 是 分别在偶数和奇数索引处的元素 , 。
S = [L[i] for i in range(0, n, 2)]
S_prime = [L[i] for i in range(1, n, 2)]
给定次数的多项式在上的求值，EXTEND 操作计算 在 上的求值。该论文的主要结果是, 有一个 的EXTEND算法。请注意, 一个朴素的算法要通过在上拉格朗日插值恢复的系数 ，然后在 上求值，需要 。
该算法的工作原理如下。 如果 是常数且 在 和 上的求值是一样的。
否则, 从 在 上的求值导出 在 上的求值, 如上节所示。然后应用 EXTEND操作两次以获得 在 上的求值 。最后，恢复 在 上的求值，如上节所示。
def extend(Q_evals, S, S_prime, matrices, inverse_matrices):
n = len(Q_evals)
nn = n // 2
if n == 1:
return Q_evals
Q0_evals = []
Q1_evals = []
q = nn - 1
for j in range(nn):
s0, s1 = S[j], S[j + nn]
y0, y1 = Q_evals[j], Q_evals[j + nn]
Mi = inverse_matrices[n][j]
q0, q1 = Mi * vector([y0, y1])
Q0_evals.append(q0)
Q1_evals.append(q1)
Q0_evals_prime = extend(Q0_evals, [psi(s) for s in S], [psi(s) for s in S_prime], matrices, inverse_matrices)
Q1_evals_prime = extend(Q1_evals, [psi(s) for s in S], [psi(s) for s in S_prime], matrices, inverse_matrices)
return [
M * vector([q0, q1])
for M, q0, q1 in zip(matrices[n], Q0_evals_prime, Q1_evals_prime)
]
运行时间的递推关系与经典 FFT 中的相同, 为 。
上求值多项式
上一节中给出的 EXTEND 操作的算法会是论文中给出的许多高效算法的构建块。
其中之一是 ENTER 操作的算法, 该算法计算次数 的多项式 在 上的求值。这与我们在经典 FFT 中所做的类似。
算法非常简单。按其低系数和高系数分解 为 。对 和 在 上调用 ENTER 得到 在 上的求值。然后调用 EXTEND 两次以获得 在 上的求值。因为 , 我们得到 在 上的求值，并且可以导出 在 上的求值 。
运行时间的递推关系现在为
因为我们必须在递归步骤中调用 EXTEND。因此运行时间 为 ，比经典 FFT 稍差。
在 Sage 之外运行
正如这篇文章 (希望) 展示的, 这些算法在像 Sage 这样的计算机代数系统中实现起来非常简单。然而, Sage 是极其慢的, 这些实现远谈不上速度。
这些算法最酷的一点是对于给定的域 , 我们可以在 Sage 中预先计算所有必要的数据：所有的集合 以及 EXTEND 算法中所使用的（逆）矩阵。
一旦我们有了这些预计算数据, 算法就只使用上的基本操作 , 即加法和乘法。
因此, 实际的算法可以很容易地用 C++ 或 Rust 等快速语言实现, 同时不必用这些语言实现所有的椭圆曲线和同源机制。
最终代码
#p 是 Secp256k1 曲线基域的大小
p = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F
F = GF(p)
# 有关如何找到 a 和 b 使得 2^12 能整除 E 的阶，请看文章。
a, b = 1, 641632526086088189863104799840287255425991771106608433941469413117059106896
E = EllipticCurve(F, [a,b])
log_n = 12
n = 2^log_n
assert E.order() % n == 0
g = E.gens()[0]
G = (g.order()//n) * g
assert G.order() == n
R = E.random_element()
H = [R + i*G for i in range(2^log_n)]
L = [h.xy()[0] for h in H]
S = [L[i] for i in range(0, n, 2)]
S_prime = [L[i] for i in range(1, n, 2)]
def precompute(log_n, S, S_prime, E):
Ss = {}
Ss_prime = {}
matrices = {}
inverse_matrices = {}
for i in range(log_n, -1, -1):
n = 1 << i
nn = n // 2
Ss[n] = S
Ss_prime[n] = S_prime
matrices[n] = []
inverse_matrices[n] = []
for iso in E.isogenies_prime_degree(2):
psi = iso.x_rational_map()
if len(set([psi(x) for x in S]))==nn:
break
v = psi.denominator()
q = nn - 1
for j in range(nn):
s0, s1 = S[j], S[j + nn]
assert psi(s0) == psi(s1)
M = Matrix(F, [[v(s0)^q,s0*v(s0)^q],[v(s1)^q, s1*v(s1)^q]])
inverse_matrices[n].append(M.inverse())
s0, s1 = S_prime[j], S_prime[j + nn]
assert psi(s0) == psi(s1)
M = Matrix(F, [[v(s0)^q,s0*v(s0)^q],[v(s1)^q, s1*v(s1)^q]])
matrices[n].append(M)
S = [psi(x) for x in S[:nn]]
S_prime = [psi(x) for x in S_prime[:nn]]
E = iso.codomain()
return Ss, Ss_prime, matrices, inverse_matrices
# 预计算 EXTEND_S,S' 所需的数据
Ss, Ss_prime, matrices, inverse_matrices = precompute(log_n-1, S, S_prime, E)
def extend(P_evals):
n = len(P_evals)
nn = n // 2
if n == 1:
return P_evals
S = Ss[n]
S_prime = Ss_prime[n]
P0_evals = []
P1_evals = []
for j in range(nn):
s0, s1 = S[j], S[j + nn]
y0, y1 = P_evals[j], P_evals[j + nn]
Mi = inverse_matrices[n][j]
p0, p1 = Mi * vector([y0, y1])
P0_evals.append(p0)
P1_evals.append(p1)
P0_evals_prime = extend(P0_evals)
P1_evals_prime = extend(P1_evals)
ansL = []
ansR = []
for M, p0, p1 in zip(matrices[n], P0_evals_prime, P1_evals_prime):
v = M * vector([p0, p1])
ansL.append(v[0])
ansR.append(v[1])
return ansL + ansR
# 生成一个随机多项式进行测试
R.<X> = F[]
P = sum(F.random_element() * X^i for i in range(1<<(log_n - 1)))
# P 在 S 上求值
P_evals = [P(x) for x in S]
# 结果包括P在S'上的求值
result = extend(P_evals)
assert result == [P(x) for x in S_prime]
