代数快速傅立叶变换

XPTY
发布于 2024-09-14 14:12
阅读 8

代数快速傅立叶变换

>原文：https://rje.li/24-07-03-alg-fft.html
>作者：Eli Riggs
>译者：Kurt Pan

FFT（“快速傅立叶变换”）是有史以来最伟大的算法之一，但大多数人只见过它的经典形式，即将时域中的离散样本转换为频域中的复数表示。它的同类算法 NTT（“数论变换”）对模数上的多项式插值，虽然更少为人所知，但仍在广泛使用，主要用于密码学。

然而，为了更特殊的目的、更陌生的领域和更不熟悉的性质，人们发明了更加奇异的“FFT”。在本文中，我将说明任何一种 FFT 背后的那一个单一算法，然后快速介绍这些奇异的 FFT，并展示如何通过在我们的共用算法中插入一小组参数来实现它们。

最终我们会看到，FFT 让我们可以在$N \log(N)$的时间内**更改一个在$N$个点上定义的函数的基**。朴素的方法是通过应用线性变换矩阵来改基，而矩阵乘法需要 ~$N^3$时间，因此这是一个巨大的优势。

在继续之前，先说几点：
- 我将一组允许我们应用 _FFT 算法_ 的参数称为**一个**FFT 。这起初有点令人困惑，但很快就会变得不那么含糊也更方便。为清楚起见，我仅将复平面上众所周知的 FFT 称为 _经典 FFT_。
- $N$始终表示$2^n$，通常是在指数变得太小的时候。
- 我使用了很多下标，因此向量索引由上标$\overrightarrow{v}^{(i)}$表示。下标$\overrightarrow{v}_n$表示标签或参数。

本文是用 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的定义

对于给定域$\mathbb{F}$和大小$N = 2^n$，一个对$\mathbb{F}$和$N$的 FFT 由其包含$2^n$个点的初始定义域$D_n$，和一个$n$层的序列：
					$\ell_n = [(\pi_n, t_n), (\pi_{n-1}, t_{n-1}), \ldots, (\pi_1, t_1)]$
完全确定。
```
Layer = namedtuple('Layer', 'π t')   
Fft = namedtuple('Fft', 'domain layers')   
```

其中，要么$n = 0$，意味着$D_0$只有一个点，$\ell_0$为空，要么：
- $\pi_n:D_n \rightarrow D_{n-1}$是一个函数，将$D_n$的点映射到一个大小为原来一半的新定义域$D_{n-1}$上，而且它是**二对一**的，意味着对于任意输出点$y \in D_{n-1}$，恰好有两个输入点$x_0,x_1\in D_n$映射到它：$\pi_{n}(x_{0}) = \pi_{n}(x_{1}) = y$。
- $t_n:D_n \rightarrow \mathbb{F}$是一个从$D_n$到$\mathbb{F}$的函数，它对映射到$D_{n-1}$中同一点的$D_n$中的任意两点进行**区分**，或者输出一个不同的值。
- $(D_{n-1},\ell_{n-1})$也是一个 FFT，其中$\ell_n = [(\pi_n, t_n), (\pi_{n-1}, t_{n-1}), \ldots, (\pi_1, t_1)]$，或者用python写为 `layers[1:]`。

根据该定义，我们得到一个由$n + 1$个由映射$\pi _n,\pi _{n-1},...,\pi _1$和除 之外的每个定义域上的区分器（通常称为 _twiddles_）生成的定义域$D _n,D_{n-1},...,D _1$组成的序列，每个定义域的大小减半，
				$D_n \xrightarrow{\pi_n} D_{n-1} \xrightarrow{\pi_{n-1}} \cdots \xrightarrow{\pi_2} D_1 \xrightarrow{\pi_1} D_0$
				 $\downarrow t_n$      $\downarrow t_{n-1}$                   $\downarrow t_1$
				 $\mathbb{F}$          $\mathbb{F}$                          $\mathbb{F}$
