写在前面
本篇Blog为了方便理解和阅读，修改了原文的一些符号和记法
Abstract
d-MSM算法主要用于计算下列等式，其中$a_1,\dots,a_d$为$d$个$l$比特标量，$P_1,\dots,P_d$为$d$个椭圆曲线$E$上的点，$[a_1]P_1$表示点的标量乘法
$$ [a_1]P_1 + \cdots + [a_d]P_d $$
这一计算在密码学算法中很常见，包括但不限于EdDSA、ECDH，但是朴素的计算方法开销很大，之前有桶方法或者Karabina算法来加速计算
之前的Blog也介绍过MSM，这里跳过
这篇paper提出了对上述两类方法的优化，第一种方法在内存可以保存$2^d-1$个点时非常有效，第二种方法则用于处理超过$2^d-1$个点时的情况
Notation
- $a,b,x\in \Z$：标量
- $C$：矩阵
- $[a,b]$：由整数$a,b$划定的整数区间
- $\mathrm{HW}(x)$：$x$的汉明重量
- $C^T$：矩阵$C$的转置矩阵
- $p$：素数
- $\mathbb K$：由$p$确定的有限域
- $E$：定义在$\mathbb K$上的椭圆曲线
- $\mathcal O$：无穷远点
- $\mathrm A$：$E$上的加法运算
- $\mathrm D$：$E$上的二倍点乘法
Case of d=2
这里先看$d=2$的情况，也就是两个点的标量乘相加$[a_1]P_1+[a_2]P_2$，有很多种方式可以处理，一一介绍
Shamir Trick
Shamir Trick$[2,3]$可以处理此类情况，具体如下
首先将$a_1,a_2$表示为二进制形式，如下
$$ (a_1)_2 = (b^{(1)}_{l-1}b^{(1)}_{l-2}\dots b^{(1)}_1b^{(1)}_0)\\ (a_2)_2 = (b^{(2)}_{l-1}b^{(2)}_{l-2}\dots b^{(2)}_1b^{(2)}_0) $$
然后将这些比特组成一个矩阵$C$
$$ C=\begin{pmatrix}b_{l-1}^{(1)}&\cdots&b_1^{(1)}&b_0^{(1)}\\b_{l-1}^{(2)}&\cdots&b_1^{(2)}&b_0^{(2)}\end{pmatrix} $$
这样一来，我们要计算的式$[a_1]P_1+[a_2]P_2$可以转换为下面这个形式
$$ \begin{aligned} \begin{aligned}[\alpha_1]P_1+[\alpha_2]P_2\end{aligned}& =\quad(P_1,P_2)C\begin{pmatrix}2^{l-1}\\\vdots\\2^1\\2^0\end{pmatrix} \\ &=\quad[2^{l-1}]([b_{l-1}^{(1)}]P_1+[b_{l-1}^{(2)}]P_2)+\cdots+[2^1]([b_1^{(1)}]P_1 \\ & +[b_1^{(2)}]P_2)+[2^0]([b_0^{(1)}]P_1+[b_0^{(2)}]P_2) \end{aligned} $$
其实就是拆分成$l$个点加法$([b^{(1)}_i]P_1+[b^{(2)}_i]P_2)$，然后再将加法结果乘上对应的2的幂$2^i$
分析一下上述方法需要的加法次数，可以看出主要取决于$a_1,a_2$的联合汉明重量（Joint Hamming Weight，JHW），也即矩阵$C$中非零元素的数量
由于$a_1,a_2$可以视为长度为$l$比特的随机标量，因此从统计的角度来看，联合汉明重量$\mathrm{JHW}(a_1,a_2)\approx \frac{3}{4}l$
注意到整个算法包含$l$个步骤，每一步中需要一次二倍点运算（$[2^i]$部分），且至少需要1次群加法运算（如果当前步骤的两个比特$(b_i^{(1)},b^{(2)}_i)$任意一个不为0，则需要一次额外的群加法）
综上，算法整体需要的计算开销为$\frac{3}{4}l\cdot \mathrm A+(l-1)\cdot \mathrm D$（需要预计算$G = P_1+P_2$）
对于标量$a$的表示方式，如果二进制表示中1的数量太多的话（汉明重量太大），可以通过特殊方式降低汉明重量，比如下面这样，其中$\bar 1 = -1$
$$ a=(11100011111011)_2=(100\bar{1}00100000\bar{1}0\bar{1}) $$
如果$\mathrm {HW}(a_1)\approx \mathrm{HW}(a_2)\approx \frac{l}{3}$，则$\mathrm{JHW}(a_1,a_2)\approx\frac{5}{9}$，但是这里引入了负数，因此预计算的点需要多一个：$H = P_1 - P_2$，此时算法的计算开销可以降低为$\frac{5}{9}l\cdot \mathrm A+l\cdot \mathrm D$
Bernsteins's method（DJB）
方法来自$[4]$
考虑蒙哥马利椭圆曲线方程
$$ y^2 = x^3+ax^2+x $$
蒙哥马利曲线有一个特点，可以直接通过相加两个点的横坐标来实现点加法（$P_1 + P_2$）或减法（$P_2 - P_1$）
基于这个特性，Bernstein设计了一个快速计算$[a_1]P_1+[a_2]P_2$的方法，每个滋补粥只需要一次倍点和两次点加法运算，但是这个方法限制了$P_1,P_2$必须为蒙哥马利曲线上的点
这个方法需要引入三个中间变量$T_1,T_2,T_3$，这三个变量在初始化时分别被设置为$T_1 = \mathcal O,T_2 = P_1,T_3 = P_2$，，之后在每个计算步骤都会更新$T_1,T_2,T_3$
下面这个图给出了一个如何计算$13P_1+17P_2$的例子（这里没有用到$(1,-1)$这个中间值）
对于长度为$l$比特的标量，Bernstein的方法需要计算$2l$次点加法和$l$次倍点
Case of d>2
如果有三个或以上的点，此时又该如何完成计算
Hankerson等人提出了一种交织的方法来计算，具体可以看$[5-6]$两篇Refs
另外还有一种桶方法，思想有点类似于桶排序，之前的Blog$[7]$有介绍过
This Paper's methods
本文给出了两种计算方法来计算MSM，假设内存中可以保存至多$2^M - 1$个点，这两种方法分别用于处理$d\leq M$（$d$不太大）和$d< M$（$d$比较大）时的情况
Multidimensional Shamir's trick
第一种方法，多维Shamir trick，用于处理$d$不太大的情况
还是和原版的Shamir trick一样，将$d$个标量拆分为二进制表示，然后组合成矩阵，如下
$$ C=\begin{pmatrix}b_{l-1}^{(1)}&\cdots&b_1^{(1)}&b_0^{(1)}\\b_{l-1}^{(2)}&\cdots&b_1^{(2)}&b_0^{(2)}\\\vdots&\vdots&\ddots&\vdots\\b_{l-1}^{(d)}&\cdots&b_1^{(d)}&b_0^{(d)}\end{pmatrix} $$
然后将待计算的$d$个点拼成一个行向量$P$，并构造一个由2的幂组成的列向量$S$，如下
$$ P=\begin{pmatrix}P_1&P_2&\cdots&P_d\end{pmatrix}\\S=\begin{pmatrix}2^{l-1}\\\vdots\\2^1\\2^0\end{pmatrix} $$
然后MSM的计算过程就可以转化为$P\cdot C\cdot S$，如下
$$ \begin{aligned} \begin{aligned}[\alpha_1]P_1+[\alpha_2]P_2+\cdots+[\alpha_d]P_d\end{aligned}& =P\cdot C\cdot S \\ &=P\begin{pmatrix}b_{l-1}^{(1)}&\cdots&b_1^{(1)}&b_0^{(1)}\\b_{l-1}^{(2)}&\cdots&b_1^{(2)}&b_0^{(2)}\\\vdots&\vdots&\ddots&\vdots\\b_{l-1}^{(d)}&\cdots&b_1^{(d)}&b_0^{(d)}\end{pmatrix}\begin{pmatrix}2^{l-1}\\\vdots\\2^1\\2^0\end{pmatrix} \\ &=[2^{l-1}]([b_{l-1}^{(d)}]P_d+\cdots+[b_{l-1}^{(2)}]P_2+[b_{l-1}^{(1)}]P_1) \\ &+[2^{l-2}]([b_{l-2}^{(d)}]P_d+\cdots+[b_{l-2}^{(2)}]P_2+[b_{l-2}^{(1)}]P_1) \\ &+[2^1]([b_1^{(d)}]P_d+\cdots+[b_1^{(2)}]P_2+[b_1^{(1)}]P_1) \\ &+[2^0]([b_0^{(d)}]P_d+\cdots+[b_0^{(2)}]P_2+[b_0^{(1)}]P_1) \end{aligned} $$
和原版的Shamir trick，需要一些预计算，但是这里不只是预计算一个点这么简单了，需要预计算一个表$T$，表项由所有的$\sum^d_{i = 1,b^{(i)}_j\neq 0} P_i$组成，$T$长下面这样
$$ \begin{aligned} T[0] &= \mathcal O \\ T[1] &= P_1 \\ T[2]&=P_{2} \\ T[3]&=P_{1}+P_{2} \\ T[4]&=P_{3}, \\ T[5]&=P_{1}+P_{3} \\ T[6]&=P_{2}+P_{3} \\ &\vdots \\ T[2^d-1]&=P_{1}+P_{2}+\cdots+P_{d} \end{aligned} $$
这样MSM的计算过程就可以转化为下面这样
$$ [a_1]P_1+[a_2]P_2+\cdots+[a_d]P_d=[2^{l-1}]T[x_{l-1}]+[2^{l-2}]T[x_{l-2}]+\cdots+[2^1]T[x_1]+[2^0]T[x_0]. $$
其中$x_j,T[x_j]$如下（$j\in [0,l-1]$）
$$ \begin{aligned}x_j=\bigoplus_{i=1}^d(b_j^{(i)}<<(i-1))\end{aligned} \\T[x_j] = \sum_{i=1,b_j^{(i)} \neq 0}^{d} P_i $$
这里分析一下复杂度，首先保存这个表$T$需要$2^d-1$的存储空间，然后预计算表$T$的这些表项需要$2^d-d-1$次点加法
对于每一次计算$j\in [0,l-1]$，如果当前步骤中的比特$b^{(i)}_j(1\leq i\leq d)$有至少一个非零，则本步骤的计算需要一次倍点和一次点加法（和上一步骤的运算相加），算法主体（不考虑预计算）需要的倍点运算次数为$l$次，需要的点加法运算与$d$个标量的联合汉明距离相关，即
$$ l\cdot \mathrm{JHW}(a_1,\dots,a_d) \approx l \cdot (1-\frac{1}{2^d}) $$
这里原文给了个简单的例子，用上面的方法计算$13P_1+17P_2+21P_3$
首先将三个标量表示为5比特二进制数（也即$l = 5$），然后排列成矩阵$C$，如下
$$ C= \begin{pmatrix}0&1&1&0&1\\1&0&0&0&1\\1&0&1&0&1\end{pmatrix} $$
然后计算表$T$如下
$$ T = \{\mathcal O,P_1,P_2,P_1+P_2,P_3,P_1+P_3,P_2+P_3,P_1+P_2+P_3\} $$
然后$13P_1+17P_2+21P_3$就可以转化为下面这样
$$ \begin{aligned} [13]P_1+[17]P_2+[21]P_3& =P\cdot C\cdot S \\ &=\begin{pmatrix}P_1,P_2,P_3\end{pmatrix}\begin{pmatrix}0&1&1&0&1\\1&0&0&0&1\\1&0&1&0&1\end{pmatrix}\begin{pmatrix}2^4\\2^3\\2^2\\2^1\\2^0\end{pmatrix} \\ &=\begin{pmatrix}P_1,P_2,P_3\end{pmatrix}\begin{pmatrix}2^3+2^2+2^0\\2^4+2^0\\2^4+2^2+2^0\end{pmatrix} \\ &=[2^4]T[6]+[2^3]T[1]+[2^2]T[5]+[2^1]T[0]+[2^0]T[7] \end{aligned} $$
预计算需要4次点加法，算法主体需要5次倍点和4次点加法，内存开销为7个点
Multidimensional Double and Add
注意到有些算法中的$d$往往很大（比如zkSNARK中的多项式往往可以到$d = 2^{22}$），此时没有任何服务器可以存储这么多点$2^{2^{22}} = 2^{4,194,304}$，不能再按照上述方法实现，因此需要一种新的方法来处理此类情况
首先还是和之前的方法，将每一个标量$a_i$表示为二进制形式，并和之前一样计算$x_j$
这里与上一节的区别在于，如果$x_j\neq 0$，则按照下列方式计算$T[x_j]$
$$ T[x_j]\leftarrow 2\cdot T[x_j] + \sum_{i=1,b_j^{(i)} \neq 0}^{d} P_i $$
这个方法的原理是这样：$x_j = 0$的概率与$\forall i \in [i,d],b^{(i)}_j = 0$的概率相等，当$d$足够大时，此时$x_j = 0$的概率接近于可忽略（具体是$\frac{1}{2^d}$，当$d = 2^{22}$时，$\mathrm{Pr}[x_j = 0] = 2^{-4,194,304}$）
这个方法的每一次计算步骤的计算开销为1次倍点运算（上式的$2\cdot T[x_j]$部分），$\mathrm{HW}(x_j) - 1$次加法（计算上式的求和部分），以及一次加法（将上式的两个计算结果加起来）
综上，总的计算开销为$l$次倍点计算和$\sum_{j = 0}^{l-1} \mathrm{HW}(x_j)$次点加法，如果标量$a_i$满足均匀分布，则$\mathrm{HW}(x_j)\approx \frac{d}{2}$，因此计算开销为
$$ \frac{d}{2}\cdot l\cdot \mathrm A +l\cdot D $$
这里还是给一个例子，计算$17P_1+25P_2+28P_3+12P_4$
这里按照上面的计算思路，有$x_0 = 3,x_1 = 0,x_2 = 12,x_3 = 14,x_4 = 7$，因此计算过程如下
$$ \begin{aligned} [a_{1}]P_{1}+[a_{2}]P_{2}+[a_{3}]P_{3}+[a_{4}]P_{4}& =[17]P_1+[25]P_2+[28]P_3+[12]P_4 \\ &=[16](P_1+P_2+P_3)+[8](P_2+P_3+P_4 \\ &+[4](P_3+P_4)+[2]\mathcal O+(P_1+P_2) \\ &=[2^4]([1]P_1+[1]P_2+[1]P_3+[0]P_4) \\ &+[2^3]([0]P_1+[1]P_2+[1]P_3+[1]P_4) \\ &+[2^2]([0]P_1+[0]P_2+[1]P_3+[1]P_4) \\ &+[2^1]([0]P_1+[0]P_2+[0]P_3+[0]P_4) \\ &+[2^0]([1]P_1+[1]P_2+[0]P_3+[0]P_4). \end{aligned} $$
然后计算开销主要分析汉明重量$\mathrm{HW}(x_i)$，不难计算得到$\mathrm{HW}(x_0) = 2,\mathrm{HW}(x_2) = 2,\mathrm{HW}(x_3) = 3,\mathrm{HW}(x_4) = 3$，因此计算开销为5次倍点运算和10次点加法，需要的存储开销为4个点
安全性分析与效率分析
TODO，可以参考原文$[1]$
References
$[1] On Computing the Multidimensional Scalar Multiplication on Elliptic Curves (iacr.org)
$[2]$ Addition chains of vectors
$[3]$ A public key cryptosystem and a signature scheme based on discrete logarithms
$[4]$ Differential addition chains
$[5]$ Elliptic Curve Cryptography
$[6]$ wnaf*, an efficient left-to-right signed digit recoding algorithm
$[7]$ MSM算法 - 三两老友杂的小岛 (snowolf0620.xyz)