Plonky3 Mersenne素数域的Reed-Solomon codes设计

1. 引言

主要见2023年Polygon Labs和Polygon Zero团队论文《Reed-Solomon codes over the circle group》。
相关开源代码实现见：

• https://github.com/Plonky3/Plonky3/blob/main/mersenne-31/src/lib.rs（Rust）

本文重点讨论基于Mersenne素数域$F=2^{31}-1$的复数扩域$\mathbb{C}(F)$的unit circle内所定义的Reed-Solomon codes：

• 1）在该unit circle内，interpolants（插值）是"real"，即$F$-valued，函数也几乎是real的，意味着几乎无额外开销的情况下，可将函数值矫正为某实数值表示。
• 2）采用实数值（real-valued）函数的标准FFT技术，可有效加速编码。

由于Mersenne域的特别高效运算，可期待对其实现“almost native” Reed-Solomon codes，其原生性与 基于高two-adicity的（但不够处理器友好运算的）素数域的Reed-Solomon codes 基本相当。

STARKs（scalable and transparent arguments of knowledge）基于：

• Reed-Solomon codes
• 和 FRI proof of proximity

STARKs的性能取决于其witness trace所表示用的有限域。StarkWare的Cairo zkVM仍使用的是密码学大素数域，Plonky2采用的是64位小素数域，其模为：
$p=2^{64}-2^{32}+1$
通常称其为Goldilocks域。

Goldilocks域的模支持在标准计算机架构中高效实现，其乘法群足够smooth，two-adic subgroup size大至$2^{32}$，可支持高效FFT Reed-Solomon编码。目前Polygon Hermez和Polygon Miden等项目均采用了Goldilocks域。

但是，从高效域实现以及算术电路的角度来看，需要更小的域。本文成稿时，RISC Zero为唯一往更小域努力的项目，其采用31位Baby Bear小素数域：
$p=2^{31}-2^{27}+1$
作为其witness的原生模。与Goldilocks域一样，Baby Bear的multiplicative group对大多数实际应用来说足够smooth：$p-1=2^{27}\cdot 3\cdot 5$。与Baby Bear类似的现有素数，如$p=2^{31}+2^{30}+1$，其也有相当的smooth multiplicative group和高效域运算。

最高效的域运算似乎是Mersenne filed，其形式为$p=2^e-1$。
特别是Mersenne素数域$p=2^{31}-1$，其可在32位架构中非常高效地运行：

• 由于$2^{32}=2(\mathrm{mod}p)$，因此形如$2^{32}\cdot x_{hi}+x_{lo}$的宽乘积编码，可 reduce为小得多的量级：$2\cdot x_{hi}+x_{lo}$。
• 但是，由于$p-1=2\cdot 3^2\cdot 7\cdot 11\cdot 31\cdot 151\cdot 331$，其乘法群缺少two-adic subgroups。而two-adic subgroups对于高效Cooley-Tukey FFT实现非常有用。

为此，本文描述了对Mersenne field的almost native Reed-Solomon codes，其利用的是基于复数扩域的circle group $S_1$的FFT编码。

扩展到复数扩域，初看似乎效率不高：

• 1）扩域运算要比基域运算昂贵得多。
• 2）code words以复数表示，会将code rate减半。

但是，对于定义于$S_1$ subgroup的实数函数（即basefield-valued），其基于该circle group的剩余部分的插值点，也几乎是实数的，即意味着其值在某线性实数子空间内，其仅依赖于该subgroup的coset。进一步，采用标准技术来改进real-valued函数的FFT，插值开销将减半。
基于这些事实，所获得的Reed-Solomon codes是“almost native”的：

• 1）对于实数函数，编码开销，与具有相同size且high two-adicity的Cooley-Tukey-based编码开销，基本相当。
• 2）所编码获得的codewords，可压缩为完全实数表示，从而具有与原生方式完全相同的承诺开销。

由于$p=2^{31}-1$具有特别高效的运算，其almost native Reed-Solomon codes性能，与Baby Bear的原生解决方案，基本相当。同时，$p=2^{31}-1$具有更快的域运算，可用于trace计算和算术哈希计算。

2. 复数扩域和实数FFT

尽管本文重点关注Mersenne域，但本节结论适于每个素数域$F=F_p$，其中$p=3\mathrm{mod}4$。即要求$\frac{p-1}{2}$是技术，或$-1$不是a quadratic residue。此时：

• 多项式$X^2+1$为irreducible over $F$

