「深入 ZKP 底层系列」

• 背景
  当下，ZKP（Zero-Knowledge Proof） 已经在加密货币、数据隐私保护、安全认证等多个领域发挥着至关重要的作用。
  Plancker Langlands 小组现也已启动 ZKP 研究，专注于通过硬件加速 MSM 与 NTT 等运算，提高 ZKP 中 proof 的生成速度。
  ZKP 的底层协议是零知识证明技术的核心和灵魂，在本系列中，我们将一起带你深入展开探索。
• 文章概述
  椭圆曲线密码学（ECC）在 PLONK 中有着重要的应用，而 ECC 上的一个重要运算即 MSM，其基础是由点乘进行构造的，更进一步可分解为点加和倍点运算，再进一步分解为：模乘、模平方、模加减、模逆这些模运算。其中，模逆运算是最为昂贵的计算，本篇文章对模逆运算上的巴雷特约减（Barrett Reduction）改进版进行了说明和公式推导。

引言

MSM 的性能直接影响了 PLONK 系统的整体性能。PLONK 是一种特定类型的 ZK-SNARK，被认为是一种比传统 ZK-SNARKs 更灵活和高效的构造，它的效率很大程度上依赖于 MSM 操作的效率。

而在实现椭圆曲线加密算法时，MSM 的效率又在很大程度上依赖于模乘运算的效率，因为模乘运算是执行标量乘法的基础，所以提高这些基础数学运算的效率对于性能的提升至关重要。

本文是 Langlands ZK  硬件加速小组近期论文阅读和分享 pipeMSM 的公式推导梳理，主要涉及一个优化算法 —— Barrett Reduction 算法（具体见 Barrett 1987年论文^[1]）。

我们当前没有发现相关的代码实现（比如在哪个相关方向或分支上），如果有读者知道这些信息，欢迎提供给我们，让我们一起 dive deeply！

正文

更新版论文，作者 Yuval Domb，来自 Ingonyama

• 2022 年发布 Fast Modular Multiplication (模乘论文第一篇)^[2]
• 更新版本：https://github.com/ingonyama-zk/papers/blob/main/multi_precision_fast_mod_mul.pdf

论文的算法推导

理论：原始问题和转化待求的问题

原始问题「公式 (1)」:

其中 ， 为素数，并利用 标准。

