secp256k1算法详解五（kG点乘多梳状算法）

1 理论基础

在椭圆曲线密码学（ECC）中，kG也称标量乘法运算，即把椭圆曲线上的基点G与标量k进行相乘的运算，结果是椭圆曲线上的另一个点R=kG，其定义为k个G连续相加的结果。该运算是椭圆曲线密钥生成、加解密、签名及验证中的核心运算，所以围绕它产生了多种加速算法，这里仅对secp256k1中前后出现的两种典型算法进行说明。

1.1 窗口查表算法（Windowed lookup）

该算法是secp256k1库早期所使用的算法，其核心思想是将标量k（以256位标量为例）按bits位划分为256/bits个窗口，然后根据预计算查找表查找每个窗口值对应的点，再将每个窗口点相加即可得到最终结果R，其数学公式如下：

以窗口大小bits=4为例，则窗口数为256/4=64个，上述公式可以表示为：

这里ki为4bit二进制数，所以每个窗口包含24=16个查找点（一个窗口对应点0G, 1G, 2G, ..., 15G），则包含16个窗口的查找表为64x16的二维数组，整个表的内存大小为64*16*sizeof(POINT)，假如POINT点x，y坐标都以32字节保存，则查找表所需的内存大小为64*16*(32+32)=64KB。这里如果将bits增大到16，则每个窗口包含216=65536个查找点，窗口个数为256/16=16个，整个查找表大小为16*65536*sizeof(POINT)=64MB。
根据以上分析可知整个kG计算复杂度为O(256/bits)点加操作，例如当bits=4时，需要查找64次表再将对应结果点进行点加即可。
1.2 多梳状算法（Multi-Comb）

为了严格的对应上源码，这里也加入了盲化部分内容，即计算R=kG时不直接进行计算，而是计算R=(k - b)*G + b*G，这样是为了防止攻击者通过分析功耗、电磁辐射等侧信道信息来推测出k，从而引入随机盲化值b，对于盲化后的计算公式，b*G可以预运算获得（对应后续源码中的ge_offset，通过查表计算出(k-b)*G后要加上该值才能得到最终值kG），所以对于任意的k，其实盲化后主要计算部分是(k - b)*G。在计算该部分时，使用了有符号数字多梳状算法，该算法中定义的comb(s, P)函数是一个关键内容：
comb(s, P) = sum( (2*s - 1) * 2^i * P ) for i=0..COMB_BITS-1
其中，s是标量s的第i个比特值（0或1），P是椭圆曲线上点。公式中2*s - 1是一个巧妙的转换：如果s = 1，则2*1 - 1；如果s = 0，则2*0 - 1 = -1，所以上式中对于每个比特位i，不是加上2^i*P就是减去2^i*P，这就是“有符号数字”的含义（系数是+1或-1），可以进一步将该式简化称一个更熟悉的形式：

1. comb(s, P) = sum( (2*s[i] - 1) * 2^i * P )                                （1）
2.            = [ sum( 2*s[i]*2^i ) - sum( 2^i ) ] * P
3.            = [ 2 * sum(s[i]*2^i) - (2^COMB_BITS - 1) ] * P
4.            = [ 2*s - (2^COMB_BITS - 1) ] * P

复制代码

1. 将盲化与梳状算法结合

最直接的想法是在计算(k - b)*G时，让其等于comb(s, G)，即有：
(k - b) * G = comb(s, G) = [ 2*s - (2^COMB_BITS - 1) ] * G
由此可知只要解出相应s，即可通过调用（1）公式计算出(k-b)*G，s可以通过以下公式进行求解：
s = (k - b + (2^COMB_BITS - 1)) / 2 (mod order)     (2)
这个公式需要对2进行模逆运算（相当于除法），在模运算中，除法是复杂且耗时的操作，需要避免。
2. 避免模除2

为了避免昂贵的模除2操作，这里采用了一个巧妙地优化，将公式中基点G除于2，不再计算comb(s, G)，而是计算comb(d, G/2)，仍利用最一开始的公式，则有：

1. comb(d, G/2) = sum( (2*d[i] - 1) * 2^i * (G/2) )                         (3)
2.              = [ 2*d - (2^COMB_BITS - 1) ] * (G/2)
3.              = [ d - (2^COMB_BITS - 1)/2 ] * G

复制代码