从而可构建复数扩域$\mathbb{C}(F)=F[X]/(X^2+1)$。该复数扩域通过连接formal root $i:=\sqrt{-1}$而构成：
C ( F ) = { x + i ⋅ y : x , y ∈ F } \mathbb{C}(F)=\{x+i\cdot y:x,y\in F\} C(F)={x+i⋅y:x,y∈F}
该复数扩域的域运算被强加了对该root的代数约束$i^2=-1$。

基于$F$的unit circle为代数集合：$S_1=\{(x,y)\in F^2:x^2+y^2=1\}$，或以复数形式表示为：
S 1 = { z ∈ C ( F ) ∗ : z ⋅ z ˉ = 1 } S_1=\{z\in\mathbb{C}(F)^*:z\cdot \bar{z}=1\} S1​={z∈C(F)∗:z⋅zˉ=1}
其中$\bar{z}$为$z=x+i\cdot y$的共轭复数$\bar{z}=x-i\cdot y$。由于共轭复数是field automorphism的，因此$S_1$在复数乘法下是闭环的，形成复数乘法群$\mathbb{C}(F{)}^{\ast }$的子群——即（unit）circle group。

Lemma 1 (Circle group):

• 令$F=F_p$为素数域，有$p=3\mathrm{mod}4$，则基于$F$的circle group $S_1$，为$\mathbb{C}(F)$中的$(p+1)$-th roots of unity group，其order为$p+1$。

由于有共轭等价为Frobenius isomorphism，即$\bar{z}=z^p$。从而有$x^2+y^2=z\cdot \bar{z}=z^{p+1}$，同时有$(p+1)|(p^2-1)$，从而可得出结论：该circle group为$(p+1)$-th roots of unity subgroup，具有的order为$p+1$。

基于$F$的几何复数平面，类似于calculus经典案例。尽管后续不需要，仍可迅速总结出复数值的极坐标表示。

Proposition 1：

• 令$F=F_p$为素数域，有$p=3\mathrm{mod}4$，由所有复数平方 Q C = { z 2 : z ∈ C ( F ) ∗ } Q_{\mathbb{C}}=\{z^2:z\in\mathbb{C}(F)^*\} QC​={z2:z∈C(F)∗}组成的乘法子群，可分解为：
  $Q_{\mathbb{C}}=Q_F\otimes S_1$
  其中：
  • $Q_F$为基域$F$内的$\frac{p-1}{2}$个quadratic residues，
  • $S_1$为基于$F$的order为$p+1$的unit circle group。

Remark 1：

• 从group orders来看，有$|Q_{\mathbb{C}}|=\frac{p^2-1}{2}=\frac{(p-1)\cdot (p+1)}{2}$。实数quadratic residues $Q_F$对应奇数$\frac{p-1}{2}$项，且unit circle group $S_1$互素数因子$p+1$。

Remark 2：

• 整个乘法群$\mathbb{C}(F{)}^{\ast }$可分解为：
  $\mathbb{C}(F{)}^{\ast }=Q_F\otimes S_{\pm 1}$
  其中：
  • S ± 1 = { z ∈ C ( F ) ∗ : z ⋅ z ˉ ∈ { ± 1 } } S_{\pm 1}=\{z\in\mathbb{C}(F)^*: z\cdot \bar{z}\in\{\pm 1\}\} S±1​={z∈C(F)∗:z⋅zˉ∈{±1}}，为由unit circle group和anti-unit circle group S − 1 = { z ∈ C ( F ) ∗ : z ⋅ z ˉ = − 1 } S_{-1}=\{z\in\mathbb{C}(F)^*: z\cdot \bar{z}=-1\} S−1​={z∈C(F)∗:z⋅zˉ=−1}组成的subgroup。其中$S_{\pm 1}$对应的group order为$|\mathbb{C}(F{)}^{\ast }|=2\cdot (p+1)\cdot \frac{p-1}{2}$中的$2\cdot (p+1)$项。

2.1 傅里叶变换和实数FFT

令$H$为任意有限域$K$的乘法子群，其order为$|H|=N$，令$g$为$H$的generator。
对函数$f:H\to K$的（discrete）Fourier transform $\hat{f}$函数基于$H$定义为：
$$\begin{array}{c}\hat{f}(g^k):=\frac{1}{N}\cdot \sum _{j=0}^{N-1}f(g^j)\cdot g^{-k\cdot j}\end{array}$$
其中$k=0,\cdots ,N-1$。这些值是$f$插值点的系数，即，存在唯一的degree小于$N$的多项式$p(X)$，$p(X)$是对$f$基于$H$对所有的$i=0,\cdots ,N-1$插值：
$$\begin{array}{c}f(g^i)=\sum _{k=0}^{N-1}\hat{f}(g^k)\cdot g^{k\cdot i}\end{array}$$

