PDF     Slide     Video

这是FAST13 中一篇关于通过Intel SSE3 指令集加速伽罗瓦域上计算的短文。纠错码广泛用于存储系统中，可分为XOR（异或）码和Reed-solomon(RS) 码。前者（LDPC, RAID5, X-code）在编码时仅进行异或运算，速度快；后者编码运算基于伽罗瓦域/有限域GF(2^w)，速度相对较慢（XOR 可看做是GF(2) 上的运算）。该文旨在提出基于处理器128 比特指令集加速编码时的乘法运算。

一、基础知识

首先啰嗦下，我们来了解存储系统如何对数据进行纠错编码，以及和字大小w 的关系，假设有D₀，D₁ 和D₂ 三个磁盘，磁盘内数据按照2^w 的大小划分数据：

D₀ = d₀₀    d₀₁    d₀₂    d₀₃   ……

D₁ = d₁₀    d₁₁    d₁₂    d₁₃   ……

D₂ = d₂₀    d₂₁    d₂₂    d₂₃   ……

当需要对这三个磁盘进行编码生成一个校验盘时，可以通过不同磁盘乘以不同系数相加后得到，这里说的乘法和加法都是在伽罗瓦域上：

C = m₀*D₀ + m₁*D₁ + m₂*D₂

从我们初中知识知道，m₀，m₁ 和m₂ 只要都不为零，那么三块磁盘和校验盘任意失效一个，通过该方程可以计算出失效磁盘的数据，RAID5 的原理就是这样的，当校验磁盘不止一块时：

C₀ = m₀₀*D₀ + m₀₁*D₁ + m₀₂*D₂

C₁ = m₁₀*D₀ + m₁₁*D₁ + m₁₂*D₂

C₂ = m₂₀*D₀ + m₂₁*D₁ + m₂₂*D₂

相乘的系数组成一个编码矩阵M，由线性代数的知识可以知道，编码矩阵M 满秩时，矩阵的各行系数互不相关，失效若干块磁盘可以通过等量的校验盘通过编码矩阵恢复出来。当编码矩阵M 是范德蒙矩阵时，这种编码就是大名鼎鼎的vandermonde-RS码。

下图给出了范德蒙矩阵中的一行系数与磁盘数据相乘生成校验盘C₀ 的示意图：

上文和图中的d_i,j 是一个个字，是编码中计算的最小单位，大小为w 比特，可以看做GF(2^w) 域中的一个元素。图中系数aⁱ  也是GF(2^w) 域中的一个元素。

通过以上分析不难看出，编码中最常遇到的计算以及耗时最多的计算是乘法（相对于加减法都是异或来说），更进一步来说是一个元素乘以一块区域。这篇文章就旨在加速这样的乘法。

二、指令_mm_shuffle_epi8(a, b)

Intel 和其他处理器提供了一些128 比特计算的指令，比如：_mm_set1_epi8(b) 返回拷贝一个比特十六次的128 位向量；_mm_and_si128(a, b) 和_mm_xor_si128(a, b)  返回128 位向量a 和b 的按位与（AND）和按位异或（XOR）的结果；_mm_srli_epi64(a, b)/_mm_slli_epi64(a, b) 将128 位向量看做两个64 比特的值，分别右移/左移b 比特……

_mm_shuffle_epi8(a, b)其中是最重要的一个指令，它将向量a 作为16 个元素（00-ff）的表、向量b 作为16个索引（00-0f）同时进行了16次查询，返回查询结果。这样就大大加速了GF(2⁴)域上的查表速度。

上图是一个例子，b 中第一个比特是0c，将其作为索引查到a 中对应0c（十进制12）中值为5a，那么对应结果v 中第一位结果为5a（注意上图仅用来描述_mm_shuffle_epi8(a, b) 指令的作用）。

GF(2⁴) 中只有16 个元素（0-15），我们为每一个元素乘以其他元素y都用这样的128 位向量来表示，y 作为保存结果的索引，这样可以为域GF(2⁴)  上每个元素构造一个乘法表。比如元素7 的乘法表如下图（table1）：

这个table1 表中的元素表示的是 ‘7 × 索引’ 的值，比如 3 × 7 = 9， 3 × f = b。下面会分别介绍如何利用该指令加速不同w 的查表运算。

三、GF(2⁴) 中 y × A

假设有这样一块数据A，已预先计算了 y = 7 时的乘法表，需要 y×A 的值：

因为索引只有4 bit 大小，那么每次只能计算A 中 8×16 = 64 比特的数据，通过掩码将A 每个字节分高位和低位，分别查表：

比如区域A 中第一个byte 是39，高位是3，低位是9，对应查值为9 和a（3 × 7 = 9 ， 9 × 7 = a），那么结果应该是9a。

区域乘法以低位为例，通过对A 进行掩码得到低位区域 l = mm_and_si128(A, mask1)，低位区域和乘法表进行shuffle 得到结果 l = mm_shuffle_epi8(l, table)。高位通过同样方法，只不过在查表时需要先移到低位查表，再移回高位，才能得到正确结果（下图方法不同，是通过移位低位乘法表，生成一个高位的查找表）。

最后将两次查找的结果异或实现了 y × A 的最终结果：

这个方法写在了开源库GF-complete 的gf_w4.c 文件中的gf_w4_single_table_sse_multiply_region()函数中。

四、GF(2⁸) 中 y × A

从上面GF(2⁴) 中的例子可以看出，虽然_mm_shuffle_epi8(a, b) 指令是128 位的，且每条指令能并行进行16 次查找，但分开来看，每次查找只能运算A 中4 bits。如果y 看做b，A 中一个元素为a，那么GF(2⁸) 中乘法可以写成：