现在，令comb(d, G/2) = (k - b) * G，则有：
(k - b) * G = [ d - (2^COMB_BITS - 1)/2 ] * G
得到：
k - b = d - (2^COMB_BITS - 1)/2
最终解出d：
d = k - b + (2^COMB_BITS - 1)/2 (mod order)　　(4)
这里(2^COMB_BITS - 1)/2 是一个常量，可以预先计算。现在，计算d只需要一次模加法和一次模减法，完全避免了模除运算。在后续的源码实现中定义偏移量scalar_offset=(2^COMB_BITS - 1)/2 - b (mod order)，则有d=k+scalar_offset，最终kG=(k-b)*G+b*G=comb(d, G/2)+ge_offset。
3. 梳状含义

在梳状算法中会将标量k拆分成多个位段，每个位段称之为块（Block），每个块中又有T个齿（Teeth），齿与齿之间又有S个间距（Spacing）。以源码中典型取值为例：

1. #define COMB_BLOCKS 11
2. #define COMB_TEETH 6
3. #define COMB_SPACING 4

复制代码

表明会将标量k分成11块，每个块有6个齿，齿间距为4，则每个Block中包含24bits数据，梳子每次可以选中两两间隔为4的6bits数据，在这样的结构下11个块可以包含264bits数据，完全覆盖256bits的标量k（标量k还需进行数据位填充）。

定义mask(b) = sum(2^((b*COMB_TEETH + t)*COMB_SPACING) for t=0..COMB_TEETH-1)，则当b取值从0到COMB_BLOCKS-1时，有以下对应关系：

1. mask(0)  = 2^0   + 2^4   + 2^8   + 2^12  + 2^16  + 2^20,
2. mask(1)  = 2^24  + 2^28  + 2^32  + 2^36  + 2^40  + 2^44,
3. mask(2)  = 2^48  + 2^52  + 2^56  + 2^60  + 2^64  + 2^68,
4. ...
5. mask(10) = 2^240 + 2^244 + 2^248 + 2^252 + 2^256 + 2^260

复制代码

由此可通过这些掩码拆分比特位d，具体来说，每个掩码会被使用COMB_SPACING次，且每次使用时的偏移量不同。

1. d = (d & mask(0)<<0) + (d & mask(1)<<0) + ... + (d & mask(COMB_BLOCKS-1)<<0) +
2.     (d & mask(0)<<1) + (d & mask(1)<<1) + ... + (d & mask(COMB_BLOCKS-1)<<1) +
3.     ...
4.     (d & mask(0)<<(COMB_SPACING-1)) + ...

复制代码

即：

1. table(b, m) = (m - mask(b)/2) * G                             （5）

复制代码

以下是计算过程伪代码：

1. table(b, m) = sum(2^i * (2*d[i] - 1) * G/2)                   （6）

复制代码

在以上伪代码中，梳子是从高位“梳到”低位的，即依次先处理每个块的高位，再依次处理每个块的低位。
4. 折半表

之前已经有结论每个块对应的m可以有2^COMB_TEETH个不同取值，所以对于table(b, m)来说需要有2^COMB_TEETH个表项才能覆盖该b块所有可能取值，但实际上只需一半的表项即可覆盖所有可能取值，即table(b, m)实际只包含2^(COMB_TEETH-1)个表项，以下是分析过程：
由m=(d & mask(b))，令m'是m所有位反转对应的值，则有m'=m XOR mask(b)，进而有table(b, m') = table(b, m XOR mask(b)) = table(b, mask(b) - m)  = (mask(b) - m - mask(b)/2)*G = -(m - mask(b)/2)*G = -table(b, m)。
即如果m'对应是m的所有位反转，那么m'对应的table表项即为m相应表项的负值，对于任意椭圆点P = (x, y)，其负值很容易求得即-P = (x, -y)，所以table(b, m)表只需保存一半的表项，另一半的位反转表项，只需对相应的表项取负即可。对于0块来说，当梳齿对应块中最低位时，有以下对应关系：

例如由m=2^0 + 2^8  + 2^16时的表项（对应点P=(x, y)），可以直接求出m'=2^4 + 2^12 + 2^20的表项（对应点-P = (x, -y)），只需要对m值对应的表项中y值取负即可得到m'值对应的表项。
2 源码详解

2.1 预计算表

函数secp256k1_ecmult_gen_compute_table用于产生第1节中的table(b, m)表项，函数源码如下：
[code]  1 static void secp256k1_ecmult_gen_compute_table(secp256k1_ge_storage* table, const secp256k1_ge* gen, int blocks, int teeth, int spacing) {  2     size_t points = ((size_t)1)