STARK/SNARK加速小技巧

1. 引言

STARK/SNARK中包含了大量的有限域运算，如：

• STARK中包含了对 roots of unity domain D = { 1 , ω , ⋯   , ω 2 k − 1 } D=\{1,\omega,\cdots,\omega^{2^k-1}\} D={1,ω,⋯,ω2k−1} 的批量 polynomial evaluation运算

除可通过CPU/GPU等硬件加速外，还有一些算法可实现加速。

• CPU加速可参看：【监控CPU使用率用top】

  • CPU指令集——AVX2

• 采用CPU加速的密码学库有：

  • https://github.com/dalek-cryptography/ed25519-dalek

• 采用GPU加速的密码学库有：【监控GPU使用率：NVIDIA GPU可用nvidia-smi；Intel GPU可用intel-gpu-tools；AMD GPU可用fglrx或RadeonTop。】

  • https://github.com/filecoin-project/bellperson 以 520 core的GPU为例，其做multiexponentiation计算，与CPU的性能对比为：

$ RUST_LOG=info cargo test --features opencl -- --exact multiexp::gpu_multiexp_consistency --nocapture
   Compiling bellperson v0.18.0 (/home/administrator/bellperson)
    Finished test [unoptimized + debuginfo] target(s) in 8.21s
     Running unittests (target/debug/deps/bellperson-8b909898d5e3834a)

running 1 test
[2021-12-16T08:53:13Z INFO  bellperson::gpu::utils] Device: Device { vendor: Nvidia, name: "Quadro P620", memory: 2097479680, pci_id: PciId(768), uuid: Some(148dcfa4-76fd-6819-dfb2-6199f912527f), opencl: Some(Device { vendor: Nvidia, name: "Quadro P620", memory: 2097479680, pci_id: PciId(768), uuid: Some(148dcfa4-76fd-6819-dfb2-6199f912527f), device: Device { id: 139725486794256 } }) }
Testing Multiexp for 1024 elements...
[2021-12-16T08:53:13Z INFO  bellperson::gpu::locks] GPU is available for Multiexp!
[2021-12-16T08:53:13Z INFO  bellperson::gpu::program] Using kernel on OpenCL.
[2021-12-16T08:53:13Z INFO  bellperson::gpu::multiexp] Multiexp: 1 working device(s) selected. (CPU utilization: 0)
[2021-12-16T08:53:13Z INFO  bellperson::gpu::multiexp] Multiexp: Device 0: Quadro P620 (Chunk-size: 894826)
[2021-12-16T08:53:13Z INFO  bellperson::multiexp] GPU Multiexp kernel instantiated!
GPU took 187ms.
CPU took 208ms.
Speedup: x1.1122994
============================
Testing Multiexp for 2048 elements...
GPU took 52ms.
CPU took 395ms.
Speedup: x7.5961537
============================
Testing Multiexp for 4096 elements...
GPU took 92ms.
CPU took 378ms.
Speedup: x4.1086955
============================
Testing Multiexp for 8192 elements...
GPU took 132ms.
CPU took 688ms.
Speedup: x5.212121
============================
Testing Multiexp for 16384 elements...
GPU took 237ms.
CPU took 1471ms.
Speedup: x6.206751
============================
Testing Multiexp for 32768 elements...
GPU took 281ms.
CPU took 2611ms.
Speedup: x9.291815
============================
Testing Multiexp for 65536 elements...
GPU took 454ms.
CPU took 6154ms.
Speedup: x13.555066
============================
test multiexp::gpu_multiexp_consistency ... ok

前序博客有：

• 十分简明易懂的FFT（快速傅里叶变换）
• Halo中的快速傅里叶（逆）变换算法(I)FFT

2. 采用FFT来加速roots of unity domain内的批量polynomial evaluation

有限域内的Fast Fourier Transform（FFT），又可称为Number Theory Transform (NTT)。FFT用于将多项式系数表示法，转换为点值表示法。

