Reed-Solomon码：从入门到放弃 Part 1

开个大坑，分几部分讲吧，也不知道能填多少。

前言#

Reed-Solomon码（下称RS码）是一种应用广泛的纠错码，应用场景包括但不局限于日常使用的QR码（即二维码）、RAID6阵列等等。其本质就是在原始数据中添加一段额外的校验数据，从而实现数据的检查与纠正的目的。

深究其数学理论可能需要群论的知识，但鄙人非数学专业出身，N年前学的数学分析也忘得差不多了，就大概讲讲如何编程实现吧。

参考： https://en.wikiversity.org/wiki/Reed%E2%80%93Solomon_codes_for_coders

目录#

• BCH码
  • BCH码的错误检测
  • BCH码的错误纠正
• Galois Field（有限域）
  • GF加减法运算
  • GF乘法运算
  • GF乘法运算（查表法）
  • GF除法运算
  • GF的幂与逆运算
• 总结

BCH码#

RS码的纠错原理与BCH码类似，先挑简单的下手。

BCH码的错误检测#

首先，我们通过一个算法，对5 bit的数据00011，生成10 bit的校验码11010 11001（后面会讲如何生成），拼接起来形成15bit的数据段：00011 11010 11001（对应十六进制为0xf59）。对其使用某个系数1 01001 10111（0x537）作除法运算（实际上使用的是bit级别的异或代替除法使用的减法运算），如：

                        00011
            _________________
10100110111 ) 000111101011001
               ^ 10100110111
               --------------
                 010100110111
                ^ 10100110111
                  -----------
                  00000000000

从而得到商00011，以及余数00000 00000（余数总共为10 bit，即上式中的后10个bit，而首个bit在执行最后一次异或后恒为0）。

这里先假设数据恒为5 bit，校验码为10 bit，用python实现为：

def bch_checksum(data: int) -> int:
    g = 0x537  # 0b10100110111
    for i in range(4, -1, -1):
        if data & (1 << (i + 10)):
            data ^= g << i
    return data

若没有任何错误（如bit翻转），这5 bit的数据与10 bit校验码拼接，再与某个系数（如0x537）作除法运算后，得到的余数为0，而非0则表示其出现错误。

依旧假设数据为5 bit，校验码为10 bit。编码的工作为将数据左移10 bit，空的位置用0填充，得到填充位置上的余数作为其校验码。这样做保证了在校验时，无错误的数据得到的余数恒为0。编码函数为：

def bch_encode(data: int) -> int:
    return (data << 10) ^ bch_checksum(data << 10)

对00011编码得到的结果与最初的15 bit数据是一致的（即其十六进制为0xf59），不妨自己动手试一下。

BCH码的错误纠正#

通常来说，高效的解码算法是很复杂的，但是在对5 bit的数据面前（也就是只有32种情况下），穷举不妨是一个非常简单粗暴但很直观管用的方法。

用最简单的话来说，我们可以对数据从0到31进行穷举，对每个数据分别进行编码。在32种情况中，再找出与手头上的数据的汉明距离最小（即bit差距最小）时，所对应的数据作为其正确的数据。当然，这种方法在出现多解的情况下，就无法还原正确的数据了。

def hamming_distance(x: int) -> int:
    dist = 0
    while x > 0:
        dist += x & 1
        x >>= 1
    return dist


def bch_decode(data: int) -> int:
    best_guess_data = -1
    best_dist = 15
    for test_data in range(32):
        test_code = bch_encode(test_data)
        test_dist = hamming_distance(data ^ test_code)
        if test_dist == 0:
            return test_data
        if test_dist < best_dist:
            best_dist = test_dist
            best_guess_data = test_data
        elif test_dist == best_dist:
            best_guess_data = -1
    return best_guess_data

从而得到：

>>> bch_decode(0b000111101011001)  # No error
3
>>> bch_decode(0b111111101011001)  # 3 bits error
3
>>> bch_decode(0b111011101011001)  # 4 bits error
-1

最后一种情况说明bit错误太多导致无法纠正其错误。