等价为「公式 (2）」：

其中， （是一个整数），使得 。

注：公式（2）可以转换成 ，其核心思想是不作除法来实现计算——大约估出一个满足条件的 ，将 代入后获得 ， 既是余数，也是取模之后的结果。

因为公式（1）和（2），所以求  的公式可以转化为：

注：普通的 mod 是用除法来做，但是硬件操作特别昂贵，本算法用位移和加法来代替除法

算法实现精简整合

理论部分：

首先，估算一个 ，让它和 的值非常接近。非常接近的意思是指误差，即  是一个常数(更准确地说 是一个上限)，即 ；

之后，取 ，（ 可以提前算出，并且 与 无关），再去计算 （  是估计出来的数值）。

实际的误差可以用    来表示， 这个值的上限是 4，也就是意味着估算值与实际值之间的差为 。

接下来就是对 具体的证明 ：

证明过程中会涉及一些表达式，为了容易理解，以下是举例：

这里

（n=4）

那么，

因为 ，

所以 的近似值

（ 向下取整，所以取整的数值一定小于 ），

（同样是向下取整，所以 ）

则

我们要计算 ，而 ，所以 可以用 n 位（bit）表示，那么看起来我们只需要取    和 就可以了。

根据上面提供的举例来看，如果 是 1110， 就是 11110， = 101100。假设误差是 ，而你不能只取低的次位进行计算，此时就需要额外的两位进行计算，即 n+2 位，也就是 的后 n+2 位减去 的后 n+2 位。

如何计算？（ 是 位表示，但是只需要 位的后 n+2 位）

，s 用 n 位表示；

，所以 用 n 位表示即可。

 如何计算？

观摩 ，

可以得到：

在椭圆曲线上， s 一定是大于 的

即需 n+1 位来表示 m 的最高位为 1。

本来我们只能用 n 位和 n+1 位去算 ，但是有一个优化，以下是优化的具体形式：

算法实现部分

流程图详解：

• 注1: 首先拆解

• 注2:    取

• 其中 m 的值介于 与 之间，故其最高位始终为 1，所以 可以用 n 位表示;
• 计算 是为了得到   的实际数值

• 注3: 的计算结果再加一次    ，从而完整的计算下列公式，去获得

• 注4: 之后就是算 ，取 n 位即可，因为我们最后要求解的值 满足该关系式 ，故 是一定小于 的，所以计算出

• 注5: 取 的 n+2 位 （因为这里可能存在误差，所以还是需要多留两位进行计算）

• 之后出来的误差去和 比，比 大就 ，从而计算出

论文里的推导逻辑

这部分不感兴趣的其实可以跳过，是论文里概括性的推导逻辑，比较简略的过程，和上文有一定的重复。

其中 ，s 为素数，并利用 标准。

等价为：

其中， （是一个整数），使得 。

注：（2）的公式可以转换成  ，（2）的核心思想是不作除法来实现计算，大约估出一个满足条件的 ，将 代入后获得 ， 既是余数，也是取模之后的结果

因为公式（1）和（2），所以求  的公式可以转化为：

注：普通的 mod 是用除法来做，但是除法的硬件开销非常高，本算法用位移和加法来代替除法

将所有变量以 d 进制来表示，其中 内的每个元素都以 n 个位来表示，有：

注：（3）里的 表示所有小于 s 整数的最大的位数， 采用向上取整，比如 ，则 ；

接下来，简单点，令 d = 2 ，所有元素以二进制来表示。

尽管本文重点关注模素数乘法（ modulo-prime multiplication），但可将其推广到任意 运算，其中 ， 可为素数或非素数的任意值。

3.1 假设 为近似已知

假设 为近似已知，将其近似值表示为   ，使得：

注：相当于  

其中   为一个已知的常量。

若 ，则有：

注：其中[ ]中括号内的值表示了位的位置（bits location）和规模（sizes）。 和 都是宽度为 n 的数，两个相乘后用 2n 宽度的数就能表示。后面类似。

注意，当 时，可推测余数 r 最大长度为  n bits，使得等式（5）中右侧值的剩余最高有效位（msb (most-significant bits）必须为 0。

通过简单的 bit 操作，可以 long addition 表示为：

其中，上横杠表示的是 bit-inversion 运算符，右侧上的 「1」 表示为初始进位位（carry bit）。

不过，对上图的长加法仔细观察可知，仅需要   和 来完成该计算，从而可节约近一半的计算量。最终的加法（adder） 为一个固定宽度加法（fixed width adder） —— 即 。这意味着可忽略最高有效位 （ms bits（msb））的任何溢出。可将其等价为 a fixed-width subtractor ——即 n ，可将其结果看成是无符号整数。

将生成以上乘积的乘数表示为 ，其中 是指该整个产品的 n 个 lsb（least-significant bits）。 和   都可通过   来生成。

此外，若 为常量， 可通过一个常量   相乘来生成。

当   时：

推导如下：

此时，用于表示等式（5）中右侧值所需的位数为：

推导如下：

注：

• （4）中因为 ，所以 一定小于 ， 向上取整比它大的数，所以就用 替换
• （5）中因为 ，所以 ，所以（4）中的 = 换成了 ≤

因此，若 ，则仅需要额外再增加 1 个位（bit）来表示。

3.2 用巴雷特约减算法求 的近似值

采用巴雷特的模块约减算法对 求近似值为：

注：（8）中 m 和 k 的关系，可以理解为：k 越大， m 的分辨率越高

其中：

关于 k 和 n 为什么要这么取

这正是巴雷特创新的点。假如做模乘的时候 s 不去模，直接相乘得到值，这个值的位数就是（k+n）位 ，通过这个公式（9）可以算出 m 的值。

推导如下：

m(k) 是一个变量为 k 的函数，最多有 k+1 位，为公式（8）的下界近似器（ lower-bound approximator)。对于有限的 k 值，该误差可被表示为：

利用前面的条件，可推导如下：

详细解释:
首先，将 m(k) 代入（10）中；
其次，通分
∵ 公式（9）中 ,  最小的整数是当 k=0 时，为 2，
而 是向下取整，所以只能取 1 ，，
∴e(k) ＜
假设当 ，则
e(k)=0, 则

其中，可检查二进制表示的左右项的最大差异来派生出该上限，之后将 代入进来，从而 的误差计算为：

若  ，则该误差最大为 1 。

3.3 参数选择 以及 error bounding

『这部分的目的是为了求出真实的 和 的差』

根据 （8）对于 的定义，以及（11）的关系，当 k = n（即 ），则对 的近似值为：

其中乘法为：

且上述公式误差的计算遵循（11）。

分两个阶段来实现以上相乘：

1）首先，假设有 ，按如下方式计算相乘：

其中最右侧项最大可能值显而易见地被确定为 2。

注：在最新的论文中，，所以 的误差为 3，再取下个整数，误差会扩大 1，即

推导如下：

注：这里 一定是＞0的，先推出 ，再转化成与 2s 的关系

2）从而 的近似变为：

注：在最新的论文中，此处的误差为 4，

推导如下：

其中最里侧相乘为 ，最外侧的常量相乘为: ，近似误差的上限是 （13）和 （15） 中最右边项的总和。

注：

• m 非常接近 ，那么两者相除就非常接近 1，再乘以一个 n 位的数，基本 n+1 位就足够表示了，因为
• 尽管如此，还是要看有没有超过上限。所以算法实现里是 n+2 ，如下图所示

3.4 总体算法

以下为硬件优化的模块化乘法器结构图(基于最新的论文)（这个结构图实现了 ， 为上面的 ），假定了 和 为已知的常量，使用 来表示  的近似值。

流程图的解释：

• 是一个 2n 的数，取了最高有效位和最低有效位，最高有效位表现为 ，用在了 上，用来计算出 ；
• 里，又把 用了一遍，只是计算 的时候用的是它的最低有效位，  的最低有效位因为误差的原因，所以多取了两位
• 这个流程图最终实现的是 ，而  也就是公式（16）:

这个流程图的大致实现思路：

1. 把 的结果拆成了高位和低位两部分
2. 用公式（16）去替换
3. 再用 得到结果
4. 由于有误差，所以再减去

有几个点要注意：

• 在上面第三步所得出的结果是比实际值更大的，因为 有误差；
• 假设 没有误差，那么因为 ，就可以算出 的值了，但实际上 是有误差的，所以要算真正的 的值，根据 来看，实际的 的值可能减少了几个 ；
• 具体要算多少个 ，可以将 的值和 对比，如果值＞，那么就减去 （可以理解为 的误差是几，就减多少次 ；一般来说减一次 即可）
• 这个图里面的计算省略了除法计算，这也是他的优化：做了三次乘法，两次位移，巧妙地用了 n 位的最高有限位

对 的计算也做了优化——做了拆分：

• 拆到处理器正好能处理的位数，从而加快计算速度
• 核心思想： 优化，取模优化为不用除法，底层的 优化，最终把 MSM 优化。

优化总结

新版论文相较于初版论文，新版优化点：

• 纠正了 误差的上限或下限为 4 而非 3，猜想的推导过程具体可见论文，这里不做详细解释；

• 证明  位足以表示 （优化了 范围的取值，n 位就够了，之前是 n+1 位）;

• 稍微修改最高位乘法器（MSB multiplier），为多精度方案做准备。

相关论文链接：