Polygon zkEVM PIL-STARK Spearbit审计报告解读（2023年6月版）

mutourend

已于 2024-01-11 17:12:04 修改

阅读量1k

点赞数 23

分类专栏： zkVM 文章标签： zkVM

于 2024-01-10 16:57:35 首次发布

本文链接：https://blog.csdn.net/mutourend/article/details/135507239

版权

zkVM 专栏收录该内容

130 篇文章 22 订阅

订阅专栏

1. 引言

前序博客有：

具体内容见：

Polygon zkEVM Cryptography Security Review (Pil-stark and Snarkjs cryptography review), June 20, 2023

本轮审计历时2天，Spearbit团队重点审计了Polygon zkEVM proving pipeline的2大部分：

1）审计了在该pipeline的final stage，对fflonk verifier合约的优化修改。
2）审计了STARK verifier中的一些定制PIL约束，包括：
- Fp3 multiplication
- 1-of-4 Fp3 selector
- Poseidon hash evaluation
这些约束用于proving pipeline的recursive composition stage中。

实际审计的代码库为：

本轮审计未发现可靠性问题：

有一些对fflonk verifier的进一步优化建议。
对定制PIL约束的一些小的文档/优化建议。

2. fflonk verifier优化

对https://github.com/iden3/snarkjs/blob/e70271b30a6392a34f664e6bececddad284c8854/test/fflonk/verifier.sol FFlonk verifier合约进行审计和优化。

定义了乘法群、其order、相关常量以及cosets：

$H =< w >, H = 256$
$H_4=<w_3>,H_3=3$
$H_4=<w_4>,H_4=4$
$H_8=<w_8>,H_8=8$
$r\leftarrow \mathbb{F},\xi\leftarrow r^{24}$
$h_0=r^3,h_1=r^6,h_2=r^8,h_3=r^8w^{1/3}$
$h_0^8=\xi,h_1^4=\xi,h_2^3=\xi,h_3^3=\xi$
$S_0=h_0H_8=r^3H_8$
$S_1=h_1H_4=r^6H_4$
$S_2=h_2H_3\cup h_3H_3=r^8H_3\cup r^8w^{1/3}H_3$

https://github.com/iden3/snarkjs/blob/e70271b30a6392a34f664e6bececddad284c8854/test/fflonk/verifier.sol FFlonk verifier合约内有：

