ZPrize：着眼于奖励

简介

这篇文章总结了使用 CUDA 优化多标量乘法的不同方法，这些方法是在 ZPrize 中提出的。这在某些证明系统 (zk-SNARKs) 中是一项重要的计算，其中需要对椭圆曲线上的大量点进行相加。这些大型总和可以分解为较小的总和，每个总和都可以由处理器并行计算，从而使 CUDA 的使用非常适合加速该过程。有关 CUDA 的简要介绍，请参见此处。ZPrize 的结果令人鼓舞，速度提高了 2 倍以上。更好的是，每个解决方案都引入了不同的技巧和策略。下面你将找到一些解决方案的概述以及指向每个存储库的链接。

目录

• 目标
• Pippenger 算法
• 基本算法
• Speakspeak 的提交
• 6block 的提交
• Mike Voronov 和 Alex Kolganov 的提交
• MatterLabs 的提交
• Yrrid 的提交

目标

给定 $n \in \mathbb{N}$，椭圆曲线点 $P_1, \dots, P_n$ 和有限域中的标量 $k_1, \dots, k_n$，计算

$P = \sum_{i=1}^n k_i P_i$

其中求和是根据普通的椭圆曲线加法来理解的。

点的数量与表示计算所需的约束数量有关（可能大于 100,000,000）。例如，当我们想使用 Kate-Zaverucha-Goldberg (KZG) 承诺方案计算多项式 $a_0 + a_1x + a_2x^2 + ... a_dx^d$ 的承诺时，就会出现此计算。

Pippenger 算法

设 $\lambda$ 是存储标量所需的位数，设 $s$ 是 1 到 $\lambda$ 之间的整数。用 $\lceil \lambda/s \rceil$ 表示 $\lambda/s$ 的上限，即大于或等于 $\lambda/s$ 的最小整数。

将每个标量 $k_i$ 用 $2^s$ 进制表示

$ki = \sum{j=0}^{\lceil \lambda/s \rceil - 1} m_{i,j} (2^s)^j$，

其中 $0 \le m_{i,j} < 2^s$。那么

$\sum_{i=1}^n k_i Pi = \sum{i=1}^n \sumj m{i,j} 2^{sj} Pi = \sum{j=0}^{\lceil \lambda/s \rceil - 1} 2^{sj} (\sum{i=1}^n m{i,j} P_i)$。

让我们以不同的方式重写内部求和。对于每个 $1 \le m < 2^s$，我们可以对所有具有 $m_{i,j} = m$ 的内部求和项进行分组并编写

$\sum{i=1}^n m{i,j} Pi = \sum{m=1}^{2^s - 1} m (\sum{i:, m{i,j} = m} P_i)$

对于没有 $i$ 使得 $m{i,j} = m$ 的元素 $m$，我们将总和 $\sum{i:, m_{i,j} = m} P_i$ 解释为 $0$。最后一步称为 分桶。

综上所述，我们得到：

$\sum_{i=1}^n k_i Pi = \sum{j=0}^{\lceil \lambda/s \rceil - 1} 2^{sj} \sum{m=1}^{2^s - 1} m \sum{i:, m_{i,j} = m} P_i \qquad (1)$

Pippenger 算法包括从最里面的总和开始计算上面的总和：

(1) 对于每个 $j$ 和 $m$，计算 $B{j,m} := \sum{i:, m_{i,j} = m} P_i$。

(2) 对于每个 $j$，按如下方式计算 $G_j := \summ B{j,m}$。对于所有 $1 \le m < 2^s$，按索引的降序计算所有部分和

$S{j,m} = B{j,2^s - 1} + \dots + B_{j,m}$。

然后计算部分和 $S{j,1} + \dots + S{j,2^s - 1}$ 的总和。这等于我们想要的总和 $G_j$。

(3) 计算

$\sum_{j=0}^{\lceil \lambda/s \rceil - 1} 2^{sj} G_j$。

在伪代码中（摘自本文）：

基本算法

该实现在 algorithms/src/msm/variable_base/ 中。它特定于 BLS12-377 曲线。对于此曲线，我们有 $\lambda = 253$。

Aleo 使用 $s = 1$ 的 Pippenger 算法。方程 (1) 简化为

$\sum_{i=1}^n k_i Pi = \sum{j=0}^{\lambda - 1} 2^j \sum{i:, 1 = m{i,j}} P_i$，

其中 $m{i,j}$ 的定义与之前相同，但在此特定情况下，$m{i,j}$ 与 $k_i$ 的第 $j$ 位一致。

对于这种特定的 $s$ 选择，Pippenger 算法的步骤 1 是微不足道的，我们得到 $Gj = B{j,1}$。

