Halo2 學習筆記——背景資料之Polynomials（3）

1. 引言

令 A ( X ) A(X) A(X)為基于 F p \mathbb{F}_p Fp​的degree-3多項式：

A ( X ) = a 0 + a 1 X + a 2 X 2 + a 3 X 3 A(X)=a_0+a_1X+a_2X^2+a_3X^3 A(X)=a0​+a1​X+a2​X2+a3​X3

其中 a 0 a_0 a0​稱為constant term。

degree n − 1 n-1 n−1的多項式，其有 n n n個系數。

通常，将變量 X X X替換為具體的值 x x x，來計算相應的結果 A ( x ) A(x) A(x)——數學上稱為“evaluating A ( X ) A(X) A(X) at a point x x x”。（注意，此處的point不要與橢圓曲線上的point混淆。）

多項式的主要屬性有：

• d e g ( A ( X ) B ( X ) ) = d e g ( A ( X ) ) + d e g ( B ( X ) ) , d e g ( A ( X ) / B ( X ) ) = d e g ( A ( X ) ) − d e g ( B ( X ) ) deg(A(X)B(X))=deg(A(X))+deg(B(X)), deg(A(X)/B(X))=deg(A(X))-deg(B(X)) deg(A(X)B(X))=deg(A(X))+deg(B(X)),deg(A(X)/B(X))=deg(A(X))−deg(B(X))
• 若 A ( X ) A(X) A(X)的degree為 n − 1 n-1 n−1，則對該多項式的 n n n個不同點進行evaluation，這些evaluation值可完全定義該多項式。換句話說，由這 n n n個不同點的evaluation值，通過多項式插值獲得唯一的degree為 n − 1 n-1 n−1的多項式 A ( X ) A(X) A(X)。

• 多項式 A ( X ) A(X) A(X)的系數表示法為： [ a 0 , a 1 , ⋯   , a n − 1 ] [a_0,a_1,\cdots,a_{n-1}] [a0​,a1​,⋯,an−1​]，其點值表示法為：

  [ ( x 0 , A ( x 0 ) ) , ( x 1 , A ( x 1 ) ) , ⋯   , ( x n − 1 , A ( x n − 1 ) ) ] [(x_0,A(x_0)), (x_1,A(x_1)), \cdots, (x_{n-1},A(x_{n-1}))] [(x0​,A(x0​)),(x1​,A(x1​)),⋯,(xn−1​,A(xn−1​))]

  其中 x 0 , x 1 , ⋯   , x n − 1 x_0,x_1,\cdots,x_{n-1} x0​,x1​,⋯,xn−1​為 n n n個不同的值。

  這兩種方式都可唯一定義同一多項式。

2. Horner’s rule

可借助Horner’s rule來對一個 n − 1 n-1 n−1 degree多項式進行高效evaluation，僅需 n − 1 n-1 n−1個乘法運算 + n − 1 n-1 n−1個加法運算：

a 0 + a 1 X + a 2 X 2 + ⋯ + a n − 1 X n − 1 = a 0 + X ( a 1 + X ( a 2 + ⋯ + X ( a n − 2 + X a n − 1 ) ) ) a_0+a_1X+a_2X^2+\cdots+a_{n-1}X^{n-1}=a_0+X(a_1+X(a_2+\cdots+X(a_{n-2}+Xa_{n-1}))) a0​+a1​X+a2​X2+⋯+an−1​Xn−1=a0​+X(a1​+X(a2​+⋯+X(an−2​+Xan−1​)))

3. FFT

FFT可将多項式的系數表示法 和 點值表示法 之間高效互相轉換。

點值表示法時，evaluate的point為 n n n-th roots of unity { w 0 , w 1 , ⋯   , w n − 1 } \{w^0,w^1,\cdots,w^{n-1}\} {w0,w1,⋯,wn−1}，其中 w w w為a primitive n n n-th root of unity。

根據roots of unity的對城鄉，FFT每輪運算将reduce the evaluation into a problem only half the size。是以，通常使用的多項式長度為some power of two，即 n = 2 k n=2^k n=2k，進而可apply the halving reduction recursively。

3.1 将FFT用于快速多項式乘法

