一、前言

伽罗瓦域上的乘法在包括加/解密编码和存储编码中经常使用，常见的AES 和Reed-Solomon 编码就使用了伽罗瓦域GF(2⁸) 中的运算。以2 或者2^w 形式的伽罗瓦域来说，加减法都是异或运算，乘法相对较复杂一些，本文就GF(2^w) 上有限域的乘法运算和优化进行分析。现代计算机是为二进制普通运算所设计，对伽罗瓦域计算最多仅有指令集上的优化，而且仅限于某些处理器，因此在更高层次上优化伽罗瓦域上的乘法显得尤为重要。

二、本原多项式

域中不可约多项式（primitive polynomial）是不能够进行因子分解的多项式，本原多项式是一种特殊的不可约多项式。当一个域上的本原多项式确定了，这个域上的运算也就确定了，本原多项式一般通过查表可得，同一个域往往有多个本原多项式。通过将域中的元素化为多项式的形式，可以将域上的乘法运算转化为普通的多项式乘法模以本原多项式的计算。比如g(x) = x³+x+1 是GF(2³)上的本原多项式，那么GF(2³) 域上的元素3*7 可以转化为多项式乘法：

3*7(in GF(2³)) → (x+1)*(x²+x+1) mod g(x) = x³+1 mod x³+x+1 = x → 3*7 = 2 (in GF(2³))

需要注意的是，系数为2 的整数倍会被约去。

三、乘二运算

无论是普通计算还是伽罗瓦域上运算，乘二计算是一种非常特殊的运算。普通计算在计算机上通过向高位的移位计算即可实现，伽罗瓦域上乘二也不复杂，一次移位和一次异或即可，并且2 是GF(2⁸) 中的生成元。从多项式的角度来看，伽罗瓦域上乘二对应的是一个多项式乘以x，如果这个多项式最高指数没有超过本原多项式最高指数，那么相当于一次普通计算的乘二计算，如果结果最高指数等于本原多项式最高指数，那么需要将除去本原多项式最高项的其他项和结果进行异或，在硬件中是这样实现的（g(x) = x⁸+x⁴+x³+x²+1）：

x₇  ← x₆ ← x₅  ← x₄ ← x₃ ← x₂  ← x₁  ← x₀

↓_____________↑___↑___↑_________↑

这里x₀到x₇ 分别代表域中一个数的比特位，每次乘二除最高位不移位，其余各位向高位移一位，最高位和指定位进行异或（由本原多项式决定）。不难知道最高位异或的那几位对应着本原多项式系数为1 的几项。

用C 语言可以写成：

   1: uint8_t c, cc;

   2: cc = (c<<1)^((c&0x80)?0x1d:0)

c&0x80 可以判断数c 最高位是否为1，如果为1 则异或ox1d。

四、一个例子

本节通过一个例子来看看在有限域内做乘法有哪些方法，以GF(2⁸) （g(x) = x⁸+x⁴+x³+x²+1）上15*15 为例：

4.1 多项式乘法

15*15 写成多项式乘法的形式:

  (x³+x²+1)*(x³+x²+1) mod g(x)

=x⁶+x⁴+x²+1

x⁶+x⁴+x²+1 化成二进制表示形式1010101，即为1+4+16+64 = 85 。

4.2 移位乘法

15 写成二进制表示1111,15*15 = (1111*1111) 为

         1 1 1 1

×       1 1 1 1

—————-

         1 1 1 1

      1 1 1 1 0

   1 1 1 1 0 0

1 1 1 1 0 0 0

——————

                               =   1 0 1 0 1 0 1     → 1+ 4+16+64=85

4.3 乘二乘法

15 写成生成元指数和异或的形式 2³ + 2² + 2 + 1，那么：

15*15 = 15*(2³ + 2² + 2 + 1)

          = 2*(2*(2*15 + 15) + 15) + 15

计算过程如下：

有限域内数                                          操作                                         二进制值                        十进制值

0                                                                                               0 0 0 0 0 0 0 0                        0

15                                                                                             0 0 0 0 1 1 1 1                       15

2*15                                                移位                                   0 0 0 1 1 1 1 0                       30