显式地构造出逆映射$\pi ^{-1}_n$将会很有用，该映射将$D_{n-1}$中的每个输出点映射到其两个输入$\{x_0, x_1\} \in D_n$上（又名其 _纤维_），以及$t(x_0)$和$t(x_1)$的值。让我们为此向 `Fft` 添加一个方法，它还将检查我们的要求是否得到满足。返回的字典具有以下形式$\{\pi_{n}(x) \to ((x_{0}, t(x_{0})), (x_{1}, t(x_{1})))\}$，其所有键的集合恰好是$D_{n-1}$。
```
@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 算法_ 应用于任意函数$f : D_n \rightarrow \mathbb{F}$来得到一个由$2^n$个系数组成的向量$\vec{c}$：
							$\mathrm{FFT}_{D_{n}, \ell_{n}}(f) = \vec{c}$
给定任意一个 FFT，我们可以导出其 _基_$\hat b_n$，它是$D_n \rightarrow \mathbb{F}$的 _函数_ 向量。然后我们就可以理解系数$\vec{c}$的含义了：它们是基函数的 _权重_，而基的加权和将重构出我们的原始函数$f$。
							$\sum_{i} \overrightarrow{c}^{(i)} \cdot \hat{b}_{n}^{(i)}(x) = f(x)$

嗯……如果你的 FFT 构成的基是具有足够的维度来重构原始函数的，那么它_应该_如此。我们假设这是真的，因为我们只对应用 FFT 算法感兴趣，而不是结果到底有多有用。

我将解释如何导出这个基，因为它是由 FFT 的结构形成的，但需要注意的是，基纯粹是概念性的：它告诉你如何解释得到的系数，但如果你只是想运行算法，它并不必要。

### FFT 算法及其基

给定$f$（$D_n$上的一个函数），FFT 算法有三个步骤：

1. 使用当前层$(\pi_n,t_n)$将$f$**分解**为两个较小的定义域$D_{n-1}$上的函数$f_0$和$f_1$ ，满足：
					$f(x) = f_{0}(\pi_{n}(x)) + t_{n}(x) \cdot f_{1}(\pi_{n}(x))$
2. 在$f_0$和$f_1$上**递归**调用$\mathrm{FFT}_{D_{n-1}, \ell_{n-1}}$以得到它们的系数。
3. 以“合理的方式”**组合$f_0$和$f_1$的系数**，并返回这些组合系数。

在这些步骤中，分解是最容易解释的，但却最难说清楚。现在假设我们会编写一个函数`decompose`，该函数输入逆映射和函数在$D_n$上的求值，返回$f_0$和$f_1$在$D_{n-1}$上的求值。

那么除`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))   
```
我知道这还没什么意义，我只是想把它写出来，因为它非常短。

我们会从头开始，使用上面列出的内容，尝试得出相同的解决方案。

我们的目标是取一个函数$f$并返回一个系数向量$\vec{c}$作为基函数$\hat b_n$的权重，使得加权和可以重构$f$。（在此过程中，我们还需要定义$\hat b_n$）。所以我们需要定义这个等式的左边：
							$\sum_{i} \overline{c}^{(i)} \cdot \hat{b}_{n}^{(i)}(x) = f(x)$
							
在$n = 0$的基础情况下，定义域中只有一个点$x_0$。$f$只是将$x_0$映射到一个单个值，因此可以用一个始终返回该值$f(x_0)$的常数函数来完整描述之。将我们的单个基函数定义为$\hat{b}_0^{(0)}(x) = 1$，可以返回$f(x_0)$作为该函数的系数或权重，从而满足我们的标准：$\sum_i \vec{c}^{(i)} \cdot \hat{b}_n^{(i)}(x) = \vec{c}^{(0)} \cdot \hat{b}_0^{(0)}(x) = f(x_0) \cdot 1 = f(x)$。

否则，继续对$n>0$进行操作。首先，将 `decompose` 函数应用于$f$，将其拆分为$f_0$和$f_1$，或者“依赖于$t(x)$的部分”和“不依赖于$t(x)$的部分”。
					$\sum_i \hat{c}^{(i)} \cdot \hat{b}_n^{(i)}(x) = f_0(\pi_n(x)) + t_n(x) \cdot f_1(\pi_n(x))$

因此，我们从每个$x \in D_n$处的$f(x)$值开始，现在有了每个$x \in D_n$处的$f_0(\pi_n(x))$和$f_1(\pi_n(x))$值，或者等效地说，每个点$y \in D_{n-1}$处的$f_0(y)$和$f_1(y)$值。

现在开始递归，但我们需要小心。我们是将对$f_0$和$f_1$应用下一个（更小的）FFT，**而不是**$f_0 \circ \pi_n$和$f_1 \circ \pi_n$。我的意思是，例如对$f_0$，更小的 FFT 将返回由较小定义域的函数组成较小的基$\hat b_{n-1}$的系数$\vec{c}_0$。为了将其代回我们的方程，必须移动内部的$\pi_n$以将定义域投影到更小的定义域，对于$f_1$也是如此：
		$\sum_i \vec{c}^{(i)} \cdot \hat{b}_n^{(i)}(x) = \left[ \sum_i \vec{c}_0^{(i)} \cdot \hat{b}_{n-1}^{(i)}(\pi_n(x)) \right] + t_n(x) \cdot \left[ \sum_i \vec{c}_1^{(i)} \cdot \hat{b}_{n-1}^{(i)}(\pi_n(x)) \right]$

再次地 - 我们将更小的 FFT 应用于$f_0$和$f_1$，为了使结果准确地重构$f_0(\pi_n(x))$，我们需要将更小的基函数应用于$\pi_n(x)$而不是$x$。

现在我们试着推导$\vec{c}$和$\hat b_n$的表达式，以使得等式始终成立。我们会试着将右侧表示为单个求和，因此首先将$t_n(x)$移入内部：
			$\sum_{i} \vec{c}^{(i)} \cdot \hat{b}_{n}^{(i)}(x) = \sum_{i} \vec{c}_{0}^{(i)} \cdot \hat{b}_{n-1}^{(i)}(\pi_{n}(x)) + \sum_{i} \vec{c}_{1}^{(i)} \cdot t_{n}(x) \cdot \hat{b}_{n-1}^{(i)}(\pi_{n}(x))$

但是……在一般情况下，我们能达到的也就只有这些了。我们有两个不同的内积式的求和，而构造$\vec{c}$和$\hat b_n$以使求和始终相等的唯一方法就是直接连接，设置$\vec{c} = \vec{c}_0||\vec{c}_1$，对于$\hat b_n$也类似，所以我们就是这么做的。

差不多吧。实际上，我们选择交错$\vec{c}_0$和$\vec{c}_1$的元素，但这种选择或任何其他选择都是任意的。我们通常选择交错的原因是$t_n$通常是 1 次函数，而$\pi_n$通常是 2 次函数。因此，当我们交错或将$f_0$中的项放入偶数索引并将$f_1$中的项放入奇数索引时：
			$(\vec{c}^{(i)}, \hat{b}_n^{(i)}(x)) = \begin{cases} (\vec{c}_0^{(k)}, \hat{b}_{n-1}^{(k)}(\pi_n(x))) & i = 2k \\ (\vec{c}_1^{(k)}, t_n(x) \cdot \hat{b}_{n-1}^{(k)}(\pi_n(x))) & i = 2k+1 \end{cases}$`