如需計算 A ( X ) ⋅ B ( X ) = C ( X ) A(X)\cdot B(X)=C(X) A(X)⋅B(X)=C(X)，若采用系數表示法時，需要進行 O ( n 2 ) O(n^2) O(n2)次運算：

A ( X ) = a 0 + a 1 X + a 2 X 2 + ⋯ + a n − 1 X n − 1 A(X)=a_0+a_1X+a_2X^2+\cdots+a_{n-1}X^{n-1} A(X)=a0​+a1​X+a2​X2+⋯+an−1​Xn−1

B ( X ) = b 0 + b 1 X + b 2 X 2 + ⋯ + b n − 1 X n − 1 B(X)=b_0+b_1X+b_2X^2+\cdots+b_{n-1}X^{n-1} B(X)=b0​+b1​X+b2​X2+⋯+bn−1​Xn−1

C ( X ) = A ( X ) ⋅ B ( X ) = a 0 ⋅ ( b 0 + b 1 X + b 2 X 2 + ⋯ + b n − 1 X n − 1 ) + a 1 X ⋅ ( b 0 + b 1 X + b 2 X 2 + ⋯ + b n − 1 X n − 1 ) + ⋯ + a n − 1 X n − 1 ⋅ ( b 0 + b 1 X + b 2 X 2 + ⋯ + b n − 1 X n − 1 ) C(X)=A(X)\cdot B(X)=a_0\cdot (b_0+b_1X+b_2X^2+\cdots+b_{n-1}X^{n-1})\\+ a_1X\cdot (b_0+b_1X+b_2X^2+\cdots+b_{n-1}X^{n-1})\\ +\cdots \\+ a_{n-1}X^{n-1}\cdot (b_0+b_1X+b_2X^2+\cdots+b_{n-1}X^{n-1}) C(X)=A(X)⋅B(X)=a0​⋅(b0​+b1​X+b2​X2+⋯+bn−1​Xn−1)+a1​X⋅(b0​+b1​X+b2​X2+⋯+bn−1​Xn−1)+⋯+an−1​Xn−1⋅(b0​+b1​X+b2​X2+⋯+bn−1​Xn−1)

若采用點值表示法，則多項式乘法僅需 O ( n ) O(n) O(n)次運算：

A ： { ( x 0 , A ( x 0 ) ) , ( x 1 , A ( x 1 ) ) , ⋯   , ( x n − 1 , A ( x n − 1 ) ) } A：\{(x_0,A(x_0)), (x_1,A(x_1)),\cdots,(x_{n-1}, A(x_{n-1}))\} A：{(x0​,A(x0​)),(x1​,A(x1​)),⋯,(xn−1​,A(xn−1​))}

B ： { ( x 0 , B ( x 0 ) ) , ( x 1 , B ( x 1 ) ) , ⋯   , ( x n − 1 , B ( x n − 1 ) ) } B：\{(x_0,B(x_0)), (x_1,B(x_1)),\cdots,(x_{n-1}, B(x_{n-1}))\} B：{(x0​,B(x0​)),(x1​,B(x1​)),⋯,(xn−1​,B(xn−1​))}

C ： { ( x 0 , A ( x 0 ) B ( x 0 ) ) , ( x 1 , A ( x 1 ) B ( x 1 ) ) , ⋯   , ( x n − 1 , A ( x n − 1 ) B ( x n − 1 ) ) } C：\{(x_0,A(x_0)B(x_0)), (x_1,A(x_1)B(x_1)),\cdots,(x_{n-1}, A(x_{n-1})B(x_{n-1}))\} C：{(x0​,A(x0​)B(x0​)),(x1​,A(x1​)B(x1​)),⋯,(xn−1​,A(xn−1​)B(xn−1​))}

其中 C C C中的值是直接multiplied pointwise的。

進行快速多項式乘法運算的主要步驟為：

• 1）Evaluate polynomials at all n n n points。
• 2）Perform fast pointwise multiplication in the evaluation representation O ( n ) O(n) O(n)。
• 3）将結果再轉換為系數表示法（通過IFFT）。

主要難點在于高效的對多項式進行evaluate和插值。

直覺地，evaluate a polynomial at n n n points需要 O ( n 2 ) O(n^2) O(n2)次運算（每個point采用Horner’s rule需 O ( n ) O(n) O(n)次運算， n n n個point對應為 O ( n 2 ) O(n^2) O(n2)次運算）：