2*15+15                                          异或                                   0 0 0 1 0 0 0 1                       17

2*(2*15+15)                                    移位                                   0 0 1 0 0 0 1 0                       34

2*(2*15+15)+15                              异或                                   0 0 1 0 1 1 0 1                       45

2*(2*(2*15+15)+15)                        移位                                    0 1 0 1 1 0 1 0                       90

2*(2*(2*15+15)+15)+15                  异或                                    0 1 0 1 0 1 0 1                       85

五、伽罗瓦域上乘法的系统优化

本来这篇文章只想写这部分，但为了能够理解这部分内容，才瞎扯了下前四节的基础内容。下面针对真实系统中伽罗瓦域上乘法的实现和优化进行讨论和分析，现有的乘法有乘法表（Multiplication Table）、对数表（Log/AntiLog Table）、移乘（shift multiplication）和裂乘（split multiplication）四种实现方式。伽罗瓦域的大小为GF(2^w)。

5.1 乘法表

乘法表是一个二维表，保存的是两个乘数伽罗瓦域上相乘的结果，表的大小为2^w · 2^w 。比如m_table[2^w ][2^w ] 即是这样的乘法表，那么当计算a 乘以b 时，查表m_table[a][b] 即可。这样实现方式最简单，只需要一次查表操作（LOOK）即可；缺点就是当w 较大时，内存使用量大，当w = 8和 16 时，分别要使用64KB 和8GB 内存，因此GF(2⁸)上乘法常用乘法表实现。

5.2 对数表

伽罗瓦域GF(2^w) 的生成元g 指的是g⁰ g¹……g^{2^w-2} 都是域上不重复的元素，生成元的周期是域的大小减一（去掉零元素）。那么Log 表建立的是元素a 与以生成元为底的对数log a 的映射，iLog 表建立元素a 与其反对数antiLog a 的映射。那么乘法可以写成对数、反对数和异或的形式：

x·y = iLog(Log x·y) = iLog(Log x +Log y)

x·y = Log(iLog x·y) = Log(iLog x +iLog y)

这样需要保存2^w  大小的对数表和反对数表，每次计算需要三次查表（LOOKUP）和一次异或（AND/XOR）。内存的使用量也大大降低，当w = 8、16 和32  时，分别要使用0.5KB、256KB 和32 GB。不难看出w=32 时，对于一般的机器内存还是难以接受的。

5.3 移乘

移位乘法是这四种方法中唯一一个不需要保存表的方法，而是采用了上节谈到的乘二乘法的方法。以x 乘以y 为例，首先计算出y乘2ⁱ (i = 1~w-1)的结果，然后将x化为二进制形式，分别取出二进制中为1 的位置p 对应的值y*2^p 进行异或得到最终结果。以GF(2⁸) （g(x) = x⁸+x⁴+x³+x²+1）上15*85 为例，首先需要计算出15*2、15*2²、……15*2⁷ 的值，分别表示为a1-a7，a0 = 15，85 二进制表示形式为01010101，那么只需要将有1 位置的0、2、4、6 对应的a0、a2、a4、a6 异或即得到结果，原理如下：

15*85 = 15*(2⁶ + 2⁴ + 2² +2⁰) = 15*2⁶ + 15*2⁴ + 15*2² + 15*2⁰  = a0 + a2 + a4 + a6

这种方法的优势就在于内存实用量少。缺点是运算较为复杂，需要w 次乘二操作，w 次判断和少于w 次的异或操作。 用伪代码表示：

   1: galois_shift_multiply(int x, int y, int w):

   2:     for i = 0 to w-1:

   3:         yy[i] = y*(2^i)

   4:     for j = 0 to w-1:

   5:         if(((1<<j)&x)==1):

   6:             prod = prod^yy[j]

5.4 裂乘

了解上面了三种方法，可以看出前两种方法占用内存较多，步骤少，实现速度快，后一种方法，步骤较多，但节省内存，权衡两种方法的优缺点，裂乘方法诞生了。裂乘的思路是保存乘法表的一部分而不是全部，不能用乘法表直接得到结果的，通过移位到乘法表范围内后再查表得到结果。