并行化策略

CUDA 并行化仅用于按如下方式修改步骤 2。

(2) 目标是计算所有 $j$ 的 $Gj = \sum{i:, 1 = m_{i,j}} P_i$。为此，执行以下步骤。

(2.a) 首先，并行计算所有 $0 \le a < \lceil n/128 \rceil$ 和所有 $j$ 的

$G{j,a} = \sum{i:, 1 = m_{i,j}, 128a \le i < 128(a+1)} P_i$

这是使用 $\lambda * \lceil n/128 \rceil$ 个线程完成的。

(2.b) 然后，通过添加所有 $a$ 的 $G_{j,a}$ 来计算每个 $j$ 的 $G_j$。每个 $j$ 都有其线程，因此此步骤需要 $\lambda$ 步。每个线程添加 $\lceil n/128 \rceil$ 个椭圆曲线点。

一旦计算出所有 $G_j$，Pippenger 算法的其余部分将在 CPU 中执行，如上一节所述。

Speakspeak

这篇文章的标题是 "cuZK: Accelerating Zero-Knowledge Proof with A Faster Parallel Multi-Scalar Multiplication Algorithm on GPUs"，可以在这里找到。论文中描述的内容与 ZPrize 提交的实际实现之间存在差异。

论文中描述的并行化策略

此处的策略是更改 Pippenger 算法的步骤 1 和步骤 2，以利用 GPU 并行化。

我们使用 Pippenger 算法部分中引入的符号。设 $t$ 是要使用的线程数。

(1) 按如下方式计算 $B_{j,m}$。对于每个 $1 \le j < \lceil \lambda/s \rceil$：

(1.a) 使用所有 $t$ 个线程并行计算所有 $i$ 的 $m_{i,j}$。

(1.b) 对于每个 $0 \le l < t$，计算

$B{j,m,l} := \sum{i \text{ such that } m_{i,j} = m \equiv l \mod t} P_i$

为此使用所有 $t$ 个线程。

(1.c) 设 $M(j)$ 是具有椭圆曲线点条目的矩阵，使得 $M(j){m,l} = B{j,m,l}$。这是一个稀疏矩阵。计算 $B_{j,m} = M(j) \cdot 1_t$，其中 $1_t$ 是长度为 $t$ 且所有条目都等于 $1$ 的向量。这可以使用现有的稀疏矩阵-向量乘法的并行算法来完成。为此使用所有 $t$ 个线程。

(2) 按如下方式计算所有 $G_j$。对于所有 $0 \le j < \lceil \lambda/s \rceil$，使用每个 $t' := t / \lceil \lambda/s \rceil$ 个线程并行执行以下操作。

(2.a) 对于给定的 $j$，为了计算 $G_j = \summ B{j,m}$，将总和分成 $t'$ 个均匀块，并在其线程中分别计算每个块。也就是说，如果我们用 $\sigma = (2^s - 1) / t'$ 表示，对于每个 $0 \le \xi < \sigma$，计算

$\sum{m = \xi \sigma}^{(\xi+1)\sigma - 1} m B{j,m}$。

可以使用与 Pippenger 的步骤 2 中相同的部分和技巧来完成此操作。在这种情况下，需要一个额外的步骤，因为上面总和中的系数序列是 $\xi \sigma, \xi \sigma + 1, \dots$，而不是 $1, 2, 3, \dots$。但这很容易通过将 $(\xi \sigma - 1)$ 乘以最大部分和来解决。

(2.b) 添加上一步的所有块。结果是 $G_j$。

最后，按照 Pippenger 中的步骤 3 进行计算。

在伪代码中：

实现中的并行化策略

提交的实际代码中的并行化策略非常简单。没有使用稀疏矩阵乘法进行优化。但是，有一些有趣的事情需要注意。

• 仅支持 BLS12-377 曲线。
• 窗口大小 $s$ 选择为 $21$。由于 $\lambda = 253 = 12 21 + 1$，因此有 13 个窗口，其中一个窗口只有二进制标量。最后一个窗口的处理方式与其他 12 个窗口不同。考虑到 $253 = 11 23$，这是一个奇怪的选择
• 用于存储输入、结果以及所有中间部分结果的所有内存都在开始时分配。这需要大量内存。在 $2^{26}$ 个基点的情况下，仅存储所有窗口的标量就需要大约 3.2GB 的 GPU RAM。
• 在大多数情况下（前 12 个窗口），内核以 12 列和 $M$ 行的块网格启动，其中 $M$ 根据任务而变化。另一方面，块是一维的，大小为 32，因此 warp 和块重合。然后，每列处理与特定窗口相关的计算。
• 大量使用 cub library 进行排序、计算运行长度、计算累积和 和过滤 列表。
• 大部分代码实际上位于 sppark/msm 目录中。原始的 sppark/msm 代码已被修改。