等式（2）右侧为$\hat{f}:H\to K$的inverse Fourier transform。

用于傅里叶变换的标准算法为Cooley-Tukey算法。该算法利用了$H$的group结构，并利用所有的：
$$\begin{array}{c}|H|\cdot \log |H|\cdot (\frac{1}{2}\cdot M+A)\end{array}$$
来计算上面的等式（1）（以及等式（2）），其中：【详情见Introduction to algorithms (3rd ed.)，或 Modern computer algebra (3rd ed.)】

• $M$：表示$K$中的乘法运算。
• $A$：表示$K$中的加法/减法运算。

当$K=\mathbb{C}(F)$时，即为有限素数域$F$（模为$p=3\mathrm{mod}4$）的复数扩域，可使用order为$8$的如下roots of unity简单形式：
$$c\cdot (1+i),i,c\cdot (-1+i),-1,c\cdot (-1-i),-i,c\cdot (1-i),1$$
其中：

• $c$为$+\frac{1}{2}$或$-\frac{1}{2}$的实数二次方根，具体取决于二者哪个是quadratic residue modulo $p$。
• 这些roots的乘法，包含更少的实数乘法，甚至没有实数乘法。从而产生了更高效的FFT更高radix变种，如：
  • radix-4和radix-8变种。（详情见 A fast Fourier transform algorithm using base 8 iterations.）
  • split-radix变种：
    • 开销要略少于radix-8，在3/3-regme中，最多为：$|H|\cdot \log |H|\cdot (1\cdot M_F+3\cdot A_F)$，其中每个复数乘法运算对应3个实数乘法运算$M_F$和3个实数加法运算$A_F$（详情见Real-valued fast Fourier transform algorithms）。
    • 比Cooley-Tukey radix-2 FFT算法少33%的乘法运算和少20%的加法运算。
  • radix 2和radix 4混合变种。

此外，所有这些不同radix FFT算法，都可利用傅里叶变换的Hermitian对称性来对real-valued函数$f:H\to F$进行加速，即对每个$w\in H$，有：
$\hat{f}(\bar{w})=\overline{\hat{f}(w)}$
该对称性可用于让算法中每一步的计算次数减半，对反傅里叶变换也可使用Hermitian对称性。

尽管本文描述的是对实数$\mathbb{R}$的复数扩域的经典情况，该算法及其开销分析也适用于有限域。具体总结为Theorem 3（详情见Real-valued fast Fourier transform algorithms）：

• 令$K=\mathbb{C}(F)$对有限素数域$F$的复数扩域，其中$F$域的素数模满足$p=3\mathrm{mod}4$，令$H$为$K$的two-adic乘法子群。对real-valued函数$f:H\to F$的傅里叶变换，以及对共轭Hermitian函数$\hat{f}:H\to K$的反傅里叶变换，均可以少于：
  $$\begin{array}{c}|H|\cdot \log |H|\cdot (\frac{1}{2}\cdot M+\frac{3}{2}\cdot A)\end{array}$$
  次数计算，其中$M$和$A$分别表示基于基域$F$的域乘法和域加法运算。

Remark 4：

• 当针对模为$p=2^{31}-1$的Mersenne域，其8-th primitive root of unity为$\frac{1}{\sqrt{2}}\cdot (1+i)=2^{15}\cdot (1+i)$。使得与8-th roots of unity的乘法运算仅包含位移、加法和减法。为此，可使用packing trick和radix-8 transform来实现类似的运算次数：
  • 1）将函数按size对半切分（对应2个cosets of an index 2 subgroup值）
  • 2）将他们组合为单个复数形式，并应用radix-8算法。
  • 3）可从步骤2）中的结果中提取2个half-sized函数的傅里叶变换，然后按Cooley-Tukey方式实现所需结果。

3. 基于$S_1$插值

如下proposition是基于unit circle group的Reed Solomon code实数表示的关键。尽管本文关注的是Mersenne素数域，但结论适用于所有满足$p=3\mathrm{mod}4$的素数域。

Proposition 2：