令多项式 f ( X ) = ∑ i = 0 d c i X i f(X)=\sum_{i=0}^{d}c_iX^i f(X)=∑i=0d​ci​Xi的degree 不高于 2 k − 1 2^k-1 2k−1，其系数 c i ∈ F p c_i\in\mathbb{F}_p ci​∈Fp​。 令 ω \omega ω为 2 k 2^k 2k-th root of unity，求以下所有polynomial evaluation 值： ( f ( ω i ) ) i = 0 2 k − 1 = ( f ( 1 ) , f ( ω ) , f ( ω 2 ) , … , f ( ω 2 k − 1 ) ) (f(\omega^i))_{i=0}^{2^k-1} = (f(1), f(\omega), f(\omega^2), \ldots, f(\omega^{2^k-1})) (f(ωi))i=02k−1​=(f(1),f(ω),f(ω2),…,f(ω2k−1))

解决方案有：

• 1）最直观的方法是，依次计算每个evaluation值。【令 N = 2 k N=2^k N=2k，算法复杂度为 O ( N 2 ) O(N^2) O(N2)。】【interpolate、evaluate、divide等naive运算可参看：https://github.com/aszepieniec/stark-anatomy/blob/master/code/univariate.py】
• 2）更明智的方法为，根据FFT的divide-and-conquer策略，将多项式分为奇数项和偶数项表示：【递归调用，令 N = 2 k N=2^k N=2k，算法复杂度为 O ( N ⋅ log ⁡ N ) O(N\cdot \log N) O(N⋅logN)。】【interpolate、evaluate、divide等借助FFT/IFFT加速运算可参看：https://github.com/aszepieniec/stark-anatomy/blob/master/code/ntt.py】 f ( X ) = f E ( X 2 ) + X ⋅ f O ( X 2 ) f(X)=f_E(X^2)+X\cdot f_O(X^2) f(X)=fE​(X2)+X⋅fO​(X2) 其中： f E ( X 2 ) = f ( X ) + f ( − X ) 2 = ∑ i = 0 d + 1 2 − 1 c 2 i X 2 i f_E(X^2)=\frac{f(X)+f(-X)}{2}=\sum_{i=0}^{\frac{d+1}{2}-1}c_{2i}X^{2i} fE​(X2)=2f(X)+f(−X)​=∑i=02d+1​−1​c2i​X2i f O ( X 2 ) = f ( X ) − f ( − X ) 2 X = ∑ i = 0 d + 1 2 − 1 c 2 i + 1 X 2 i f_O(X^2)=\frac{f(X)-f(-X)}{2X}=\sum_{i=0}^{\frac{d+1}{2}-1}c_{2i+1}X^{2i} fO​(X2)=2Xf(X)−f(−X)​=∑i=02d+1​−1​c2i+1​X2i 从而有： f ( ω i ) = f E ( ω 2 i ) + ω i ⋅ f O ( ω 2 i ) f(\omega^i)=f_E(\omega^{2i})+\omega^i\cdot f_O(\omega^{2i}) f(ωi)=fE​(ω2i)+ωi⋅fO​(ω2i)

递归调用，多项式系数 c i c_i ci​为参数values，则调用ntt返回的即为 ( f ( ω i ) ) i = 0 2 k − 1 = ( f ( 1 ) , f ( ω ) , f ( ω 2 ) , … , f ( ω 2 k − 1 ) ) (f(\omega^i))_{i=0}^{2^k-1} = (f(1), f(\omega), f(\omega^2), \ldots, f(\omega^{2^k-1})) (f(ωi))i=02k−1​=(f(1),f(ω),f(ω2),…,f(ω2k−1)) evaluation值：【代码见：https://github.com/aszepieniec/stark-anatomy/blob/master/code/ntt.py】

def ntt( primitive_root, values ):
    assert(len(values) & (len(values) - 1) == 0), "cannot compute ntt of non-power-of-two sequence"
    if len(values)