该存储库包括代码主要部分的演练。这是一个摘要。

设 $n = 2^{26}$ 是基点的数量，设 $N$ 是 $2 \text{Cores} / (12 32)$，其中 Cores 是 GPU 的内核数。通常将使用 12 x $N$ 个 32 线程块的网格启动内核。$N$ 中的因子 2 在至少一步中使用二进制归约算法来将大小为 $12 N 32$ 的数组的所有点相加时才有意义。

对于 Pippenger 的步骤 (1)。

(1.a) 为了计算所有 $i, j$ 的 $m_{i,j}$，启动了一个具有 $N$ x 12 个 32 线程块的内核。所有标量都被分区，并且每个线程负责计算其分区中系数 $ki$ 的所有 $j$ 的 $m{i,j}$。分区的大小约为 $n / (32 * N)$。

(1.b) 依次对于每个窗口 $j$，使用 cub::DeviceRadixSort::SortPairs 对标量集合 $m{i,j}$ 进行排序。让我们将排序后的窗口 $j$ 的第 $i$ 个标量表示为 $m'{i,j}$。

(1.c) 依次对于每个窗口 $j$，使用 cub::DeviceRunLengthEncode::Encode 在先前排序的标量 $m'_{i,j}$ 上计算窗口中每个标量 $1 \le m < 2^s$ 的出现次数。

(1.d) 出于下一步所需的技术原因，使用 cub::DeviceScan::InclusiveSum 计算出现次数的累积和。

(1.e) 启动一个内核来计算存储桶。该内核获得一个大小为 $N$ x 12 个 32 线程块的网格。总共 12 列的第 $j$ 列处理窗口 $j$ 的存储桶。索引范围 1 到 $n$ 平均分为子范围，并且每个线程处理与其范围内的 $i$ 的 $m'_{i,j}$ 对应的唯一标量的存储桶。范围略有扩大和缩小，以便线程获得不重叠的标量集。

这结束了所有 $0 \le j < 12$ 和 $1 \le m < 2^s$ 的存储桶 $B_{j,m}$ 的计算。

对于 Pippenger 的步骤 (2)。

(2.a) 在 $N$ x 12 个 32 线程块的网格上启动一个内核。每个线程计算总和 $G_j = \summ B{j,m}$ 的均匀块，就像论文中描述的那样。与之前一样，网格的每一列处理一个不同的窗口。

(2.b) 对于每个窗口 $j$，使用二进制归约算法将其块相加。这有效地计算了所有 $0 \le j < 12$ 的 $G_j$。

然后，在 CPU 中执行 Pippenger 的步骤 (3)。

6block 提交

此解决方案的主要贡献是步骤 (1) 和步骤 (2) 的不同方法。如果我们暂时忘记 GPU 并行化，则两个步骤都按如下方式一步执行。