• 令$F=F_p$为素数域，其中$p=3\mathrm{mod}4$，且$H$为复数扩域$\mathbb{C}(F)$的unit circle group $S_1$的子群，$H$的order为偶数 $|H|\ge 2$。对于每个基于$H$的$F$-valued函数$f$，其基于coset $\tau \cdot H,\tau \in S_1$的插值$p(X)=c_0+{\sum }_{k=1}^{|H|-1}{c}_k\cdot X^k$的值，在由$c_0和\tau$决定的$\mathbb{C}(F)$实数线性子空间中，即，对于所有的$x\in \tau \cdot H$，有：
  $$\begin{array}{c}p(x)-c_0\in \varphi (\tau )\cdot F\end{array}$$
  其中：
  $$\begin{array}{c}\varphi (\tau )={\tau }^{\frac{|H|}{2}}\end{array}$$
  此外，mapping $\varphi$定义了injective homomorphism ϕ : S 1 / H → S 1 / { ± 1 } \phi:S_1/H\rightarrow S_1/\{\pm 1\} ϕ:S1​/H→S1​/{±1}，其image等于the “projective” cyclic subgroup C 2 ⋅ β / { ± 1 } C_{2\cdot \beta}/\{\pm 1\} C2⋅β​/{±1}，其中$\beta =|S_1|/|H|$。

Remark 5：

• 对于具有奇数order的子群$H$，以上结论也成立，只是需更小心的证明，因为$\sqrt{\tau }$可能在anti-unit circle $S_{-1}$中，从而破坏了conjugate-equals-inverse关系$\sqrt{\tau }\cdot \overline{\sqrt{\tau }}=-1$。本文不使用该特例情况，因此忽略相应细节。