以GF(2³²) 为例，x、y 分别是域内元素（uint32_t 类型），每个元素的高八位、低八位和中间两个八位数表示，那么x 和y 可以写成如下形式：

x = x₀ + 2⁸·x₁ + 2^2*8·x₂ + 2^3*8·x₃

y = y₀ + 2⁸·y₁ + 2^2*8·y₂ + 2^3*8·y₃

x₀、x₁、x₂ 和x₃ 分别是x 按八比特位分割得到的四个uint8_t 类型数据。那么x*y 化为：

x*y = ( x₀ + 2⁸·x₁ + 2^2*8·x₂ + 2^3*8·x₃ )*(y₀ + 2⁸·y₁ + 2^2*8·y₂ + 2^3*8·y₃ )

      = x₀y₀

        + 2⁸ ·(x₀y₁ + x₁y₀)

        + 2^2*8 ·(x₀y₂ + x₁y₁ + x₂y₀)

        + 2^3*8 ·(x₀y₃ + x₁y₂ + x₂y₁ + x₃y₀)

        + 2^4*8 ·(x₁y₃ + x₂y₂ + x₃y₁)

        + 2^5*8 ·(x₂y₃ + x₃y₂) 

        + 2^6*8 ·x₃y₃

结果共计7 项，各项系数分别为1、2⁸ 、2^2*8  、2^3*8 、2^4*8 、2^5*8 和2^6*8 ，设α 是uint32_t 类型数中只有低八位有值的任意数字，那么根据x*y 可以归结为α * α 、α * (α<<8) 、 α * (α<<16) 、α * (α<<24) 、(α<<8) * (α<<24) 、(α<<16) * (α<<24) 和 (α<<8) * (α<<24) 七种运算。每种运算用一个乘法表表示，每个乘法表为2⁸*2⁸*4 Bytes 大小。下面是Jerasure 库中生成这七个表的C 语言函数：

   1: int galois_create_split_w8_tables()

   2: {

   3:   int p1, p2, i, j, p1elt, p2elt, index, ishift, jshift, *table;

   4:  

   5:   for (i = 0; i < 7; i++) {

   6:     galois_split_w8[i] = new int [( 1 << 16 )];

   7:   }

   8:     

   9:   for (i = 0; i < 4; i += 3) {

  10:     ishift = i * 8;

  11:     for (j = ((i == 0) ? 0 : 1) ; j < 4; j++) {

  12:       jshift = j * 8;

  13:       table = galois_split_w8[i+j];

  14:       index = 0;

  15:       for (p1 = 0; p1 < 256; p1++) {

  16:         p1elt = (p1 << ishift);

  17:         for (p2 = 0; p2 < 256; p2++) {

  18:           p2elt = (p2 << jshift);

  19:           table[index] = galois_shift_multiply(p1elt, p2elt, 32);

  20:           index++;

  21:         }

  22:       }

  23:     }

  24:   }

  25:   return 0;

  26: }

相应的乘法函数为：

   1: int galois_split_w8_multiply(int x, int y)

   2: {

   3:   int i, j, a, b, accumulator, i8, j8;

   4:  

   5:   accumulator = 0;

   6:  

   7:   i8 = 0;

   8:   for (i = 0; i < 4; i++) {

   9:     a = (((x >> i8) & 255) << 8);

  10:     j8 = 0;

  11:     for (j = 0; j < 4; j++) {

  12:       b = ((y >> j8) & 255);

  13:       accumulator ^= galois_split_w8[i+j][a|b];

  14:       j8 += 8;

  15:     }

  16:     i8 += 8;

  17:   }

  18:   return accumulator;

  19: }

当w=32，表大小分为2⁸*2⁸*4 Bytes ，7张表总计使用内存1.792MB。每次乘法需要6次移位，16次查表。

六、总结

伽罗瓦域上乘法步骤越简单，速度也就越快。第五节中介绍了系统实现伽罗瓦域上的四种算法，一些文献也还介绍了许多其他算法，往往是以上几种算法的变形。通常内存占用小的往往速度越慢，占用内存大的速度快。现实中采用合理的算法达到时间和空间使用率的权衡。