于是索引将对应于基函数的次数，假设对$\hat b_{n-1}$也是如此。换句话说，假设$\hat b_{n-1}^{(k)}(x)$的次数为$k$，则$\hat b_{n-1}^{(k)}(\pi_n(x))$的次数将为$2k$，而$t_n(x)·\hat b_{n-1}^{(k)}(\pi_n(x))$的次数将为$2k + 1$，与我们的外部索引$i$完全对应。

掌握了以上所有数学知识后，你现在应该可以理解我们算法的最后两行了：
```
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 的每一层，我们将函数$f$（定义在$D_n$上）分解为两个函数$f_0$和$f_1$，每个函数定义在$D_{n-1}$上：
				$f(X) = f_0(\pi_n(X)) + t_n(X) \cdot f_1(\pi_n(X))$

为此我们需要利用$\pi_n$和$t_n$的特殊性质。考虑$D_n$中映射到$D_{n-1}$中的同一点的两个点$x_0$和$x_1$，
				$f(x_0) = f_0(\pi_n(x_0)) + t_n(x_0) \cdot f_1(\pi_n(x_0))$ (1)  
				$f(x_1) = f_0(\pi_n(x_1)) + t_n(x_1) \cdot f_1(\pi_n(x_1))$ (2)

我们知道$f(x_0)$和$f(x_1)$的值，因为$f$（函数或其所有求值的向量）是作为算法的输入给出的，并且我们知道$t_n(x_0)$和$t_n(x_1)$，因为$t_n$是我们 FFT 的一个参数。所以这里有 4 个未知数，这是无解的……但是等等！$x_0$和$x_1$映射到同一点，这意味着$\pi_n(x_0) = \pi_n(x_1)$，我们将其简单地写为$\pi_n(x)$。重写为：
					$f(x_0) = f_0(\pi_n(x)) + t_n(x_0) \cdot f_1(\pi_n(x))$   (3)
					$f(x_1) = f_0(\pi_n(x)) + t_n(x_1) \cdot f_1(\pi_n(x))$   (4)
	
现在我们有 2 个方程和 2 个未知数，可以求解了。如果你仔细看看，你可能会认出这是$y = mx + b$，其中$m = f_1(\pi_n(x))$和$b = f_0(\pi_n(x))$，所以你可能已经知道了如何手动求解之。在我看来，将其视为矩阵乘法会更容易，我们使用 2x2 求逆公式：
					$\begin{bmatrix} f(x_0) \\ f(x_1) \end{bmatrix} = \begin{bmatrix} 1 & t_n(x_0) \\ 1 & t_n(x_1) \end{bmatrix} \begin{bmatrix} f_0(\pi_n(x)) \\ f_1(\pi_n(x)) \end{bmatrix}$                                         (5)
					$\begin{bmatrix} f_0(\pi_n(x)) \\ f_1(\pi_n(x)) \end{bmatrix} = \begin{bmatrix} 1 & t_n(x_0) \\ 1 & t_n(x_1) \end{bmatrix}^{-1} \begin{bmatrix} f(x_0) \\ f(x_1) \end{bmatrix}$                                     (6)
							$= \left( \frac{1}{t_n(x_1) - t_n(x_0)} \right) \begin{bmatrix} t_n(x_1) & -t_n(x_0) \\ -1 & 1 \end{bmatrix} \begin{bmatrix} f(x_0) \\ f(x_1) \end{bmatrix}$           (7)
							$= \left( \frac{1}{t_n(x_1) - t_n(x_0)} \right) \begin{bmatrix} t_n(x_1)f(x_0) - t_n(x_0)f(x_1) \\ f(x_1) - f(x_0) \end{bmatrix}$         (8)
这为我们提供了$f_0(\pi_n(x))$和$f_1(\pi_n(x))$的显式公式，仅和可以访问的$f$和$t_n$相关。你可以看到为什么我们需要$t_n$来区分$\pi_n$的纤维，从而保证$t_n(x_0) \neq t_n(x_1)$，否则我们的矩阵将不可逆，从而导致出现除以零。

或者我们也可以让 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()   
```
							$\begin{pmatrix} f_{x_1}t_{x_0} - f_{x_0}t_{x_1} & f_{x_0} - f_{x_1} \\ t_{x_0} - t_{x_1} & t_{x_0} - t_{x_1} \end{pmatrix}$