由于$f$是基于$H$的实数函数，其插值$p(X)={\sum }_{k=0}^{|H|-1}{c}_k\cdot X^k$的系数对其共轭是Hermitian对称的，即$c_0={⟨c|}_0$，对于所有的$k=1,\cdots ,|H|-1$，有$c_{|H|-k}={\bar{c}}_k$。
接下来考虑基于$\tau \cdot H$的$p(X)-c_0$值，其中$\tau \in S_1$，有：
$p(\tau \cdot X)-c_0={\sum }_{k=1}^{|H|-1}{c}_k\cdot {\tau }^k\cdot X^k={\sum }_k{d}_k\cdot X^k$
其中：
$d_k=\{\begin{array}{cc}{c}_k\cdot {\tau }^k& \text{for }1\le k\le |H|-1\\ 0& \text{otherwise.}\end{array}$

将该函数扩大${\tau }^{-\frac{|H|}{2}}$，有系数$d_k^{\prime }:={\tau }^{-\frac{|H|}{2}\cdot d_k}$，其满足Hermitian对称性，即对所有的$k=1,\cdots ,|H|-1$，有：
$d_{|H|-k}^{\prime }={\tau }^{-\frac{|H|}{2}}\cdot {\tau }^{|H|-k}\cdot c_{|H|-k}={\tau }^{\frac{|H|}{2}}\cdot {\tau }^{-k}\cdot {\bar{c}}_k=\overline{{\tau }^{-\frac{|H|}{2}}\cdot {\tau }^k\cdot c_k}={\overline{d^{\prime }}}_k$

从而可得出结论：

• 对于每个$x\in H$，可得出结论${\tau }^{-\frac{|H|}{2}}\cdot (p(\tau \cdot x)-c_0)\in F$，从而证明了Proposition 2中的第一条claim。

对于Proposition 2中的第二条claim：

• 事实上，对于每个$x\in H$，有 x ∣ H ∣ 2 ∈ { ± 1 } x^{\frac{|H|}{2}}\in\{\pm 1\} x2∣H∣​∈{±1}，从而有 ( τ ⋅ x ) ∣ H ∣ 2 ∈ τ ∣ H ∣ 2 ⋅ { ± 1 } (\tau \cdot x )^{\frac{|H|}{2}}\in \tau ^{\frac{|H|}{2}}\cdot \{\pm 1\} (τ⋅x)2∣H∣​∈τ2∣H∣​⋅{±1}。由于$\varphi$指向由$S_1/H$到 S 1 / { ± 1 } S_1/\{\pm 1\} S1​/{±1}的isomorphism，该isomorphism的kernel是trivial的，因为 τ ∣ H ∣ 2 ∈ { ± 1 } \tau^{\frac{|H|}{2}}\in\{\pm 1\} τ2∣H∣​∈{±1}等价为${\tau }^{|H|}=1$。由于$S_1$是cyclic的，因此，the projective unit circle S 1 / { ± 1 } S_1/\{\pm 1\} S1​/{±1}和$\varphi$的images order为$\beta =|S_1/H|$的唯一子群，其等价为 C 2 ⋅ β / { ± 1 } C_{2\cdot \beta}/\{\pm 1\} C2⋅β​/{±1}子群。

4. Mersenne域的almost native Reed-Solomon codes

对Mersenne域$F=F_p$，其中素数模满足$p=2^e-1$，应用如上结论。

根据Lemma 1，unit circle group $S_1$是$\mathbb{C}(F{)}^{\ast }$的purely two-adic子群：
$|S_1|=p+1=2^e$

其可用于第2节中的FFT算法。在STARKs上下文中，使用witness多项式的Lagrange表示，如AIR（ algebraic intermediate representations）或Plonkish算术化，这就需要面对如下问题：

• 已知witness多项式基于某witness domain $H$的值，对其插值到更大的sampling domain $D$，有$|D|=\beta \cdot |H|$，其中$\beta$为整数，并对所获得的值通过Merkle hash进行承诺。若$H$和$D$均在unit circle group $S_1$内，则可获得almost native Reed-Solomon code：
  • 对real-valued函数做mixed-radix FFT（见Theorem 3），计算插值及其所需的coset值，其开销与等同size具有high two-adicity域 的开销几乎相同。
    • 首先计算real witness函数$w:H\to F$的傅里叶变换，然后，对每个coset KaTeX parse error: Undefined control sequence: \codt at position 5: \tau\̲c̲o̲d̲t̲ ̲H\subset D，将shifted factors ${\tau }^{k-\frac{|H|}{2}},k=1,\cdots ,|H|-1$与$w$的傅里叶变换相乘。
    • 然后对Hermititian对称函数应用反傅里叶变换，以获得矫正后的基于$\tau \cdot H$的${\tau }^{-\frac{|H|}{2}}\cdot (w(X)-c_0)$的实数值。
  • codeword的native表示中，包含实数coset evaluations：
    $$\begin{array}{c}{\tau }^{-\frac{|H|}{2}}\cdot (w(x)-c_0{)}_{x\in \tau \cdot H}\end{array}$$
    其中$\tau \cdot H\subset D$，常量项$c_0={\sum }_{x\in H}w(x)$。
    • 不对$w(x){|}_{x\in D}$进行承诺，改为仅对（7）中的校正值进行承诺，并额外对scalar $c_0$进行announce或commit。
    • 在很多应用中，甚至可假设$c_0=0$，通常，并不是整个domain $H$都会被witness数据所消费，可将未用值调整使该domain之和为0。

在非零值设置下，当考虑sampling domain $D=S_1$的size为$2^{31}$，对应的blow-up factor $\beta =2$时，对Mersenne素数$p=2^{31}-1$可支持size到$2^{30}$的witness domain $H$的almost native Reed-Solomon codes。

当考虑零知识属性时，$H$的最大size 为$2^{29}$，采用的disjoint sampling domain $D$ within in $S_1$ of double the size of $H$。

在标准计算机架构下，Mersenne算术运算可balance掉，基于复数扩域的mixed-radix FFT所引入的加法运算此时。具体benchmark为：

对应上面的等式（3）和（4），基于M31（Mersenne31域）的mixed-radix real FFT，在Apple M1 ARM处理器上和Intel Ice Lake x86处理器上的开销分别为：
$|H|\cdot \log |H|\cdot (\frac{1}{2}\cdot 0.31+\frac{3}{2}\cdot 0.19)\approx 0.44\cdot |H|\cdot \log |H|$
$|H|\cdot \log |H|\cdot (\frac{1}{2}\cdot 0.31+\frac{3}{2}\cdot 0.09)\approx 0.31\cdot |H|\cdot \log |H|$
个clock cycles，而对Baby Bear域应用Cooley-Tukey算法，在Apple M1 ARM处理器上和在Intel Ice Lake x86处理器上开销分别为：
$|H|\cdot \log |H|\cdot (\frac{1}{2}\cdot 0.5+\frac{1}{2}\cdot 0.19)\approx 0.44\cdot |H|\cdot \log |H|$
$|H|\cdot \log |H|\cdot (\frac{1}{2}\cdot 0.44+\frac{1}{2}\cdot 0.09)\approx 0.31\cdot |H|\cdot \log |H|$

因此，基于Mersenne域的almost native Reed-Solomon codes，其性能与具有相同size和高two-adic域 的性能相当，同时，提供了额外更快的域运算。