将a 的高四位和低四位乘以值b 的值保存在两个不同的表中，每次将查表的两个结果异或得到最终结果。

相关代码在gf_w8.cc 中gf_w8_split_multiply_region_sse() 函数中：

   1:  581       va = _mm_load_si128 ((__m128i *)(sptr));

   2:  582       t1 = _mm_and_si128 (loset, va);

   3:  583       r = _mm_shuffle_epi8 (mtl, t1);

   4:  584       va = _mm_srli_epi64 (va, 4);

   5:  585       t1 = _mm_and_si128 (loset, va);

   6:  586       r = _mm_xor_si128 (r, _mm_shuffle_epi8 (mth, t1));

   7:  587       va = _mm_load_si128 ((__m128i *)(dptr));

   8:  588       r = _mm_xor_si128 (r, va);

   9:  589       _mm_store_si128 ((__m128i *)(dptr), r);

  10:  590       dptr += 16;

  11:  591       sptr += 16;

loset 是通过 loset = _mm_set1_epi8 (0x0f); 得到的低位掩码，mtl 和mth 分别是低位和高位的乘法表。方法和GF(2⁴) 中类似。

五、GF(2¹⁶) 中 y × A

同理因为索引只有4bit，那么在GF(2¹⁶) 中乘法会分为16/4 = 4部分，每个部分分别是4 bit 索引对16bit 数据，但与GF(2⁸) 和GF(2⁴) 不同的是，每个子乘的结果不是4bits 和8 bits 大小，而是16 bits 大小，于是每个子乘的结果就需要两个表来保存，一张表保存子乘后的高位，一张表保存子乘后的低位。

如果内存数据是在一段区域中连续存放的（文中把这种映射方式称为标准映射方式 stdmap）：

那么对应的是八张乘法表，八次乘法对应着八次查表（分别是a₀ 查表高位/低位结果，a₁ 查表高位/低位结果，a₂ 查表高位/低位结果，a₃ 查表高位/低位结果 ）：

从图中不难看出，每次查表生成4bit 结果，但GF(2¹⁶) 中一次shuffle 运算16 次查表只有8 次查表有效，另外八次查表结果都是0，这是因为乘法表中有一般的位置是以0 填充的，浪费了大量的空间，并增加了查表的次数。究其原因，还是数据映射到内存区域是连续映射的，导致每连续的16 bits 的数据必须分别串行查表。为了提高shuffle 查表的效率，文中提出了变换映射（altmap）的内存映射方法。

将连续的256 bits 内存数据按照高八位和低八位分为高位向量和低位向量两部分：

每一向量部分都分别对应着不同的乘法表，虽然还是8 次查表，但却得到了256 bits 的结果（相应的乘法表中也没有了无效的零元素）。

幸运的是从标准映射（stdmap） 到变换映射（altmap） 同样通过七条SSE 指令完成，而转换回去则只需要两条指令。

六、GF(2³²) 中 y × A

在 w = 32 中的乘法计算和 GF(2¹⁶) 域中计算类似，按4 bits 划分子乘，每个子乘结果是32 bits，于是需要4 张表来保存每个子乘结果：

和 GF(2¹⁶)  中同理，GF(2³²) 中的 y × A 运算也需要altmap 的帮助来加速运算，因为如果内存还是连续的映射的话，对应的乘法表中会有75% 的0（浪费了75% 的空间和效率）。

七、实验结果

实验环境：i7-3770（3.4 GHz，支持最大3.9GHz 睿频，256KB L2 缓存，8MB L3 缓存，支持1333MHz/1600MHz 内存），16 GB 内存（应该和内存频率也有一定关系）。

实验主要测试了传统算法和文章介绍的基于指令集优化的伽罗瓦域计算在区域乘法（y × A）性能。

图中绿色◇（w = 4,8）是查乘法表（Multiplicaiton Tables）的乘法速度，红色× （w = 16）是查对数表（Logarithm Tables）乘法的速度，紫色三角形（w = 32）是裂乘（Split Tables）的乘法速度，黑色＋（w = 4,8,16）是通过Anvin*2 乘2 （By-two）再选择性异或得到结果的乘法速度，这些算法在之前的文章中都有介绍。橙色○（w = 4,8,16）是将区域划分为16 bits 的字，然后将其作为索引查询乘法结果，乘法表预先计算得到，这样对于w = 4 时，相当于每次进行了4次查表，对于w = 8，相当于每次进行了2 次查表。

16-Bit Table 的思想是想使用更多的内存保存更大的乘法表。比如当w = 4 时，乘法表大小为2⁴ * 2⁴ = 256 字节，这对于内存来说也太小气了，于是可以将每两个字和两个字的乘法（8*8）预先得到一个乘法表，不考虑缓存空间的情况下速度提升了两倍，因为每次可以查询两个字的乘法。16-Bit Table 对于每个y 要使用2¹⁶*w/8  = bytes 大小乘法表（w = 4时，32KB；w = 8时，64KB；w = 16时，128KB）。这样当每次y × A 的区域计算时，事先计算出对应y 的乘法表，对于w = 4 时，可以将所有可能的y（16 个）全部计算出来，需要512 KB，而对于w =8或16 时，计算出全部y 对应的乘法表不太现实，所以在区域大于相应的64KB 和128KB 时，16-Bit Table 对于w = 8 和 w =16 才是适用的，这也就是为什么图中超过64KB（w = 8 ）和128KB（w =16）时才有相应的乘法速度（这个乘法速度包括了生成y × A 16-Bit Table 的时间）。