注意到$f_1(\pi_n(x)) = \frac{f(x_1) - f(x_0)}{t_n(x_1) - t_n(x_0)}$的表达式对应于计算升降的“斜率”$m$，而$f_0(\pi_n(x))$的表达式对应于“y轴截距”$b$，如果你将其重新带入“$y = mx + b$”，则会得到该值。

如果我们预先计算定义域内所有点处标量的逆 ，就可以通过对的2 次乘法 + 2 *（1 次减法和 1 次乘以逆标量）= 4 次乘法和 2 次减法。

实践中，你会希望选择一个巧妙的twiddle $t_n$和巧妙的定义域来减少这种计算。例如，如果我们知道域中的所有点都有$\left( \frac{1}{t_n(x_1) - t_n(x_0)} \right)$，就像在乘法的情况下一样，那么它的开销会低得多：
`_.substitute(t_x1 = -t_x0).simplify_rational()   `
							$\left( \frac{1}{2} f_{x_0} + \frac{1}{2} f_{x_1}, \frac{f_{x_0} - f_{x_1}}{2 t_{x_0}} \right)$
这可以通过 1 次加法、1 次减法和 2 次乘法来计算。如果将缩放1/2延迟到最后，则可以将其减少为 1 次乘法，再加上一次最终的缩放，最终开销更低（$N$次总缩放乘法，而不是$N \log N$）。