当然，在实际场景中，对达到GB、TB级别的数据量来说，穷举是行不通的，所以会采用更为复杂的算法对其进行解码，如Berlekamp-Massey。

Galois Field（有限域）#

先假设所有操作数都是8 bit，并且通过加减乘除运算之后得到的结果仍然保持为8 bit。当然，对于结果超过2^8（256）的数都需要对其进行某种取余操作，确保其结果在2^8以内。

Galois Field（下称GF）或者有限域（Finite Field）的定义可以自己去百度，简而言之，一个有限域就是能满足特定运算法则，对其中的每个数字都能进行加减乘除（除以0除外）运算的数字的集合。GF(2^8)则是阶（order）为256的有限域，能表达0~255范围内的整数值。

在数值与多项式的形式转换上，举个例子来说：在GF(2^8)中，多项式数值170（二进制形式10101010）代表的多项式为$1x^7+0x^6+1x^5+0x^4+1x^3+0x^2+1x+0=x^7+x^5+x^3+x$，等价于每个系数都是GF(2)的一元七次多项式。

GF加减法运算#

加减法通过异或运算得出，道理很简单，对2取余的加减法跟异或的运算是一样的。

0101 + 0110 = 0101 - 0110 = 0101 XOR 0110 = 0011

其代表的含义是：
$$(x^2+1)+(x^2+x)=2x^2+x+1=0x^2+x+1=x+1$$

对应的代码：

def gf_add(x: int, y: int) -> int:
    return x ^ y


def gf_sub(x: int, y: int) -> int:
    return x ^ y

GF乘法运算#

与加减法同样，乘法与多项式相乘的过程是类似的，举个例子，10001001乘00101010的过程如下：

$$
\begin{aligned}
(x^7+x^3+1)(x^5+x^3+x)&=x^7(x^5+x^3+x)+x^3(x^5+x^3+x)+1(x^5+x^3+x) \
&=x^{12}+x^{10}+2x^8+x^6+x^5+x^4+x^3+x \
&=x^{12}+x^{10}+x^6+x^5+x^4+x^3+x
\end{aligned}
$$

用异或操作则表示如下：

       10001001
*      00101010
       --------
      10001001
^   10001001
^ 10001001
  --------
  1010001111010

这种运算称为无进位（carryless）的运算，因此在代码中用cl_前缀表示。

对应的代码实现：

def cl_mul(x: int, y: int) -> int:
    z = 0
    i = 0
    while (y >> i) > 0:
        if y & (1 << i):
            z ^= x << i
        i += 1
    return z

当然，上述运算后得到的结果的长度为13 bit，超过了原本的8 bit。因此需要对10011101取余，保留余数，把结果的最高位1降到8 bit。当然，这里的取余操作也是通过异或来完成，从结果的高位到低位以此进行，位数为1时对该数进行异或：

  1010001111010
^ 100011101
  ---------
  0010110101010
  ^ 100011101
    ---------
    00111011110
    ^ 100011101
      ---------
      011000011

因此其代码为：

def cl_div(dividend: int, divisor: int) -> int:
    d1 = dividend.bit_length()
    d2 = divisor.bit_length()  # int占用了多少bit，例如0b101占了3bit
    if d1 < d2:
        return dividend  # 不需要相除
    for i in range(d1 - d2, -1, -1):
        if dividend & (1 << (i + d2 - 1)):
            dividend ^= divisor << i
    return dividend

将两步合二为一，对于GF的乘法，有：

def gf_mul(x: int, y: int, prim: int = 0) -> int:
    # prim不为0时，对结果进行取余，否则返回的是原结果
    z = cl_mul(x, y)
    if prim > 0:
        z = cl_div(z, prim)
    return z

至于选择100011101（0x11d）作为取余的数，数学上的解释挺复杂的，但直观来说就是它代表了8次多项式是不可约的，因此被称为本原多项式（primitive polynomial）或者是不可约多项式（irreducible polynomial）。