[ A ( 1 ) A ( w ) A ( w 2 ) ⋮ A ( w n − 1 ) ] = [ 1 1 1 ⋯ 1 1 w w 2 ⋯ w n − 1 1 w 2 w 2 ⋅ 2 ⋯ w 2 ⋅ ( n − 1 ) ⋮ ⋮ ⋮ ⋮ 1 w n − 1 w 2 ( n − 1 ) ⋯ w ( n − 1 ) 2 ] ⋅ [ a 0 a 1 a 2 ⋮ a n − 1 ] \begin{bmatrix} A(1) \\ A(w) \\ A(w^2) \\ \vdots \\ A(w^{n-1}) \end{bmatrix} = \begin{bmatrix} 1& 1 & 1 & \cdots & 1\\ 1& w & w^2 & \cdots & w^{n-1} \\ 1& w^2 & w^{2\cdot 2} & \cdots & w^{2\cdot (n-1)} \\ \vdots & \vdots & \vdots & &\vdots \\ 1 & w^{n-1} & w^{2(n-1)} & \cdots & w^{(n-1)^2} \end{bmatrix} \cdot \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎡​A(1)A(w)A(w2)⋮A(wn−1)​⎦⎥⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎢⎡​111⋮1​1ww2⋮wn−1​1w2w2⋅2⋮w2(n−1)​⋯⋯⋯⋯​1wn−1w2⋅(n−1)⋮w(n−1)2​⎦⎥⎥⎥⎥⎥⎤​⋅⎣⎢⎢⎢⎢⎢⎡​a0​a1​a2​⋮an−1​​⎦⎥⎥⎥⎥⎥⎤​

可将以上多項式表示為：

A ^ = V w ⋅ A \hat{\mathbf{A}}=\mathbf{V}_w\cdot \mathbf{A} A^=Vw​⋅A

其中 A ^ \hat{\mathbf{A}} A^可稱為 A \mathbf{A} A的離散傅裡葉變換（DFT）， V w \mathbf{V}_w Vw​稱為Vandermonde matrix。

3.2 The (radix-2) Cooley-Tukey algorithm

将size為 n n n的DFT分為2個交錯的size為 n / 2 n/2 n/2的DFT。

對于多項式 A ( X ) = a 0 + a 1 X + a 2 X 2 + ⋯ + a n − 1 X n − 1 A(X)=a_0+a_1X+a_2X^2+\cdots+a_{n-1}X^{n-1} A(X)=a0​+a1​X+a2​X2+⋯+an−1​Xn−1，可将其分為奇數項和偶數項：

A e v e n = a 0 + a 2 X + ⋯ + a n − 2 X n / 2 − 1 A_{even}=a_0+a_2X+\cdots+a_{n-2}X^{n/2-1} Aeven​=a0​+a2​X+⋯+an−2​Xn/2−1

A o d d = a 1 + a 3 X + ⋯ + a n − 1 X n / 2 − 1 A_{odd}=a_1+a_3X+\cdots +a_{n-1}X^{n/2-1} Aodd​=a1​+a3​X+⋯+an−1​Xn/2−1

然後基于此，有 A ( X ) = A e v e n ( X 2 ) + X A o d d ( X 2 ) A(X)=A_{even}(X^2)+XA_{odd}(X^2) A(X)=Aeven​(X2)+XAodd​(X2)

對 A ( X ) A(X) A(X) evaluate at points w n i , w n n 2 + i w_n^i, w_n^{\frac{n}{2}+i} wni​,wn2n​+i​，其中 i ∈ [ 0 , ⋯   , n 2 − 1 ] i\in[0,\cdots,\frac{n}{2}-1] i∈[0,⋯,2n​−1]，有：

A ( w n i ) = A e v e n ( ( w n i ) 2 ) + w n i A o d d ( ( w n i ) 2 ) A(w_n^i)=A_{even}((w_n^i)^2)+w_n^iA_{odd}((w_n^i)^2) A(wni​)=Aeven​((wni​)2)+wni​Aodd​((wni​)2)