常量值：
- $w,w^{1/3}$
- ${w_3,w_3^2\}$
- ${w_4,w_4^2,w_4^3\}$
- $\{w_8,w_8^2,\cdots,w_8^7\}$
fflonk_verify.js中的computeChallenges，会计算集合 $S_0,S_1,S_2$ 中的每个元素。
fflonk_prove.js中的computeLiS0，为对langrange多项式： $L_i^{(S_0)}(X)=\frac{1}{8h_0^7w_8^{7i}}\cdot \frac{X^8-h_0^8}{X-h_0w_8^i}$ ，起算其分母在challenge $y$ 处的evaluation值： $\text{LiS0dens }\{8h_0^7w_8^{7i}\cdot (y-h_0w_8^i)\}_{i\in[8]}$
fflonk_prove.js中的computeLiS1，为对langrange多项式： $L_i^{(S_1)}(X)=\frac{1}{4h_1^3w_4^{3i}}\cdot \frac{X^4-h_1^4}{X-h_1w_4^i}$ ，起算其分母在challenge $y$ 处的evaluation值： $\text{LiS1dens }\{4h_1^3w_4^{3i}\cdot (y-h_1w_4^i)\}_{i\in[4]}$
fflonk_prove.js中的computeLiS2的当前算法复杂，将在
Polygon zkEVM Cryptography Security Review (Pil-stark and Snarkjs cryptography review), June 20, 2023 的5.1.4节展示优化版本。
计算 $\text{LiS0dens }$ 、 $\text{LiS1dens }$ 、 $\text{LiS2dens }$ 这些分母项的倒数。
fflonk_verify.js中的computeR0：
- 对多项式 $C_0(X)=q_L(X^8)+X\cdot q_R(X^8)+X^2\cdot q_O(X^8)+X^3\cdot q_M(X^8)+X^4\cdot q_C(X^8)+X^5\cdot S_{\sigma_1}(X^8)+X^6\cdot S_{\sigma_2}(X^8)+X^7\cdot S_{\sigma_3}(X^8)$ ，根据已知evaluation值： $q_L(h_0^8)=q_L(\xi),q_R(h_0^8)=q_R(\xi),q_O(h_0^8)=q_O(\xi),S_{\sigma_1}(h_0^8)=S_{\sigma_1}(\xi),S_{\sigma_2}(h_0^8)=S_{\sigma_2}(\xi),S_{\sigma_3}(h_0^8)=S_{\sigma_3}(\xi)$ ，按 $C_0(e)=q_L(\xi)+e\cdot q_R(\xi)+e^2\cdot q_O(\xi)+e^3\cdot q_M(\xi)+e^4\cdot q_C(\xi)+e^5\cdot S_{\sigma_1}(\xi)+e^6\cdot S_{\sigma_2}(\xi)+e^7\cdot S_{\sigma_3}(\xi)$ 方式，对于每个 $e\in S_0=h_0H_8$ ，来计算evaluation值集合： $\{C_0(e)| e\in S_0=h_0H_8\}$ 。
- 计算分子： $y^8-h_0^8$ 。
- 为计算 $r_0(X)=\sum_{i=0}^{7}C_0(h_0w_8^i)L_i^{(S_0)}(X)=\sum_{i=0}^{7}C_0(h_0w_8^i)\frac{1}{8h_0^7w_8^{7i}\cdot (X-h_0w_8^i)}\cdot (X^8-h_0^8)$ 在 $y$ 处的evaluation值，其中 $\{\frac{1}{8h_0^7w_8^{7i}\cdot (y-h_0w_8^i)}\cdot (y^8-h_0^8)\}_{i\in[8]}$ 可根据 $\text{LiS0dens}$ 的倒数值和分子 $y^8-h_0^8$ 值，来计算，然后将其与 $C_0(X)$ 在 $S_0$ 各点值相乘求和，即可计算出 $r_0(y)$ 值。
fflonk_verify.js中的computeR1，采用类似的策略，计算 $r_1(X)=\sum_{i=0}^{3}C_1(h_1w_4^i)L_i^{(S_1)}(X)$ 在 $y$ 处的evaluation值，其中 $C_1(X):=a(X^4)+X\cdot b(X^4)+X^2\cdot c(X^4)+X^3\cdot T_0(X^4)$ 。已知evaluation值 $a(h_1^4)=a(\xi),b(h_1^4)=b(\xi),c(h_1^4)=c(\xi),T_0(h_1^4)=T_0(\xi)$ ，即可计算出 $r_1(y)$ 值。
fflonk_verify.js中的computeR1，采用类似的策略，来计算 $r_1(y)$ ，将在
Polygon zkEVM Cryptography Security Review (Pil-stark and Snarkjs cryptography review), June 20, 2023 的5.1.4节展示优化版本。

2.1 优化建议：优化fflonk verifier中的 $r_0,r_1$ 计算

之前的 $r_0(y),r_1(y),r_2(y)$ 计算，需要 $O(n^2)$ 次域运算，其中 $n$ 为相应的domain size。本节，提供了一种策略，将该计算量降低为 $O(n\log n)$ ，并降低了常量因子的子。相比于https://github.com/iden3/snarkjs/blob/e70271b30a6392a34f664e6bececddad284c8854/test/fflonk/verifier.sol中的现有方案，本解决方案，需要更少的乘法、加法和模运算：

1）使用基于cosets的不同形式的lagrange多项式，可避免在computeChallenges中计算其coset values。此外，去掉了在旧的拉格朗日多项式表示中计算复杂分母所需的多余乘法和模运算。这是特别相关的，因为可将底层group元素作为常量访问。
2）当直接使用拉格朗日插值来evaluate $r (y)$ 时，需要拉格朗日多项式在 $y$ 处的完整evaluation值。与之不同，本方案直接使用barycentric evaluation，可降低重复计算。更值得注意的是，在之前的方案中，可将求和中重复的乘因子移除出去。
3）大量的计算量节约，源自，将（evaluate $r (y)$ 所需的） $C (X)$ 多项式evaluations，看成是某快速傅里叶转换。将该快速傅里叶转换展开为具体运算，可将乘法次数由 $n^2$ 降低为 $n\log n$ 。
4）对 $r_2$ 多项式的处理方式有所不同，有望进一步降低某些运算次数。