对于本原多项式的性质，最重要的一点，就是它的唯一性，它能保证GF中的每个数都有$\alpha$的某一次幂对应，其中$\alpha$被称为是GF中的本原元（primitive element）。举个例子，对于GF(2^8)而言，${\alpha^0,\alpha^1,\ldots,\alpha^{254}}$分别对应1~255中的任意一个值，而$\alpha^{255}=1$，也就是GF中的255次单位根。

由不同的本原元或者不同的本原多项式所构建的GF是有所区别的，这点还请留意，在此处就不再展开了。

GF乘法运算（查表法）#

一般情况下使用$\alpha=2$（对应的多项式为$x$）足矣。在$\alpha=2$下，通过左移运算、与本原多项式进行异或（当数值超过256时）可以快速得到表中每个元素的数值。如：$\alpha^8$本应为1 0000 0000，但大于255的数需要用前面的cl_div进行取余，因此最终的数为00011101。

以下是GF(2^8)中的前几个值：

$$\begin{matrix}
\alpha^0=00000001 & \alpha^4=00010000 & \alpha^8=00011101 & \alpha^{12}=11001101 \
\alpha^1=00000010 & \alpha^5=00100000 & \alpha^9=00111010 & \alpha^{13}=10000111 \
\alpha^2=00000100 & \alpha^6=01000000 & \alpha^{10}=01110100 & \alpha^{14}=00010011 \
\alpha^3=00001000 & \alpha^7=10000000 & \alpha^{11}=11101000 & \alpha^{15}=00100110
\end{matrix}$$

通过这样的方法继续生成后续的数，不难发现$\alpha$的任意次幂之间都不会重复出现，直到$\alpha^{255}=00000001$。

由此，乘法运算可以通过查表（包括正向查表和反向查表）完成：

$$10001001 * 00101010=\alpha^{74} * \alpha^{142}=\alpha^{74+142}=\alpha^{216}=11000011$$

当然，$\alpha$也可以取其他值（比如在GF(2^8)中能取的值有2、3、6、7、9，在$\alpha=3$时对应的本原多项式可以为0x11b，不过这里不再展开讨论，并且后面默认使用$\alpha=2$）。

使用查表进行乘法运算的代码如下：

from typing import Tuple, List


def init_table(prim: int = 0x11d) -> Tuple[List[int], List[int]]:
    table_exp = [0] * 512
    table_log = [0] * 256
    x = 1
    for i in range(255):
        table_exp[i] = x
        table_log[x] = i
        x = gf_mul(x, 2, prim)
    for i in range(255, 512):
        table_exp[i] = table_exp[i - 255]
    return table_exp, table_log

table_exp, table_log = init_table()
def gf_mul_lut(x: int, y: int) -> int:
    if x == 0 or y == 0:
        return 0
    return table_exp[table_log[x] + table_log[y]]

GF除法运算#

根据除法为乘法的逆运算这一性质（即x = gf_div(gf_mul(x, y), y)），可以得到：

def gf_div(x: int, y: int) -> int:
    if y == 0:
        raise ZeroDivisionError()
    if x == 0:
        return 0
    return table_exp[(table_log[x] + 255 - table_log[y]) % 255]

GF的幂与逆运算#

def gf_pow(x: int, power: int) -> int:
    return table_exp[(table_log[x] * power) % 255]


def gf_inv(x: int) -> int:
    return table_exp[255 - table_log[x]]  # = gf_div(1, x)

总结#

代码采用的是wiki里面的，轮子就不重复造了。当然这里的代码也只是教学用，效率肯定不高。有人已经用python写了个pip包（也是按照这篇参考文章写的），直接pip install reedsolo就可以用了，具体用法可以自行研究。我也用过这个包，编码速度是一流的，但是解码速度实在不敢恭维（后面如果写到性能分析的话会详细说说性能瓶颈的）。

如果有机会的话，后面几个part的初步构思如下：

• Part 2：多项式运算、RS编码和解码
• Part 3：拓展GF(2^8)到GF(2^p)，以及对数据进行padding
• Part 4：性能分析以及RS码编码/解码速度与每个参数间的关系
• 番外：CPython的DLL调用吐槽