因此，这里是`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))   
```

在每一层，分解都是线性于$N$，然后我们对一半大小的问题进行两次递归。层数为$\log N$，因此总运行时间为$N \log N$。

### 矩阵表示

由于应用 FFT 是在基变换，或函数的线性变换，因此可用应用于函数求值的矩阵来表示。矩阵乘法是$n^3$（量级）的，因此性能要差得多，但我们可以将这种表示用于其他目的，例如计算当用作一个线性码时 FFT 的最小距离。

我们矩阵的列将来自把 FFT 应用于“独热”向量$(0,...,1,...,0)$，可以将其视为将求值的拉格朗日基映射到 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)       
	])   
```

## 例子

每个这些例子都值得写一篇自己的文章，但目前我只想给出一个参考，并展示它们对于其$D$、$\pi$和$t$的选择是有效的。
```
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 在$\mathbb{C}$的乘法子群中起作用，其中$e^{2\pi i/k}$是第$k$个单位根。由于精度问题，这实际上有点难展示。如果我们用浮点数，$\pi$映射就不是精确的二到一的，我们的结果也不会准确。因此，我将展示 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)
```

---

相关文章：  
[FRI邻近性测试](http://mp.weixin.qq.com/s?__biz=MzA4NjU4NzM4Nw==&mid=2247487370&idx=1&sn=56de0ec2e266770d45488c2a26baaf1d&chksm=9fc72d1ea8b0a40837d9ad7dd38fe4bea517482a30df81d5a7a6ba2917308c6515fba54d9e4f&scene=21#wechat_redirect)  
[Kyber/Dilithium中的数论变换](http://mp.weixin.qq.com/s?__biz=MzA4NjU4NzM4Nw==&mid=2247483885&idx=1&sn=02361a43b9a0938ef6c3c0dbdc14b392&chksm=9fc72379a8b0aa6f87ec45e39339444552d484f353ef268f1f1469c610136cdbd81fc90b7aae&scene=21#wechat_redirect)

原文：https://rje.li/24-07-03-alg-fft.html 作者：Eli Riggs 译者：Kurt Pan

最终我们会看到，FFT 让我们可以在$N \log(N)$的时间内更改一个在$N$个点上定义的函数的基。朴素的方法是通过应用线性变换矩阵来改基，而矩阵乘法需要 ~$N^3$时间，因此这是一个巨大的优势。

在继续之前，先说几点：

我将一组允许我们应用 FFT 算法 的参数称为一个FFT 。这起初有点令人困惑，但很快就会变得不那么含糊也更方便。为清楚起见，我仅将复平面上众所周知的 FFT 称为 经典 FFT。
$N$始终表示$2^n$，通常是在指数变得太小的时候。
我使用了很多下标，因此向量索引由上标$\overrightarrow{v}^{(i)}$表示。下标$\overrightarrow{v}_n$表示标签或参数。

本文是用 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的定义

对于给定域$\mathbb{F}$和大小$N = 2^n$，一个对$\mathbb{F}$和$N$的 FFT 由其包含$2^n$个点的初始定义域$D_n$，和一个$n$层的序列： $\ell_n = [(\pi_n, tn), (\pi{n-1}, t_{n-1}), \ldots, (\pi_1, t_1)]$ 完全确定。

Layer = namedtuple('Layer', 'π t')   
Fft = namedtuple('Fft', 'domain layers')

其中，要么$n = 0$，意味着$D_0$只有一个点，$\ell_0$为空，要么：

$\pi_n:Dn \rightarrow D{n-1}$是一个函数，将$Dn$的点映射到一个大小为原来一半的新定义域$D{n-1}$上，而且它是二对一的，意味着对于任意输出点$y \in D_{n-1}$，恰好有两个输入点$x_0,x_1\in Dn$映射到它：$\pi{n}(x{0}) = \pi{n}(x_{1}) = y$。
$t_n:D_n \rightarrow \mathbb{F}$是一个从$Dn$到$\mathbb{F}$的函数，它对映射到$D{n-1}$中同一点的$D_n$中的任意两点进行区分，或者输出一个不同的值。
$(D{n-1},\ell{n-1})$也是一个 FFT，其中$\ell_n = [(\pi_n, tn), (\pi{n-1}, t_{n-1}), \ldots, (\pi_1, t_1)]$，或者用python写为 layers[1:]。

根据该定义，我们得到一个由$n + 1$个由映射$\pi n,\pi {n-1},...,\pi _1$和除之外的每个定义域上的区分器（通常称为 twiddles）生成的定义域$D n,D{n-1},...,D _1$组成的序列，每个定义域的大小减半， $D_n \xrightarrow{\pin} D{n-1} \xrightarrow{\pi_{n-1}} \cdots \xrightarrow{\pi_2} D_1 \xrightarrow{\pi_1} D_0$ $\downarrow tn$ $\downarrow t{n-1}$ $\downarrow t_1$ $\mathbb{F}$ $\mathbb{F}$ $\mathbb{F}$ 显式地构造出逆映射$\pi ^{-1}n$将会很有用，该映射将$D{n-1}$中的每个输出点映射到其两个输入${x_0, x_1} \in D_n$上（又名其纤维），以及$t(x_0)$和$t(x1)$的值。让我们为此向 Fft 添加一个方法，它还将检查我们的要求是否得到满足。返回的字典具有以下形式${\pi{n}(x) \to ((x{0}, t(x{0})), (x{1}, t(x{1})))}$，其所有键的集合恰好是$D_{n-1}$。

@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 算法 应用于任意函数$f : Dn \rightarrow \mathbb{F}$来得到一个由$2^n$个系数组成的向量$\vec{c}$： $\mathrm{FFT}{D{n}, \ell{n}}(f) = \vec{c}$ 给定任意一个 FFT，我们可以导出其基$\hat b_n$，它是$D_n \rightarrow \mathbb{F}$的函数向量。然后我们就可以理解系数$\vec{c}$的含义了：它们是基函数的权重，而基的加权和将重构出我们的原始函数$f$。 $\sum{i} \overrightarrow{c}^{(i)} \cdot \hat{b}{n}^{(i)}(x) = f(x)$

FFT 算法及其基

给定$f$（$D_n$上的一个函数），FFT 算法有三个步骤：

使用当前层$(\pi_n,tn)$将$f$分解为两个较小的定义域$D{n-1}$上的函数$f_0$和$f1$ ，满足： $f(x) = f{0}(\pi{n}(x)) + t{n}(x) \cdot f{1}(\pi{n}(x))$
在$f_0$和$f1$上递归调用$\mathrm{FFT}{D{n-1}, \ell{n-1}}$以得到它们的系数。
以“合理的方式”组合$f_0$和$f_1$的系数，并返回这些组合系数。

在这些步骤中，分解是最容易解释的，但却最难说清楚。现在假设我们会编写一个函数decompose，该函数输入逆映射和函数在$D_n$上的求值，返回$f_0$和$f1$在$D{n-1}$上的求值。

那么除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))

我知道这还没什么意义，我只是想把它写出来，因为它非常短。

我们会从头开始，使用上面列出的内容，尝试得出相同的解决方案。

我们的目标是取一个函数$f$并返回一个系数向量$\vec{c}$作为基函数$\hat b_n$的权重，使得加权和可以重构$f$。（在此过程中，我们还需要定义$\hat bn$）。所以我们需要定义这个等式的左边： $\sum{i} \overline{c}^{(i)} \cdot \hat{b}_{n}^{(i)}(x) = f(x)$

在$n = 0$的基础情况下，定义域中只有一个点$x_0$。$f$只是将$x_0$映射到一个单个值，因此可以用一个始终返回该值$f(x_0)$的常数函数来完整描述之。将我们的单个基函数定义为$\hat{b}_0^{(0)}(x) = 1$，可以返回$f(x_0)$作为该函数的系数或权重，从而满足我们的标准：$\sum_i \vec{c}^{(i)} \cdot \hat{b}_n^{(i)}(x) = \vec{c}^{(0)} \cdot \hat{b}_0^{(0)}(x) = f(x_0) \cdot 1 = f(x)$。

否则，继续对$n>0$进行操作。首先，将 decompose 函数应用于$f$，将其拆分为$f_0$和$f_1$，或者“依赖于$t(x)$的部分”和“不依赖于$t(x)$的部分”。 $\sum_i \hat{c}^{(i)} \cdot \hat{b}_n^{(i)}(x) = f_0(\pi_n(x)) + t_n(x) \cdot f_1(\pi_n(x))$

因此，我们从每个$x \in D_n$处的$f(x)$值开始，现在有了每个$x \in D_n$处的$f_0(\pi_n(x))$和$f_1(\pin(x))$值，或者等效地说，每个点$y \in D{n-1}$处的$f_0(y)$和$f_1(y)$值。

现在开始递归，但我们需要小心。我们是将对$f_0$和$f_1$应用下一个（更小的）FFT，而不是$f_0 \circ \pi_n$和$f_1 \circ \pi_n$。我的意思是，例如对$f0$，更小的 FFT 将返回由较小定义域的函数组成较小的基$\hat b{n-1}$的系数$\vec{c}_0$。为了将其代回我们的方程，必须移动内部的$\pi_n$以将定义域投影到更小的定义域，对于$f_1$也是如此： $\sum_i \vec{c}^{(i)} \cdot \hat{b}_n^{(i)}(x) = \left[ \sum_i \vec{c}0^{(i)} \cdot \hat{b}{n-1}^{(i)}(\pi_n(x)) \right] + t_n(x) \cdot \left[ \sum_i \vec{c}1^{(i)} \cdot \hat{b}{n-1}^{(i)}(\pi_n(x)) \right]$

再次地 - 我们将更小的 FFT 应用于$f_0$和$f_1$，为了使结果准确地重构$f_0(\pi_n(x))$，我们需要将更小的基函数应用于$\pi_n(x)$而不是$x$。

现在我们试着推导$\vec{c}$和$\hat b_n$的表达式，以使得等式始终成立。我们会试着将右侧表示为单个求和，因此首先将$tn(x)$移入内部： $\sum{i} \vec{c}^{(i)} \cdot \hat{b}{n}^{(i)}(x) = \sum{i} \vec{c}{0}^{(i)} \cdot \hat{b}{n-1}^{(i)}(\pi{n}(x)) + \sum{i} \vec{c}{1}^{(i)} \cdot t{n}(x) \cdot \hat{b}{n-1}^{(i)}(\pi{n}(x))$

差不多吧。实际上，我们选择交错$\vec{c}_0$和$\vec{c}_1$的元素，但这种选择或任何其他选择都是任意的。我们通常选择交错的原因是$t_n$通常是 1 次函数，而$\pi_n$通常是 2 次函数。因此，当我们交错或将$f_0$中的项放入偶数索引并将$f_1$中的项放入奇数索引时： $(\vec{c}^{(i)}, \hat{b}_n^{(i)}(x)) = \begin{cases} (\vec{c}0^{(k)}, \hat{b}{n-1}^{(k)}(\pi_n(x))) & i = 2k \ (\vec{c}_1^{(k)}, tn(x) \cdot \hat{b}{n-1}^{(k)}(\pi_n(x))) & i = 2k+1 \end{cases}$`