A ( w n n 2 + i ) = A e v e n ( ( w n n 2 + i ) 2 ) + w n n 2 + i A o d d ( ( w n n 2 + i ) 2 ) = A e v e n ( ( − w n i ) 2 ) − w n i A o d d ( ( − w n i ) 2 ) 【 ← ( 根 據 n e g a t i o n l e m m a ) 】 = A e v e n ( ( w n i ) 2 ) − w n i A o d d ( ( w n i ) 2 ) A(w_n^{\frac{n}{2}+i})=A_{even}((w_n^{\frac{n}{2}+i})^2)+w_n^{\frac{n}{2}+i}A_{odd}((w_n^{\frac{n}{2}+i})^2)\\ =A_{even}((-w_n^i)^2)-w_n^iA_{odd}((-w_n^i)^2) 【\leftarrow (根據negation lemma)】\\ =A_{even}((w_n^i)^2)-w_n^iA_{odd}((w_n^i)^2) A(wn2n​+i​)=Aeven​((wn2n​+i​)2)+wn2n​+i​Aodd​((wn2n​+i​)2)=Aeven​((−wni​)2)−wni​Aodd​((−wni​)2)【←(根據negationlemma)】=Aeven​((wni​)2)−wni​Aodd​((wni​)2)

由上可看出，兩者存在一定的對稱性。僅需evaluate A e v e n A_{even} Aeven​和 A o d d A_{odd} Aodd​ over half the domain ( w n 0 ) 2 , ( w n ) 2 , ⋯   , ( w n n 2 − 1 ) 2 (w_n^0)^2,(w_n)^2,\cdots,(w_n^{\frac{n}{2}-1})^2 (wn0​)2,(wn​)2,⋯,(wn2n​−1​)2，其中 i = [ 0 , ⋯   , n 2 − 1 ] i=[0,\cdots, \frac{n}{2}-1] i=[0,⋯,2n​−1]（根據halving lemma）。

即意味着，可将length- n n n DFT 轉換為 2個length- n 2 \frac{n}{2} 2n​ DFTs。

當 n = 2 k n=2^k n=2k時（為power of two，若不夠可zero-padding），然後遞歸使用divide-and-conquer政策。根據Master Theorem，相應的evaluation算法僅需 O ( n log ⁡ 2 n ) O(n\log_2n) O(nlog2​n)次運算，也稱為Fast Fourier Transform (FFT)。

3.3 Inverse FFT

至此，已完成polynomial evaluation 和 multiply pointwise，現在需要将 C ( X ) = A ( X ) ⋅ B ( X ) C(X)=A(X)\cdot B(X) C(X)=A(X)⋅B(X)由 點值表示法 轉換為 系數表示法。僅需，在點值表示法 上進行FFT運算：

• 1）将Vandermonde matrix中的 w i w^i wi替換為 w − i w^{-i} w−i
• 2）将最終結果乘以 1 / n 1/n 1/n

即：

A = 1 n V w − 1 ⋅ A ^ \mathbf{A}=\frac{1}{n}\mathbf{V}_{w^{-1}}\cdot \hat{\mathbf{A}} A=n1​Vw−1​⋅A^

IFFT與FFT具有類似的格式，詳細參見 The Fast Fourier Transform and Polynomial Multiplication。

4. The Schwartz-Zippel lemma

The Schwartz-Zippel lemma表達的是“different polynomials are different at most points”：

令 p ( x 1 , x 2 , ⋯   , x n ) p(x_1,x_2,\cdots,x_n) p(x1​,x2​,⋯,xn​)為nonzero polynomial of n n n variables with degree d d d，令 S S S為a finite set of numbers with at least d d d elements in it。若從 S S S中随機選出 α 1 , α 2 , ⋯   , α n \alpha_1,\alpha_2,\cdots,\alpha_n α1​,α2​,⋯,αn​，則有：

Pr [ p ( α 1 , α 2 , ⋯   , α n ) = 0 ] ≤ d ∣ S ∣ \text{Pr}[p(\alpha_1,\alpha_2,\cdots,\alpha_n)=0]\leq \frac{d}{|S|} Pr[p(α1​,α2​,⋯,αn​)=0]≤∣S∣d​

對于degree為 d d d的單變量非零多項式 p ( X ) p(X) p(X)，即意味着最多有 d d d個根。

The Schwartz-Zippel lemma用于polynomial equality testing：

