ECFFT算法
- XPTY
- 发布于 2024-10-16 13:20
- 阅读 8
本文是关于 Eli Ben-Sasson、Dan Carmon、Swastik Kopparty 和 David Levit 最近发表的一篇论文的。
作者在这篇论文中提出了经典 FFT 算法的适用于所有有限域的令人惊叹的新推广。本文将概述该算法及其在 Sage 中的简单实现。我强烈建议阅读该论文以了解更多详细信息和背景。
经典FFT算法
令$p$为质数, 且$n = 2^k$是大小为$\mid n \mid p-1,\ \langle w \rangle = H < \mathbb{F}p^*$的子群 。经典FFT算法可用于在$O(n \log n)$中对次数$< n$的多项式$P(X) = \sum{i=0}^{n} a_i X^i$在上$H$求值。注意, 朴素的对$P$在$H$上每一点的求值算法需要$O(n^2)$次操作。
FFT 的工作原理是将$P$写作 $P(X) = P_0(X^2) + X P_1(X^2)$
其中$P_0$,$P_1$是$< n/2$次多项式$P$的偶数和奇数项系数 。
因此, 给定$P_0$和$P_1$在$H^2$,上的求值， 我们可以用$O(n)$次操作恢复出$P$在$H$上的求值。
需要注意的关键在于$H^2$的大小是$H$的一半，因为$H = -H$。因此,如果我们用$F$表示 FFT 的运行时间, 我们有以下的递推关系 $F(n) = 2 F\left(\frac{n}{2}\right) + O(n)$
由此我们可以推断出$F(n) = O(n \log n)$。
ECFFT算法
ECFFT 算法的目标是将 FFT 算法推广到不具有大小为$n$的乘法子群的有限域上 , 即$n\nmid p-1$的域$\mathbb{F}_p$。
这个想法既聪明又简单。概述如下:
- 找到一条椭圆曲线$E(\mathbb{F}_p)$使得$n | #E(\mathbb{F}_p)$。令$G <E(\mathbb{F}_p)$ 是大小为的$n$子群，$H$是$G$的一个陪集 。
- 令$\phi:E \rightarrow E'$是一个次数为2的同源 ，使得$\phi(H)$的大小是$H$的一半。
- 使用$\phi$将$P$分解为两个较小的多项式$P_0,P_1$，使用椭圆曲线$E'$和陪集$\phi(H)$，将 ECFFT 应用于$P_0,P_1$上。
我现在将更详细地解释这些步骤。用 Sage 给出 Secp256k1 曲线的基域的实现, 即$\mathbb{F}_p$
$p = 115792089237316195423570985008687907853269984665640564039457584007908834671663$
p = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F F = GF(p) 寻找曲线
我们需要找到一条曲线$E(\mathbb{F}_p)$使得$n | #E(\mathbb{F}_p)$。论文里给出了这样做的算法。然而，暴力搜索对于较小的值$n$来说也足够有效了。
这是一个通过暴力搜索寻找曲线的脚本
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) 对$n = 2^{12}$，这可以在大约一个小时内找到一条曲线 。这是我找到的曲线$E : y^2 = x^3 + x + 641632526086088189863104799840287255425991771106608433941469413117059106896$
选择陪集
令$G < E$是大小为$n$的子群 。我们选择一个$G$的随机陪集$H$。
gen = E.gens()[0]
G = (gen.order()//n) * gen
R = E.random_element()
H = [R + i*G for i in range(n)] 最后, 令$L$是$H$在$\mathbb{F}_p$上的投影 。$L$将是我们将在其上执行 ECFFT的大小为的$n$集合。
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) 使用$E$的$x$-坐标，$\phi$由2次有理函数$\psi(X) = \mu(X)/\nu(X)$出。
psi = phi.x_rational_map() u, v = psi.numerator(), psi.denominator() 由此我们可以得到一条新的椭圆曲线$E'$, 这是$\phi$的到达域,一个新的陪集$\phi(H)$, 以及一个新的$\psi(L)$给出的$\mathbb{F}_p$的子集。
上分解多项式
在经典 FFT 中, 我们能根据两个较小多项式$P_0,P_1$在较小的集合$H^2$上的求值写出$P$在$H$上的求值 。
我们现在想要通过替换来执行相同的操作，将$H$替换成$L$,$H^2$替换成$\psi(L)$, 以及平方操作替换成$\psi$。
令$P$是一个次数$< n$的多项式 。那么存在次数$< n/2$的多项式 ，使得 $P({X}) = (P{0}(\psi({X})) + {X}P{1}(\psi({X})))v({X})^{n/2-1}$ 证明见论文附录 A。证明的思想是，从$< n/2$次多项式对到次数$< n$多项式的线性映射$(P_0,P_1) \mapsto P$是单射, 同时也是双射, 因为它的定义域和到达域具有与$\mathbb{F}_p$-向量空间同样的维度。
计算多项式$P_0,P_1$不像经典 FFT 那样简单。然而, 可以很容易从$P$在$L$上的求值到$P_0$和$P_1$在$\psi(L)$上的求值，反之亦然。
实际上, 给定$s_0,s_1\in L$使得$\psi(s_0) = \psi(s_1) = t \in \psi(L)$, 通过上面的公式写出$P(X)$在$X = s_0,s_1$处的求值，我们得到以下线性关系（令$q = n/2 - 1$） $\begin{bmatrix} P(s_0) \ P(s_1) \end{bmatrix} = \begin{bmatrix} v(s_0)^q & s_0v(s_0)^q \ v(s_1)^q & s_1v(s_1)^q \end{bmatrix} \begin{bmatrix} P_0(t) \ P_1(t) \end{bmatrix}$ 可以看出矩阵是可逆的, 因此我们可以很容易地从$(P(s_0),P(s_1))$到$(P_0(t),P_1(t))$，反之亦然。这给出了我们稍后将使用的以下有效的一一对应关系：
$P$在$L$上的求值$\longleftrightarrow$$P_0,P_1$在$\psi(L)$上的求值
更重要的是, 矩阵不依赖于$P$, 因此我们可以预先计算它，并将其重用于 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)] 给定次数$< n/2$的多项式$Q$在$S$上的求值，EXTEND 操作计算$Q$在$S'$上的求值。该论文的主要结果是, 有一个$O(n\log n)$的EXTEND算法。请注意, 一个朴素的算法要通过在$S$上拉格朗日插值恢复$Q$的系数 ，然后在$S'$上求值，需要$O(n^2)$。
该算法的工作原理如下。 如果$#S = #S' = 1,Q$是常数且$S$在$S$和$S'$上的求值是一样的。
否则, 从$Q$在$S$上的求值导出$Q_0,Q_1$在$\psi(S)$上的求值, 如上节所示。然后应用 EXTEND操作两次以获得$Q_0,Q_1$在$\psi(S')$上的求值 。最后，恢复$Q$在$S'$上的求值，如上节所示。
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 中的相同, 为$O(n\log n)$。
上求值多项式
上一节中给出的 EXTEND 操作的算法会是论文中给出的许多高效算法的构建块。
其中之一是 ENTER 操作的算法, 该算法计算次数$< n$的多项式$P$在$L$上的求值。这与我们在经典 FFT 中所做的类似。
算法非常简单。按其低系数和高系数分解$P(X)$为$U(X) + X^{n/2}V(X)$。对$U$和$V$在$S$上调用 ENTER 得到$U,V$在$S$上的求值。然后调用 EXTEND 两次以获得$U,V$在$S'$上的求值。因为$L = S\cup S'$, 我们得到$U,V$在$L$上的求值，并且可以导出$P$在$L$上的求值 。
运行时间的递推关系现在为 $F(n) = 2 F(n/2) + O(n \log n)$
因为我们必须在递归步骤中调用 EXTEND。因此运行时间 为$F(n) = O\left( n \log^2 n \right)$，比经典 FFT 稍差。
在 Sage 之外运行
正如这篇文章 (希望) 展示的, 这些算法在像 Sage 这样的计算机代数系统中实现起来非常简单。然而, Sage 是极其慢的, 这些实现远谈不上速度。
这些算法最酷的一点是对于给定的域$\mathbb{F}_p$, 我们可以在 Sage 中预先计算所有必要的数据：所有的集合$L,S,S'$以及 EXTEND 算法中所使用的（逆）矩阵。
一旦我们有了这些预计算数据, 算法就只使用$\mathbb{F}_p$上的基本操作 , 即加法和乘法。
因此, 实际的算法可以很容易地用 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]`
复制链接
Twitter分享
微信扫码分享