于是索引将对应于基函数的次数，假设对$\hat b{n-1}$也是如此。换句话说，假设$\hat b{n-1}^{(k)}(x)$的次数为$k$，则$\hat b_{n-1}^{(k)}(\pi_n(x))$的次数将为$2k$，而$tn(x)·\hat b{n-1}^{(k)}(\pi_n(x))$的次数将为$2k + 1$，与我们的外部索引$i$完全对应。

掌握了以上所有数学知识后，你现在应该可以理解我们算法的最后两行了：

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 的每一层，我们将函数$f$（定义在$D_n$上）分解为两个函数$f_0$和$f1$，每个函数定义在$D{n-1}$上： $f(X) = f_0(\pi_n(X)) + t_n(X) \cdot f_1(\pi_n(X))$

为此我们需要利用$\pi_n$和$t_n$的特殊性质。考虑$Dn$中映射到$D{n-1}$中的同一点的两个点$x_0$和$x_1$， $f(x_0) = f_0(\pi_n(x_0)) + t_n(x_0) \cdot f_1(\pi_n(x_0))$ (1)
$f(x_1) = f_0(\pi_n(x_1)) + t_n(x_1) \cdot f_1(\pi_n(x_1))$ (2)

现在我们有 2 个方程和 2 个未知数，可以求解了。如果你仔细看看，你可能会认出这是$y = mx + b$，其中$m = f_1(\pi_n(x))$和$b = f_0(\pi_n(x))$，所以你可能已经知道了如何手动求解之。在我看来，将其视为矩阵乘法会更容易，我们使用 2x2 求逆公式： $\begin{bmatrix} f(x_0) \ f(x_1) \end{bmatrix} = \begin{bmatrix} 1 & t_n(x_0) \ 1 & t_n(x_1) \end{bmatrix} \begin{bmatrix} f_0(\pi_n(x)) \ f_1(\pi_n(x)) \end{bmatrix}$ (5) $\begin{bmatrix} f_0(\pi_n(x)) \ f_1(\pi_n(x)) \end{bmatrix} = \begin{bmatrix} 1 & t_n(x_0) \ 1 & t_n(x_1) \end{bmatrix}^{-1} \begin{bmatrix} f(x_0) \ f(x_1) \end{bmatrix}$ (6) $= \left( \frac{1}{t_n(x_1) - t_n(x_0)} \right) \begin{bmatrix} t_n(x_1) & -t_n(x_0) \ -1 & 1 \end{bmatrix} \begin{bmatrix} f(x_0) \ f(x_1) \end{bmatrix}$ (7) $= \left( \frac{1}{t_n(x_1) - t_n(x_0)} \right) \begin{bmatrix} t_n(x_1)f(x_0) - t_n(x_0)f(x_1) \ f(x_1) - f(x_0) \end{bmatrix}$ (8) 这为我们提供了$f_0(\pi_n(x))$和$f_1(\pi_n(x))$的显式公式，仅和可以访问的$f$和$t_n$相关。你可以看到为什么我们需要$t_n$来区分$\pi_n$的纤维，从而保证$t_n(x_0) \neq t_n(x_1)$，否则我们的矩阵将不可逆，从而导致出现除以零。