已知2個多變量多項式， p 1 ( x 1 , ⋯   , x n ) p_1(x_1,\cdots, x_n) p1​(x1​,⋯,xn​)的degree為 d 1 d_1 d1​， p 2 ( x 1 , ⋯   , x n ) p_2(x_1,\cdots,x_n) p2​(x1​,⋯,xn​)的degree為 d 2 d_2 d2​，可随機選擇 α 1 , ⋯   , α n ← S \alpha_1,\cdots,\alpha_n\leftarrow S α1​,⋯,αn​←S，其中 ∣ S ∣ ≥ ( d 1 + d 2 ) |S|\geq (d_1+d_2) ∣S∣≥(d1​+d2​)，測試 p 1 ( α 1 , ⋯   , α n ) − p 2 ( α 1 , ⋯   , α n ) = 0 p_1(\alpha_1,\cdots,\alpha_n)-p_2(\alpha_1,\cdots,\alpha_n)=0 p1​(α1​,⋯,αn​)−p2​(α1​,⋯,αn​)=0是否成立。若 p 1 , p 2 p_1,p_2 p1​,p2​這兩個多項式完全相同，必然成立。若兩者不同，則該等式成立的機率不高于 m a x ( d 1 , d 2 ) ∣ S ∣ \frac{max(d_1,d_2)}{|S|} ∣S∣max(d1​,d2​)​。

5. Vanishing polynomial

對于order n n n multiplicative subgroup H \mathcal{H} H with primitive root of unity w w w，對于所有的 w i ∈ H , i ∈ [ n − 1 ] w^i\in\mathcal{H},i\in[n-1] wi∈H,i∈[n−1]，有 ( w i ) n = ( w n ) i = ( w 0 ) i = 1 (w^i)^n=(w^n)^i=(w^0)^i=1 (wi)n=(wn)i=(w0)i=1，換句話說：