(1') 对于每个窗口 $j$，首先对所有标量 $m{i,j}$ 进行排序。分别用 $m'{i,j}$ 和 $P'_i$ 表示标量和点的排序列表。对于从 $n$ 到 $1$ 的每个 $i$，其中 $n$ 是基点的数量，计算

$t_{i-1} := t_i + P'i$ $s{i-1} := (m'{i,j} - m'{i-1,j}) t_i + s_i$

其中 $t_n = P'_n$ 和 $s_n = O$。那么 $Gj$ 等于 $m'{0,j} t_0 + s_0$。其余方法与 Pippenger 中的相同。

并行化策略。

设 $n$ 是基点的数量 $2^{26}$。使用的窗口大小为 $21$。将 253 位分组为 11 个大小为 21 的窗口和一个大小为 22 的附加窗口。

(1'.a) 在最多 128 个线程的一维块的一维网格上启动一个内核。块的精确大小是使用 CUDA 占用率计算函数 cudaOccupancyMaxPotentialBlockSize 计算的。线程总数等于点的数量。然后，每个线程负责计算单个标量 $ki$ 的所有 $m{i,j}$。

以下步骤按顺序为每个窗口 $j$ 执行。

(1'.b) 在固定 $j$ 的情况下，对列表 $(m{i,j}){i=1}^n$ 进行排序 _。为此使用了 cub 函数 cub::DeviceRadixSort::SortPairs。此函数对键值对进行排序。在这种情况下，排序的对是 $(m{i,j}, i)$，以跟踪排序列表中哪个基点对应于哪个标量。让我们分别用 $m'{i,j}$ 和 $P'_i$ 表示标量和点的排序列表。

(1'.c) 然后，在最多 128 个线程的一维块的一维网格上启动一个内核。同样，块的精确大小是使用 cudaOccupancyMaxPotentialBlockSize 计算的。范围 1 到 $n$ 平均分配，并且每个线程负责计算 $i$ 在其范围内的总和 $\sum m'_{i,j} P'_i$。这是通过计算如上所述的 $s_0$ 和 $t0$ 来完成的。这将产生一定数量的部分结果（每个线程一个），我们在此处将其表示为 $B{k,j}$。所有 $k$ 的所有这些元素的总和等于 $G_j$。

(1'.d) 将所有 $k$ 的结果 $B_{k,j}$ 复制到 CPU 并按顺序添加以获得 $G_j$。然后，将其加倍 $2^{1}$ 次以获得 $2^{21} Gj$（最后一个 $2^{22} G{12}$ 的窗口除外）。发生这种情况时，在 GPU 中处理后续窗口的步骤 (1'.b) 和 (1'.c)。

处理完所有窗口后，在 CPU 中执行最终总和。

关于此解决方案的一些有趣的事项。

• 代码非常清晰。它使用 OOP 和许多 c++11 功能和标准库。
• 在步骤 (1'c) 中启动的内核的简单实现可能会严重受到 warp 差异的影响。这是因为在构建 $t_i$ 和 $si$ 中存在大量分支。例如，如果一个 $m'{i,j}$ 等于 $m'_{i-1,j}$，则无需执行任何操作来计算 $si$。为了克服此问题，每个线程都使用其遇到的非零差异 $m{i,j} - m_{i-1,j}$ 填充最大大小为 10 的缓冲区。所有椭圆曲线运算都将推迟到 warp 中的一个线程填满其缓冲区为止。此时，warp 中的所有线程都会刷新其挂起的椭圆曲线运算。为此，使用了 warp 投票函数 __any_sync（请参见此处）。

Mike Voronov 和 Alex Kolganov

主要地，此提交通过使用有符号标量来减少 Pippenger 步骤 (1) 中的 bucket 数来改进基线（更多详细信息如下）。尽管它也声称使用仔细的线程平铺来并行计算 Pippenger 算法步骤 (2) 中的部分和，但几乎没有相关文档。

主要贡献在于步骤 (1)。此解决方案使用窗口大小 $s = 23$。

(1.a) 为了计算 $m_{i,j}$，启动了一个大小为 256 的一维块的一维网格。线程总数等于基点的数量，即 $2^{26}$。每个线程负责计算单个 $ki$ 的子标量 $m{i,j}$。由于窗口大小为 $2^3$，因此所有子标量 $m{i,j}$ 都满足 $0 \le m{i,j} < 2^3$。如果子标量 $m{i,j}$ 大于 $2^2$，则从该子标量中减去 $2^3$，将其减小到范围 $-2^2 \le m{i,j} < 0$，并将 $2^3$ 结转到下一个窗口。负子标量的符号转移到基点，因为否定椭圆曲线点的成本很低。因此，它们最终得到范围 $0 \le m_{i,j} < 2^2$ 中的子标量，并且如果最后一个窗口需要结转标量，则可能会得到一个额外的窗口。

然后将窗口分为两组：具有奇数和偶数索引的窗口。窗口以异步方式处理，并且可以处理两个以上的窗口，config.numAsync 变量管理流计数。但是对于 A40，两个流足以利用所有计算资源。唯一的例外是最后一个窗口，由于它是上一步的溢出窗口，因此在 CPU 中单独处理该窗口，因此该窗口小得多。每个组都获取其线程流，并按顺序遍历其窗口以计算存储桶和最终窗口总和 $G_j$，如下所示。

(1.b) 对于每个窗口 $j$，它对子标量 $m_{i,j}$ 进行排序，并预先计算排序列表中每个子标量出现的开始和结束索引，以及每个子标量的出现次数。这是使用 thrust 库中的 thrust::sort_by_key 和 thrust::inclusive_scan 函数完成的。然后，它启动一个具有大小为 32 的一维块的一维网格的内核，以使用上述预先计算的信息来计算存储桶。

(2.a) 所有窗口都在不同的流中并行计算（使用了两个流，但是可以使用更多流，具体取决于 GPU 内存）。

(2.b) 然后对存储桶进行排序，以便首先运行具有最多点的存储桶。这使 GPU warp 可以运行收敛的工作负载并最大程度地减少尾部效应。此解决方案编写自定义算法来实现此目的。

然后，在 CPU 中执行 Pippenger 的步骤 (3)。

从此解决方案中要注意的其他事项。

• thrust 库的使用。
• 它始终使用大小为 256 或 32 的一维块的一维网格。
• 在步骤 (1.b) 中，使用两个流来并行化均匀和奇数窗口。这是因为两个流足以利用所有计算资源。
• 看起来它是在内存受限的 GPU 卡上开发和运行的。在步骤 (1.b) 中，每个组都重用已分配的数组来存储存储桶。它还重用已分配的数组，用于计算子标量的排序列表以及与每个子标量的第一次和最后一次出现相关联的索引的不同不相关的中间步骤。

未奏效的经过测试的想法。来自提交的 README：

• NAF - 有符号标量更好，因为 NAF 将存储桶的数量减少 2/3 倍，但有符号标量减少 2 倍
• NAF + 有符号标量 - 此变体的主要缺点是在第二步中，总和的计数是两倍
• Karatsuba 乘法 + Barrett 归约 - 事实证明，CIOS 基线 Montgomery 更好
• 仿射求和 + Montgomery 技巧 - 事实证明，速度比基线求和慢

MatterLabs

此解决方案预先计算每个输入点 $P_i$ 和 $j = 0, 1, 2, 3$ 的 $2^{69j} P_i$（这与文档中所述的不同）。这些都是 EC 中的点。我们可以通过以下方式重写第一个总和：

$\sum{i=1}^n \sum{j=0}^3 k_{ij} 2^{69j} P_i$

其中每个 $k_{ij} < 2^{69}$。

由于每个 $2^{69j} Pi$ 都属于 EC 并且已计算，因此我们可以将总和重写为 $\sum{m=1}^{3n} k_m P_m$（对于其他 $k$ 和其他 $P$），其中每个 $k_m < 2^{69}$。这使我们可以将每个 $k_m$ 分成三个 23 位窗口用于 Pippenger 算法。

Arkworks 库用于在测试中表示有限域、椭圆曲线和大整数。但是，对于 MSM 算法本身，这些结构由作者实现。它们使用优化的操作版本在运行设备代码时在 GPU 上运行它们。

开发了一个名为 Bellman-CUDA 的新库。它用于对有限域执行操作、排序（使用 CUB 实用程序）、运行长度编码（也使用 CUB），等等，从而利用 GPU。此目的可能是为了将来用更高效的算法替换对 CUB 的调用。

这些窗口以较小的部分处理。首先处理所有窗口的第一个块，然后处理所有窗口的第二个块，依此类推。这允许在其他标量部分仍从主机异步传输到设备内存时进行处理。

对于每个窗口块，并行生成每个存储桶的元组索引：$(i,j)$，其中 $i$ 是存储桶的系数，$j$ 是该存储桶中的 EC 点。根据第一个分量对它们进行排序（并行排序），以便我们拥有要在每个存储桶中求和的 EC 点。然后，根据存储桶拥有的点数对它们进行长度编码和排序（也并行排序）。这样，将首先处理具有更多点的存储桶，以实现 GPU 硬件的有效使用。之后，生成一个偏移量列表（使用 CUB 实现的并行算法来计算互斥总和）以了解每个存储桶的开始和结束位置。例如：

$2P_1 + 5P_2 + 4P_3 + 1P_4 + 2P_5 \rightarrow (2,1), (5,2), (4,3), (1,4), (2,5) \rightarrow [(1,4), (2,1), (2,5), (4,3), (5,2)], [1,2,4,5], [1,2,1,1] \rightarrow [4,1,5,3,2], [0,1,3,4,5], [1,2,4,5], [1,2,1,1]$

然后，这些存储桶并行聚合。

FF 和 EC 例程已得到优化：

• 基于 Montgomery 的乘法
• 最大程度地减少了 FF 操作中的校正步骤
• 对 EC 点累加器使用 XYZZ 表示形式
• 使用快速平方

流和内存管理

创建以下流：

• stream
• stream_copy_scalars
• stream_copy_bases
• stream_copy_finished
• stream_sort_a
• stream_sort_b

第一个是主流。诸如 initialize_buckets、compute_bucket_indexes、run_length_encode、exclusive_sum 和 sort_pairs 之类的内核在该流中运行。

stream_copy_scalars 等待 event_scalars_free。

stream_copy_scalars 处理标量的异步复制并排队 event_scalars_loaded。

stream_copy_bases 等待 event_scalars_loaded 和 event_bases_free。此流还处理基数的异步复制并排队 event_bases_loaded。

stream 等待 event_scalars_loaded，处理内核 compute_bucket_indexes，并排队 event_scalars_free。

stream 处理索引的排序以及索引和运行长度的异步内存分配，以及 exclusive_sum 内核和偏移量的内存分配。

stream 将 event_sort_inputs_ready 排队。stream_sort_a 和 stream_sort_b 等待该事件以处理 GPU 上的对排序。

stream_sort_a 将 event_sort_a 排队，stream_sort_b 将 event_sort_b 排队。stream 在处理聚合存储桶的内核之前，等待该事件以及 event_bases_loaded。stream 将基的（异步）内存释放进行排队。

在窗口块处理的最后一个循环中，stream_copy_finished 等待 event_scalars_loaded 和 event_bases_loaded。

内存被释放，并且流（除了处理存储桶减少和窗口拆分内核的 stream 之外）被销毁。

存储桶聚合算法

在计算出每个存储桶后，使用此算法，并且该算法是聚合存储桶的并行化策略的基础。它是 Pippenger 中经典的部分和技巧的替代方法。在下面的内容中，我们假设已经计算出每个存储桶，并且剩下的问题是将每个窗口中的所有点相加。

符号

让我们确定一些符号。设 $W = (B0, \dots, B{2^b - 1})$ 是 $2^b$ 个椭圆曲线点的元组。让我们将这样的元组称为 $b$ 位窗口。对于每个窗口 $W$，我们都关联一个椭圆曲线点 $P_W$，其定义为

$P_W := B_1 + 2B2 + \dots + (2^b - 1) B{2^b - 1}$

我们将长度相同的 $m$ 个此类窗口 $C = (W0, \dots, W{m-1})$（长度为 $2^b$）的元组称为 窗口配置。我们说窗口配置的形状为 $(m,b)$。每个窗口配置都有一个关联的椭圆曲线点，由下式定义 PC:=PW0+2bPW1+22bPW2+⋯+2(m−1)bPWm−1PC:=PW0+2bPW1+22bPW2+⋯+2(m−1)bPWm−1

在 MSM 的上下文中，每个 $B_i$ 都是一个存储桶，而 $PC$ 是所需的最终结果。

减少过程

让我们假设每个桶都已经计算完毕，并令 $C$ 为对应的窗口配置。MatterLabs 的解决方案实现了一种算法，通过迭代地将形状为 $(m,b)$ 的窗口配置 $Ci$ 缩减为形状为 $(2m, \lceil b/2 \rceil)$ 的另一个窗口配置 $C{i+1}$ 来获得 $PC$。在每一步中，点 $PCi$ 不一定等于 $PC{i+1}$，但可以通过平移一些标量从 $C_{i+1}$ 获得。详见下文。该过程从形状为 $(3,23)$ 的配置 $C$ 开始，以形状为 $(96,1)$ 的配置 $D$ 结束。此时，$PC$ 从 $D$ 计算得出。

窗口分割

缩减包括分割配置的每个窗口。让我们描述一下单个 $b$ 位窗口 $W$ 的分割过程。我们从中构造两个新的 $\lceil b/2 \rceil$ 位窗口 $\hat{W}_0$ 和 $\hat{W}_1$，使得

PW=P^W0+2⌈b/2⌉P^W1.PW=PW^0+2⌈b/2⌉PW^1.

$$ PW = P^{\hat{W}_0} + 2^{\lceil b/2 \rceil} P^{\hat{W}_1}. $$

这个构造背后的想法如下。$W$ 的每个分量都具有 $B_r$ 的形式，其中 $0 \leq r < 2^b$。我们可以写成 $r = a + b2^k$，其中 $0 \leq a, b < 2^k$。然后，$B_r$ 被放入两个新的桶中，即窗口 $\hat{W}_0$ 的第 $a$ 个分量和窗口 $\hat{W}_1$ 的第 $b$ 个分量。

情况 $b$ 偶数：

写成 $b = 2k$。令 $W$ 为一个 $b$ 位窗口。定义新的 $k$ 位窗口 $\hat{W}_0$ 和 $\hat{W}_1$ 如下。

用 $(B0, \dots, B{2^b-1})$ 表示 $W$ 的分量。那么

^W0:=(2k−1∑i=0Bi2k,2k−1∑i=0Bi2k+1,…,2k−1∑i=0Bi2k+2k−1),^W1:=(2k−1∑i=0Bi,2k−1∑i=0Bi+2k,…,2k−1∑i=0Bi+(2k−1)2k)W^0:=(∑i=02k−1Bi2k,∑i=02k−1Bi2k+1,…,∑i=02k−1Bi2k+2k−1),W^1:=(∑i=02k−1Bi,∑i=02k−1Bi+2k,…,∑i=02k−1Bi+(2k−1)2k)

$$
\hat{W}_0 := \left( \sum_{i=0}^{2^{k-1}} B_{i2^k}, \sum_{i=0}^{2^{k-1}} B_{i2^k+1}, \dots, \sum_{i=0}^{2^{k-1}} B_{i2^k+2^{k-1}} \right),
$$
$$
\hat{W}_1 := \left( \sum_{i=0}^{2^{k-1}} B_i, \sum_{i=0}^{2^{k-1}} B_{i+2^k}, \dots, \sum_{i=0}^{2^{k-1}} B_{i+(2^{k-1})2^k} \right)
$$

情况 $b$ 奇数：

让我们写成 $b = 2k-1$。这种情况与上述情况类似。和以前一样，令 $W$ 为一个 $b$ 位窗口。$\hat{W}_0$ 和 $\hat{W}_1$ 的定义与之前遵循相同的逻辑。

但这里有一个陷阱。如果 $r$ 使得 $0 \leq r < 2^b$，并且我们写成 $r = a + b2^k$，其中 $0 \leq a, b < 2^k$，那么 $b$ 必然最多为 $2^{k-1}$。因此，$\hat{W}_1$ 的后半部分坐标将为空。这是因为 $W_n$ 的桶 $B_r$ 都不会分配给这些坐标。因此我们得到

^W0:=(2k−1−1∑i=0Bi2k,2k−1−1∑i=0Bi2k+1,…,2k−1−1∑i=0Bi2k+2k−1),^W1:=(2k−1∑i=0Bi,2k−1∑i=0Bi+2k,…,2k−1∑i=0Bi+(2k−1−1)2k,O,…,O).W^0:=(∑i=02k−1−1Bi2k,∑i=02k−1−1Bi2k+1,…,∑i=02k−1−1Bi2k+2k−1),W^1:=(∑i=02k−1Bi,∑i=02k−1Bi+2k,…,∑i=02k−1Bi+(2k−1−1)2k,O,…,O).

$$
\hat{W}_0 := \left( \sum_{i=0}^{2^{k-1}-1} B_{i2^k}, \sum_{i=0}^{2^{k-1}-1} B_{i2^k+1}, \dots, \sum_{i=0}^{2^{k-1}-1} B_{i2^k+2^{k-1}} \right),
$$
$$
\hat{W}_1 := \left( \sum_{i=0}^{2^{k-1}} B_i, \sum_{i=0}^{2^{k-1}} B_{i+2^k}, \dots, \sum_{i=0}^{2^{k-1}} B_{i+(2^{k-1}-1)2^k}, O, \dots, O \right).
$$

在上面的定义中，有 $2^{k-1}$ 个坐标具有条目 $O$，即无穷远点。

窗口配置的缩减和系数偏移

在配置 $C$ 的每个窗口上执行上述过程，我们得到所需形状的新配置 $D$。我们不会总是得到 $PC=PD$。

令 $C = (W_0, W_1, \dots, W_n)$ 为一个形状为 $(m,b)$ 的窗口配置。对于每个窗口 $Wn$，令 $\hat{W}{2n}$ 和 $\hat{W}_{2n+1}$ 为从分割 $W_n$ 获得的两个 $\lceil b/2 \rceil$ 位窗口。令 $D = (\hat{W}_0, \hat{W}1, \dots, \hat{W}{2m-1})$。这是一个形状为 $(2m, \lceil b/2 \rceil)$ 的窗口配置。

如果 $b$ 是偶数，那么很容易看出 $PC=PD$。

如果 $b$ 是奇数，那么通常情况下 $PC$ 与 $PD$ 不同。例如，$PC$ 的前 2 项是 $W_0 + 2^b W_1$。另一方面，$PD$ 的前四项是 $\hat{W}_0 + 2^k \hat{W}_1 + 2^{2k} \hat{W}_2 + 2^{3k} \hat{W}_3$。这等于 $W_0 + 2^{2k} W_1 = W_0 + 2^{b+1} W_1$。因此，$PD$ 中 $W_1$ 的系数有一个额外的因子 2。

尽管如此，$PC$ 等于

P^W0+2k−f1P^W1+22k−f2P^W2+⋯+2(2m−1)k−f2m−1P^W2m−1,PW^0+2k−f1PW^1+22k−f2PW^2+⋯+2(2m−1)k−f2m−1PW^2m−1,

$$
P_{\hat{W}_0} + 2^{k-f_1} P_{\hat{W}_1} + 2^{2k-f_2} P_{\hat{W}_2} + \cdots + 2^{(2m-1)k-f_{2m-1}} P_{\hat{W}_{2m-1}},
$$

其中 $f_i = \lfloor i/2 \rfloor$。我们称这些为系数偏移。

通常，我们可以定义 $f_i$ 对于所有 $i$ 都是 $0$，如果 $b$ 是偶数；对于所有 $i$，如果 $b$ 是奇数，则 $f_i = \lfloor i/2 \rfloor$。

算法

我们从一个形状为 $(m,b) = (3,23)$ 的窗口配置 $C_0$ 开始。归纳地，对于每个 $i$，对 $Ci$ 执行缩减步骤以获得一个新的窗口配置 $C{i+1}$，并累积系数偏移。经过 4 步后，我们获得形状为 $(96,1)$ 的 $C_5$ 和累积的系数偏移 $f_i$。

从 $C_5$ 和 $f_i$ 我们可以计算 $PC_0$。

并行化策略

当将形状为 $(m,b)$ 的窗口配置拆分为形状为 $(2m,k)$ 的窗口配置时，其中 $k = \lceil b/2 \rceil$，每个新桶都是 $2^k$ 个元素（或者当 $b$ 为奇数时，是 $2^{k-1}$ 个元素）的和。为了计算这些，启动以下内核。

1. 首先，启动一个具有 $2m2^{k+l}$ 个线程的内核，其中 $l \leq \lfloor b/2 \rfloor$。每个新桶的和的 $2^k$ 项被分成 $2^l$ 个均匀的组。然后每个线程计算组中各项的和。这些部分和是顺序计算的。
2. 启动第二个具有 $2m2^k$ 个线程的内核。每个线程负责一个桶。它使用二元缩减算法通过添加先前内核获得的 $2^l$ 个部分和来计算它。

Yrrid

此解决方案为每个输入点 $P_i$ 预计算 $2^{2 \cdot 23j} P_i$，其中 $j=1,...,6$。 这些都是 EC 中的点。 我们可以这样重写第一个求和：

n∑i=15∑j=0kij22⋅23jPi∑i=1n∑j=05kij22⋅23jPi

$$
\sum_{i=1}^n \sum_{j=0}^5 k_{ij} 2^{2 \cdot 23j} P_i
$$

其中每个 $k_{ij} < 2^{2 \cdot 23}$。

由于每个 $2^{46j} P_i$ 属于 EC 并且已经被计算，我们可以将总和重写为 ∑6nm=1kmPm∑m=16nkmPm（对于其他 $k$ 和其他 $P$），其中每个 $k_m < 2^{46}$。 这允许我们

$$
\sum_{m=1}^{6n} k_m P_m
$$

将每个 $k_m$ 分成两个 23 位窗口，用于 Pippenger 算法。

该算法使用的另一个优化是：窗口值具有一个符号位和一个 22 位标量值。 如果标量很大，我们可以否定该点并将标量更改为 $s' = m-s$，其中 $m$ 是域的阶数。 新的标量 $s'$ 将具有一个高位清除位。 这是因为 $s'(-P_i) = (m-s)(-P_i) = -s - P_i = sP_i$。

然后对桶进行排序，以便首先运行具有最多点的桶。这允许 GPU 线程束 运行收敛的工作负载并最大限度地减少尾部效应。 该解决方案编写自定义算法来实现桶排序，而不是使用 CUB 库。

桶总和使用 XYZZ EC 表示并行计算（为每个桶分配一个线程）。 可以在此处找到此曲线表示的运算。

FF 和 EC 例程已得到优化：

• 基于 Montgomery 的乘法
• 最小化 FF 运算中的校正步骤
• 对 EC 点累加器使用 XYZZ 表示
• 使用快速平方

总结

多标量乘法 (MSM) 是许多证明系统中的关键操作之一，例如具有 Kate 多项式承诺方案的 Marlin 或 Plonk。 由于操作的性质，我们可以利用 GPU 来减少其计算时间。 ZPrize 竞赛旨在改进 BLS12-377 曲线的 $2^{26}$ 点的 MSM 的当前 5.86 秒的基线。 有 6 个不同的提案，每个提案都具有基于 Pippenger 算法的独特功能：优化窗口大小、预计算某些点（用内存换速度）、用于椭圆曲线加法的不同坐标系、自同态、并行缩减算法、点否定、整数的非邻接形式、更好的有限域算术。 最佳解决方案实现了 2.52 秒（加速 2.3 倍），但我们认为仍有更多优化空间。 我们会低于 1 秒吗？ 也许你有答案...

• 原文链接： blog.lambdaclass.com/eye...
• 登链社区 AI 助手，为大家转译优秀英文文章，如有翻译不通的地方，还请包涵～