或者我们也可以让 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()

                        $\begin{pmatrix} f_{x_1}t_{x_0} - f_{x_0}t_{x_1} & f_{x_0} - f_{x_1} \\ t_{x_0} - t_{x_1} & t_{x_0} - t_{x_1} \end{pmatrix}$

注意到$f_1(\pi_n(x)) = \frac{f(x_1) - f(x_0)}{t_n(x_1) - t_n(x_0)}$的表达式对应于计算升降的“斜率”$m$，而$f_0(\pi_n(x))$的表达式对应于“y轴截距”$b$，如果你将其重新带入“$y = mx + b$”，则会得到该值。

如果我们预先计算定义域内所有点处标量的逆，就可以通过对的2 次乘法 + 2 *（1 次减法和 1 次乘以逆标量）= 4 次乘法和 2 次减法。

实践中，你会希望选择一个巧妙的twiddle $t_n$和巧妙的定义域来减少这种计算。例如，如果我们知道域中的所有点都有$\left( \frac{1}{t_n(x_1) - t_n(x0)} \right)$，就像在乘法的情况下一样，那么它的开销会低得多： `.substitute(t_x1 = -t_x0).simplifyrational() ` $\left( \frac{1}{2} f{x0} + \frac{1}{2} f{x1}, \frac{f{x0} - f{x1}}{2 t{x_0}} \right)$ 这可以通过 1 次加法、1 次减法和 2 次乘法来计算。如果将缩放1/2延迟到最后，则可以将其减少为 1 次乘法，再加上一次最终的缩放，最终开销更低（$N$次总缩放乘法，而不是$N \log N$）。

因此，这里是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))

在每一层，分解都是线性于$N$，然后我们对一半大小的问题进行两次递归。层数为$\log N$，因此总运行时间为$N \log N$。

矩阵表示

我们矩阵的列将来自把 FFT 应用于“独热”向量$(0,...,1,...,0)$，可以将其视为将求值的拉格朗日基映射到 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)       
    ])

例子

每个这些例子都值得写一篇自己的文章，但目前我只想给出一个参考，并展示它们对于其$D$、$\pi$和$t$的选择是有效的。

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

注意这里的基正是你所期望的那些。

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.&lt;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中的参数。

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)

相关文章：
FRI邻近性测试
 Kyber/Dilithium中的数论变换