Z H ( X ) = X n − 1 = ( X − w 0 ) ( X − w 1 ) ( X − w 2 ⋯ ( X − w n − 1 ) Z_H(X)=X^n-1=(X-w^0)(X-w^1)(X-w^2\cdots(X-w^{n-1}) ZH​(X)=Xn−1=(X−w0)(X−w1)(X−w2⋯(X−wn−1)

H \mathcal{H} H中的每個元素為 Z H ( X ) Z_H(X) ZH​(X)的根。稱 Z H ( X ) Z_H(X) ZH​(X)為the vanishing polynomial over H \mathcal{H} H，因為其evaluate to zero on all elements of H \mathcal{H} H。

vanishing polynomial 用于驗證polynomial constraints時将特别實用。如，為了驗證 A ( X ) + B ( X ) = C ( X ) A(X)+B(X)=C(X) A(X)+B(X)=C(X) over H \mathcal{H} H，可改為驗證 A ( X ) + B ( X ) − C ( X ) A(X)+B(X)-C(X) A(X)+B(X)−C(X)為 some multiple of Z H ( X ) Z_H(X) ZH​(X)。換句話說，若将constraints除以vanishing polynomial，仍然可産生某多項式 A ( X ) + B ( X ) − C ( X ) Z H ( X ) = H ( X ) \frac{A(X)+B(X)-C(X)}{Z_H(X)}=H(X) ZH​(X)A(X)+B(X)−C(X)​=H(X)，則可說明 A ( X ) + B ( X ) − C ( X ) = 0 A(X)+B(X)-C(X)=0 A(X)+B(X)−C(X)=0 over H \mathcal{H} H。

6. Lagrange basis functions

多項式通常以monomial basis（如 X , X 2 , ⋯   , X n X,X^2,\cdots,X^n X,X2,⋯,Xn）來表示，但是，當機遇multiplicative subgroup of order n n n時，采用Lagrange basis表示更自然。

對于order- n n n multiplicative subgroup H \mathcal{H} H with primitive root of unity w w w，該subgroup的Lagrange basis為a set of functions { L i } i = 0 n − 1 \{\mathcal{L}_i\}_{i=0}^{n-1} {Li​}i=0n−1​，其中：

L i ( w j ) = { 1 if  i = j , 0 otherwise. \mathcal{L}_i(w^j)=\left\{\begin{matrix} 1 & \text{if } i=j,\\ 0 & \text{otherwise.} \end{matrix}\right. Li​(wj)={10​if i=j,otherwise.​

可将 L i ( w j ) \mathcal{L}_i(w^j) Li​(wj)更精簡的表示為 δ i j \delta_{ij} δij​，其中 δ \delta δ為the Kronecker delta function。

至此，可将多項式表示為a linear combination of Lagrange basis functions：

A ( X ) = ∑ i = 0 n − 1 a i L i ( X ) , X ∈ H A(X)=\sum_{i=0}^{n-1}a_i\mathcal{L}_i(X), X\in\mathcal{H} A(X)=∑i=0n−1​ai​Li​(X),X∈H

即等價為說， p ( X ) p(X) p(X) evaluates to a 0 a_0 a0​ at w 0 w^0 w0, to a 1 a_1 a1​ at w 1 w^1 w1, to a 2 a_2 a2​ at w 2 w^2 w2等等。

當基于multiplicative subgroup時，Lagrange basis function具有很友善的sparse表示形式：

L i ( X ) = c i ⋅ ( X n − 1 ) X − w i \mathcal{L}_i(X)=\frac{c_i\cdot (X^n-1)}{X-w^i} Li​(X)=X−wici​⋅(Xn−1)​

其中 c i c_i ci​為barycentric weight。詳細可參看Barycentric Lagrange Interpolation*。

當 i = 0 i=0 i=0時，有 c = 1 / n ⇒ L 0 ( X ) = 1 n ( X n − 1 ) X − 1 c=1/n\Rightarrow \mathcal{L}_0(X)=\frac{1}{n}\frac{(X^n-1)}{X-1} c=1/n⇒L0​(X)=n1​X−1(Xn−1)​。

對于evaluation point set { x 0 , x 1 , ⋯   , x n − 1 } \{x_0,x_1,\cdots,x_{n-1}\} {x0​,x1​,⋯,xn−1​}，若不假設 x i x_i xi​為multiplicative subgroup，仍可采用Lagrange 多項式 L i \mathcal{L}_i Li​來表示：

L i ( X ) = ∏ j ≠ i X − x j x i − x j , i ∈ [ 0.. n − 1 ] \mathcal{L}_i(X)=\prod_{j\neq i}\frac{X-x_j}{x_i-x_j}, i\in[0..n-1] Li​(X)=∏j​=i​xi​−xj​X−xj​​,i∈[0..n−1]

對于每個 X = x j ≠ x i X=x_j\neq x_i X=xj​​=xi​，都有零分子項 ( x j − x j ) (x_j-x_j) (xj​−xj​)，進而使整個值為0。若 X = x i X=x_i X=xi​，有 x i − x j x i − x j \frac{x_i-x_j}{x_i-x_j} xi​−xj​xi​−xj​​，結果為1。進而可實作desired Kronecker delta behaviour L i ( x j ) = δ i j \mathcal{L}_i(x_j)=\delta_{ij} Li​(xj​)=δij​ on the set { x 0 , x 1 , ⋯   , x n − 1 } \{x_0,x_1,\cdots,x_{n-1}\} {x0​,x1​,⋯,xn−1​}。

6.1 Lagrange interpolation

已知多項式的點值表示：

A : { ( x 0 , A ( x 0 ) ) , ( x 1 , A ( x 1 ) ) , ⋯   , ( x n − 1 , A ( x n − 1 ) ) } A: \{(x_0,A(x_0)), (x_1,A(x_1)), \cdots, (x_{n-1},A(x_{n-1}))\} A:{(x0​,A(x0​)),(x1​,A(x1​)),⋯,(xn−1​,A(xn−1​))}

可基于Lagrange basis來建構其 系數表示：

A ( X ) = ∑ i = 0 n − 1 A ( x i ) L i ( X ) A(X)=\sum_{i=0}^{n-1}A(x_i)\mathcal{L}_i(X) A(X)=∑i=0n−1​A(xi​)Li​(X)

其中 X ∈ { x 0 , x 1 , ⋯   , x n − 1 } X\in\{x_0,x_1,\cdots,x_{n-1}\} X∈{x0​,x1​,⋯,xn−1​}。

參考資料

[1] Halo2 背景資料之Polynomials