2.1.1 lagrange多项式及基于Cosets的Barycentric Evaluation

令group $H =< w >$ 的order为 $n$ ，某coset表示为 $c$ 。用于 $r_0,r_1$ 计算的策略是通用的。
将 $r (X)$ 多项式看成是degree小于 $n$ 的多项式，其基于coset $cH$ 的evaluation值为 $r_0,r_1,\cdots,r_{n-1}$ ，将 $r (X)$ 的以lagrange basis唯一表示为：
$r(X)=\sum_{i=0}^{n-1}r_i\cdot L_i^{(cH)}(X)$

注意，可将对coset $cH$ 的lagrange多项式，表示为对 $H$ 的具有shifted $X^{'} = X / c$ 的lagrange多项式：
$r(X)=\sum_{i=0}^{n-1}r_i\cdot L_i^{(cH)}(X)=\sum_{i=0}^{n-1}r_i\cdot L_i^{(H)}(\frac{X}{c})=\sum_{i=0}^{n-1}\frac{r_iw^i}{n}\cdot \frac{(X')^n-1}{X'-w^i}=\frac{(X')^n-1}{n}\cdot \sum_{i=0}^{n-1}\frac{r_i\cdot w^i}{X'-w^i}$

对于 $y\in\mathbb{F}$ ，为计算：
$r(y)=\frac{(y')^n-1}{n}\cdot \sum_{i=0}^{n-1}\frac{r_i\cdot w^i}{y'-w^i}$

需计算以下值：

1） $y^{'} = y / c$ ：对 $c$ 做一次倒数运算，然后与 $y$ 做一次乘法运算获得。
2） $\frac{(y')^n-1}{n}$ ：将 $1/ n$ 看成是已知常量值，需一次减法， $\log(n)+1$ 次乘法运算。
3） $\{r_i\cdot w^i\}_{i\in [n]}$ ：为特例情况，将在下一小节展开。
4） $\{\frac{1}{y'-w^i}\}_{i\in [n]}$ ：将 $w^i$ 看成是已知常量值，需要 $n$ 次减法，和 $n$ 次倒数运算，当然，所有的倒数运算可batch处理，可进一步降低计算量。

2.1.2 高效计算分子 $\{r_i\cdot w^i\}_{i\in [n]}$

可使用快速傅里叶变换来高效计算分子 $\{r_i\cdot w^i\}_{i\in [n]}$ 。
某多项式 $C (X)$ 形如：
$C(X):=f_0(X^n)+X\cdot f_1(X^n)+X^2\cdot f_2(X^n)+\cdots +X^{n-1}\cdot f_{n-1}(X^n)$

忽略 $r_2(X)$ 特例情况，FFlonk verifier，对于 $i\in [n]$ ，有：
$r_i\cdot w^i=C(c\cdot w^i)\cdot w^i$

其中：
$C(c\cdot w^i)\cdot w^i=f_0(c^n)\cdot w^i+(c\cdot w^i)\cdot f_1(c^n)\cdot w^i+\cdots +(c\cdot w^i)^{n-1}\cdot f_{n-1}(c^n)\cdot w^i=f_0(c^n)\cdot w^i+(c \cdot f_1(c^n))\cdot (w^i)^2+\cdots +(c^{n-1} \cdot f_{n-1}(c^n))$

其等价为对多项式：
$C'(X)=c^{n-1}f_{n-1}(c^n)+f_0(c^n)X+cf_1(c^n)X^2+\cdots+c^{n-2}f_{n-2}(c^n)X^{n-1}$
在 $w^i$ 处的evaluation值。

因此，为计算该分子，需：

1）计算系数向量： $\text{coeffs}[c^{n-1}f_{n-1}(c^n),f_0(c^n),cf_1(c^n),\cdots,c^{n-2}f_{n-2}(c^n)]$ 。假设计算 $r_0,r_1$ 时已知 $f_0(c^n),f_1(c^n),\cdots,f_{n-1}(c^n)$ ，则计算该系数向量需 $2 n$ 次乘法运算。
2）以 $n\log n$ 次乘法运算来计算基于 $H$ 的快速傅里叶变换：
$C'(H)=FFT(H,coeffs)=\{r_i\cdot w^i\}_{i\in [n]}$
从而可获得相应的分子。具体来说，在solidity中可将该快速傅里叶变换展开为具体的运算。
这将避免当前，在verifier中的重复计算coset中每个domain点的 $C (X)$ 值。

2.2 优化FFlonk verifier中的 $r_2$ 计算

2.2.1 对union of two cosets的lagrange多项式

令group $H =< w >$ 的order为 $n$ ，某coset表示为 $c_1,c_2$ 。定义集合 $S=c_1H\cup c_2H$ 为2个cosets集合，定义如下常量：
$s_1(c_1/c_2)^n-1,s_2(c_2/c_1)^n-1$
定义 $X_1'=X/c_1,X_2'=X/c_2$ ，对于 $i\in [2n]$ ，定义对 $S$ 的lagrange多项式为：

若 $i\in [0,n-1]$ ，则有： $L_i^{(S)}(X)=\frac{w^i}{n}\cdot \frac{(X_1')^n-1}{(X_1'-w^i)}\cdot \frac{(X_2')^n-1}{s_1}$
若 $i\in [n,2n-1]$ ，则有： $L_i^{(S)}(X)=\frac{w^i}{n}\cdot \frac{(X_2')^n-1}{(X_2'-w^i)}\cdot \frac{(X_1')^n-1}{s_2}$

2.2.2 对 $S$ 的barycentric evaluation

令 $r (X)$ 为degree小于 $n$ 的多项式，其基于 $S$ 的evaluation值为 $r_0,r_1,\cdots,r_{2n-1}$ ，可基于 $S$ 的lagrange basis，将 $r (X)$ 唯一表示为：
$r(X)=\sum_{i\in[2n]}r_i\cdot L_i^{(S)}(X)=\sum_{i\in[n]}r_i\cdot L_i^{(S)}(X) + \sum_{i\in[n,2n-1]}r_i\cdot L_i^{(S)}(X)\\ =\sum_{i=0}^{n-1}\frac{w^i}{n}\cdot \frac{(X_1')^n-1}{(X_1'-w^i)}\cdot \frac{(X_2')^n-1}{s_1} + \sum_{j=n}^{2n-1}\frac{w^j}{n}\cdot \frac{(X_2')^n-1}{(X_2'-w^j)}\cdot \frac{(X_1')^n-1}{s_2}\\=\frac{[(X_1')^n-1]\cdot [(X_2')^n-1]}{n}\cdot (\frac{1}{s_1}\cdot \sum_{i=0}^{n-1}\frac{r_i\cdot w^i}{(X_1'-w^i)}+\frac{1}{s_2}\cdot \sum_{j=n}^{2n-1}\frac{r_j\cdot w^j}{(X_2'-w^j)})$

对于 $y\in\mathbb{F}$ ，为计算：
$r(y)=\frac{[(y_1')^n-1]\cdot [(y_2')^n-1]}{n}\cdot (\frac{1}{s_1}\cdot \sum_{i=0}^{n-1}\frac{r_i\cdot w^i}{(y_1'-w^i)}+\frac{1}{s_2}\cdot \sum_{j=n}^{2n-1}\frac{r_j\cdot w^j}{(y_2'-w^j)})$
需计算以下值：

1） $y_1'=y/c_1,y_2'=y/c_2$ ：策略与2.1.1节的类似。
2） $\frac{[(y_1')^n-1]\cdot [(y_2')^n-1]}{n}$ ：策略与2.1.1节的类似。
3） $s_1=(c_1/c_2)^n-1,s_2=(c_2/c_1)^n-1$ ：计算 $c_1/c_2)^n$ ，对其求倒数获得 $c_2/c_1)^n$ ，然后各自减一。
4） $1/s_1,1/s_2$ ：对 $s_1,s_2$ 求倒数。
5） $\{r_i\cdot w^i\}_{i\in[n]},\{r_j\cdot w^j\}_{j\in[n,2n-1]}$ ：需采用2.1.2节的类似技巧来高效计算分子，其包含对 $H$ 的2次傅里叶变换。
6） $\{1/(y_1'-w^i)\}_{i\in[n]},\{1/(y_2'-w^j)\}_{j\in[n,2n-1]}$ ：策略与2.1.1节的类似。

2.3 基于smooth domain的FFT

令 $\mathbb{F}$ 为FFT-friendly域，其order $|\mathbb{F}|=p$ ，使得 $(p - 1)$ 可整除powers of small primes。如，FFlonk verifier solidity合约内的scalar域，其order为：
$p = 21888242871839275222246405745257275088548364400416034343698204186575808495617$
有：
$p-1)=2^{28}*3^2*13*29*983*11003*237073*405928799*1670836401704629* 13818364434197438864469338081$

令 $H\subseteq\mathbb{F}^{\times}$ 为某smooth乘法子群，其order $∣ H ∣ = n$ 可分解为powers of small primes。
Cooley-Tukey FFT算法或（对有限域通用的）Mixed-radix FFT算法，可基于 $H$ 对degree小于 $n$ 的多项式进行高效evaluate。其domain size无需为a power of two，可将其看成是某FFT算法。

其核心思想为：

不同于每个FFT step中的对半（奇偶）切分，而是将每个FFT step切分为 $r$ 个chunks，其中 $r$ 为当前迭代的素因子。通常将 $r$ 称为radix。如，标准FFT的radix总为2。

在线上SageMath脚本执行环境中执行如下代码：

p = 21888242871839275222246405745257275088548364400416034343698204186575808495617
# F is a field whose size is p
F = GF(p)
# define x to be a variable
R.<x> = PolynomialRing(F)

# g is the multiplicative generator with order (p-1)
g = F.multiplicative_generator()

# Implementation of Mixed-Radix FFT / Cooley-Tukey FFT for finite fields
# vals are the coeffs
# w is the multiplicative generator
# For n = p1^e1 * p2^e2 * ..., factors = [(p1, e1), (p2, e2), ...]
def fft_helper(vals, w, factors):
	n = len(vals)
	# Base case
	if n == 1:
		return vals
		
	# r is the current radix
	r = factors[0][0]
	
	# decrementing exponent
	# [(p1, e1), (p2, e2), ...]
	# goes to [(p1, e1-1), (p2, e2), ...]
	# or [(p2, e2), ...] if e1 = 1
	factors = list(factors)
	factors[0] = (factors[0][0], factors[0][1]-1)
	if factors[0][1] == 0:
		factors = factors[1:]
		
	# Recursive Step
	# acc = [fft0, fft1, fft2, ...]
	acc = list()
	for i in range(r):
		acc.append(fft_helper(vals[i::r], w^r, factors))
		
	# Allocate list for evaluations
	evals = [0 for _ in range(n)]
	
	# Accumulate recursive steps
	shift = n // r
	# acc = [fft0, fft1, fft2, ...]
	# Iterate as [t0 = (fft0[0], fft1[0], fft2[0]), t1 = (fft0[1], fft1[1], fft2[1]),...]
	for i, t in enumerate(zip(*acc)):
		for j in range(r):
		evals[i + shift * j] = sum((w^(i+shift*j))^k * t[k] for k in range(r))
	return evals
	
# w is the multiplicative generator
# vals are the coeffs
def fft(w, vals):
	factors = list(factor(w.multiplicative_order()))
	return fft_helper(vals, w, factors)
	
# Domain Size
n = 2^2 * 3^2 * 13
# w is generator for group of size n
w = g^((p-1)/n)

coeffs = [F.random_element() for _ in range(n)]
poly = sum([c*x^i for i,c in enumerate(coeffs)])
fft_evals = fft(w, coeffs)
true_evals = [poly(w^i) for i in range(n)]
print(fft_evals == true_evals)

以上代码中，使用的是FFlonk verifier中的scalar域，并取其subgroup size $n=2^2*3^2*13$ ，对应的radix为 $2, 3, 13$ 。
该代码可调整为任意subgroup size和radix，但若已知某特定subgroup因式分解，则可定制该代码为特定radix。当subgroup order $n$ 为a power of $2$ 时，该算法与标准FFT算法一致。

3. 定制PIL约束

在本轮审计中，验证了STARK verifier circom descripto中所使用的定制模板中相应定制PIL约束的正确性。
该proving pipeline流程为：

1）使用包含定制模板的circom descriptor来描述某STARK verifier。
2）将STARK verifier circom descriptor编译为R1CS约束，其中定制模板不生成约束。
3）将Circom-generated R1CS约束，编写为某STARK verifier PIL descriptor的PIL约束。
4）添加定制PIL约束，以完善STARK verifier PIL descriptor，以捕获circom定制模板。

本轮审计主要关注以上proving pipeline流程的第4）步，检查3个定制模板的定制PIL约束：

Fp3乘法（circuits.gl/cmul.circom）
1-of-4 Fp3 selector（circuits.gl/treeselector4.circom）
Poseidon hash evaluation（circuits.gl目录下的poseidon.circom和poseidon_constants.circom）

本轮审计还包含，（指向committed witness多项式数量的）与C12和C18 PIL descriptors相关的PIL约束文件：

src/compressor目录下的compressor12.pil.ejs和compressor12_setup.js
src/compressor目录下的compressor18.pil.ejs和compressor18_setup.js

所审计的pil-stark版本为：

https://github.com/0xPolygonHermez/pil-stark/tree/bb2404acd221083f0efa0b3ab8fa50bb421e686e

3.1 Fp3乘法（CMUL）

未发现问题。

3.2 1-of-4 Fp3 selector（TreeSelector4）

注意，后面的“keys”约束，是指binary约束。

https://github.com/0xPolygonHermez/pil-stark/blob/bb2404acd221083f0efa0b3ab8fa50bb421e686e/src/compressor/compressor18.pil.ejs#L326，326、335、344、353中的注释，声明了：

仅当key bit多项式a[12]和a[13]设置了合适的bits，相应的“keys”多项式才为1。

若限制了a[12]和a[13]为 ${0,1\}$ ，则以上声明成立。但该约束并不存在。事实上，TreeSelector4函数中也未限制a[12]和a[13]。

更准确来说，现有PIL约束提供了如下bijection：

keys1=0 a[12]=1或a[13]=1
keys2=0 a[12]=0或a[13]=1
keys3=0 a[12]=1或a[13]=0
keys4=0 a[12]=0或a[13]=0

这强化了：不存在a[12]和a[13]的设置，使得所有 keys{1,2,3,4} 都可设置为0，这意味着该selector应用于至少一个分支。

建议：

为不失通用性，将注释替换为：“keys1 will be non-zero unless a[12]=0 or a[13]=0”。

3.3 Poseidon hash evaluation（Poseidon/CustPoseidon）

3.3.1 可能的可靠性问题：对CustPoseidon Selection Key无binary enforcement

在https://github.com/0xPolygonHermez/pil-stark/blob/bb2404acd221083f0efa0b3ab8fa50bb421e686e/src/compressor/compressor18.pil.ejs#L93的93-100行，witness多项式a[8]用于选择custom Poseidon input order，但其假设了a[8]值为二进制值 ${0,1\}$ 。若未约束a[8]，prover可设置a[8]为其它值。不清楚在该上下文中，利用该漏洞，是否会引起Poseidon的碰撞问题。

建议：

当读取custom Poseidon输入时，添加约束a[8]为二进制值的多项式约束：
$POSE I D ONC U STF I RST * (1 - a [8]) * a [8] = 0$

3.3.2 优化建议：POSEIDONM多项式可用于capture multiple rounds

在https://github.com/0xPolygonHermez/pil-stark/blob/bb2404acd221083f0efa0b3ab8fa50bb421e686e/src/compressor/compressor18.pil.ejs#L196第196行，POSEIDONM、POSEIDONAFTERPART、POSEIDONFIRST、POSEIDONCUSTFIRST用作selector，用于将poseidonM多项式中的值传递给下一行。而POSEIDONAFTERPART、POSEIDONFIRST和POSEIDONCUSTFIRST具有其它重要功能，POSEIDONM仅用于当Round 30 -> output row的情况。

建议：