英文原文:Reed–Solomon codes for coders
翻譯參照: Felix021
參考:AN2407.pdf
WIKI:里德-所羅門碼
實現(xiàn)：Pypi ReedSolo

譯注：最近看到了RS碼谤草，發(fā)現(xiàn)還挺有意思的，找了一些資料學習了下宰译，發(fā)現(xiàn)對于程序員來說景埃，從這篇看起會比較容易衷掷。看完以后想著翻譯一下試試，看看自己到底看懂了多少，于是就有了這篇闽寡。本文有部分錯誤，以及一些排版不對的地方尼酿，有興趣的還是看原文更好:)

Reed-Solomon糾錯碼（以下簡稱RS碼）廣泛用于數(shù)據(jù)存儲（如CD）和傳輸應用中爷狈。然而，在這些應用中裳擎，碼字是藏在了電子設(shè)備里涎永，所以無法一窺它們的模樣以及它們是如何生效的。有些復雜的條形碼設(shè)計也采用了RS碼鹿响，能夠暴露出所有的細節(jié)羡微，對于想要獲得這種技術(shù)如何生效的第一手技術(shù)的愛好者，這是一種很有趣的方式惶我。
在這篇文章里妈倔，我是試圖從程序員的視角（而不是數(shù)學家的視角）來介紹RS碼的基本原理。我會用以當下流行的QR碼作為例子來介紹绸贡。我選擇了Python（主要是因為寫出來的代碼看起來整潔美觀）盯蝴，但是我也會介紹一些不那么顯而易見的Python特性，以便那些不熟悉Python的人也能看懂听怕。里頭涉及到的數(shù)學知識對讀者有一定要求捧挺，并且一般是大學才教授的，但是應當能讓對高中代數(shù)掌握較好的人看懂尿瞭。

1. QR碼結(jié)構(gòu)

這一節(jié)詳細介紹QR碼的結(jié)構(gòu)松忍。本節(jié)的信息不完整，這是有意為之筷厘，只介紹了一個小的21x21的QR碼（也被稱為version 1）的常見特征鸣峭。有關(guān)二維碼的更多信息，請參考附錄(appendix)酥艳。
這是一個用于當例子的QR碼摊溶。它由深色和淺色的方格組成，在條形碼領(lǐng)域被稱為“模塊”（module）充石。在角落的3個方形定位器模式是QR碼的典型可見特征莫换。原圖

1.1 掩碼

之所以需要掩碼處理，是為了避免數(shù)據(jù)區(qū)域中出現(xiàn)諸如類似定位器模式的形狀，或者是大片的空白區(qū)域等拉岁，可能會使掃描器混淆坷剧、錯亂。掩碼處理逆轉(zhuǎn)某些模塊（白色變成黑色喊暖，黑色變成白色）惫企，保留其他模塊不變。
參考下面的圖示陵叽，紅色區(qū)域使用一個固定的掩碼模式編碼狞尔，保存了數(shù)據(jù)區(qū)域（黑白部分）的掩碼格式信息。當QR碼被創(chuàng)建的時候巩掺，編碼器通過嘗試偏序，選擇使得不期望出現(xiàn)特征出現(xiàn)最少的那個掩碼模式。被選擇的掩碼模式信息會被保存在格式信息（紅色區(qū)域）中胖替，使得解碼器知道該用哪個研儒。淺灰色的區(qū)域是不包含任何信息的固定模式。此外独令，在定位器模式中殉摔，還包含由交錯的黑白模塊組成的標尺（timing patterns 參考解釋）。原圖

Input   101101101001011
Mask  ^ 101010000010010
Output  000111101011001

1.2 格式信息

格式信息有另一份可辨別的副本挽懦，因此即使其中一份被毀壞翰意，也仍然有機會被識別。副本被分成兩個部分信柿，分別放在另外兩個定位器的邊上冀偶，同樣也是逆時針方向閱讀（沿著左下角定位器往上，然后是右上角定位器邊緣從左往右）渔嚷。
格式信息的前2 bits 給出了用于數(shù)據(jù)的糾錯級別进鸠。這個尺寸的QR碼包含26字節(jié)(bytes,1 byte = 8 bits )信息，其中一些用于保存原數(shù)據(jù)形病，一些用于保存校驗碼客年，如下表所示霞幅。左邊第一列只是給糾錯級別起了個簡單的名字。

格式信息中的接下來3 bits用于指定對數(shù)據(jù)區(qū)域使用的掩碼模式量瓜。模式使用6x6放歌的方式定義司恳，根據(jù)需要不斷重復以覆蓋整個區(qū)域。模式如下所示绍傲，包含了對應的數(shù)學公式指明（掩碼中的）每個模塊是否是黑色（i和j分別是行列從0開始的編號扔傅，從左上角開始算起）原圖

1.3 數(shù)據(jù)

下圖展示了經(jīng)過反掩碼操作后的放大了的圖樣。不同的區(qū)域被標記出來了划提，包括數(shù)據(jù)區(qū)域的邊緣枫弟。原圖

原始信息: 40 d2 75 47 76 17 32 06 27 26 96 c6 c6 96 70 ec
錯誤糾正: bc 2a 90 13 6b af ef fd 4b e0

1.4 解碼

最后的步驟是將信息解碼成可讀格式。前4 bits 指明了信息是如何編碼的哄辣。QR碼使用幾種不同的編碼方案请梢，使得不同類型的信息可以更高效地被存儲，可參見下表的總結(jié)力穗。在方案指示器之后是“長度”字段毅弧，告訴解碼器保存了多少個字符。字符的長度取決于指定的編碼方案当窗。

長度字節(jié)數(shù)只對小的QR碼有效

我們的樣例數(shù)據(jù)開頭是 0100 够坐，表明接下來是每個字符8 bits 。接下來8 bits 是長度字段崖面，00001101（10進制的13）元咙，表明有13個字符。之后才是數(shù)據(jù)字符巫员。前兩個是00100111和01010100（對于ASCII字符 ' 和 T）蛾坯。有興趣的讀者可以自行解碼后續(xù)的字符。（譯注：用微信二維碼掃掃就行了疏遏。脉课。救军。）
在數(shù)據(jù) bit 之后是另外一個4 bit 模式指示器，可以跟前一個不同倘零，從而允許在一個QR碼中混合多個編碼方案唱遭。如果沒有其他數(shù)據(jù)了，用 0000 來標記結(jié)尾（注意呈驶，標準允許忽略這個標記拷泽，如果存儲空間不夠的話）。

2. BCH碼

格式信息是用BCH碼編碼的袖瞻，能夠糾正一定數(shù)量的 bit 錯誤司致。BCH碼是RS碼的普遍形式（所有的RS碼都是BCH碼）。在QR二維碼中聋迎，BCH碼只用于格式信息脂矫，比數(shù)據(jù)信息中用到的RS碼要簡單得多，因此我們可以先從這里開始霉晕。

2.1 BCH錯誤檢測

用于檢測編碼信息的過程類似整數(shù)除法庭再，但是使用異或來代替減法。格式串應該能夠被稱為“生成子”（generator）的碼“整除”牺堰。QR碼使用 10100110111 這個生成子拄轻。下面使用前述格式串 000111101011001 演示了這個過程：

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

下面這個Python函數(shù)實現(xiàn)了這個過程

def qr_check_format(fmt):
    g = 0x537  # = 0b10100110111 in python 2.6+
    for i in range(4,-1,-1):
        if fmt & (1 << (i+10)):    # 判斷是否為1
            fmt ^= g << i
    return fmt

Python注記1：range函數(shù)可能對于非python程序猿來說不夠明確。它產(chǎn)生一個數(shù)字序列 [4, 3, 2, 1, 0]伟葫。源于C的語言中恨搓，它類似于 for (i = 4; i >= 0; i--); 在pascal類語言中則類似于 for i := 4 downto 0 。
Python注記2：&操作符是“按位與”筏养，<<操作符是左移位奶卓，與C類語言一致。

譯注：這個函數(shù)的返回值是參數(shù) fmt 除以生成子的余數(shù)

這個函數(shù)也可以被用于編碼 5-bit 的格式信息：

encoded_format = (format<<10) ^ qr_check_format(format<<10)

譯注：由于format左移了10個bit撼玄，因此這里用 ^ 和用 + 是等價的夺姑。實際上因為這里的 + 是按位加（其實也就是異或了），所以它也等同于 - 掌猛，這一點對于理解它很重要盏浙。如果記格式信息為F，生成子為G荔茬，(F<<10)/G的商為Q废膘、余數(shù)為R (即F<<10 == QG + R)，則最終的編碼信息 C = (F << 10) ^ ((F << 10) mod G) = (QG + R) - R = Q*G慕蔚，從而C應當是能夠被G整除的丐黄。如果收到的C不能被G整除，說明傳輸出錯了孔飒。

讀者也許會對修改此函數(shù)讓它能除以不同數(shù)字感興趣灌闺。例如艰争，更大些的QR碼包含6 bits 的版本信息和12個錯誤校驗碼，并使用生成子 1111100100101 桂对。
使用數(shù)學的規(guī)范格式甩卓，這些二進制數(shù)字可以用一個系數(shù)為"整數(shù)模2"（譯注:參考鏈接第三段內(nèi)容理解）的多項式來描述。數(shù)字中的每一個 bit 是多項式中對應一項的系數(shù)（譯注：該項的冪等于該bit在數(shù)字中的位置）蕉斜。例如：

10100110111 = 1 x^10 + 0 x^9 + 1 x^8 + 0 x^7 + 0 x^6 + 1 x^5
         + 1 x^4 + 0 x^3 + 1 x^2 + 1 x^1 + 1
      = x^10 + x^8 + x^5 + x^4 + x^2 + x^1 + 1

2.2 BCH糾錯

如果qr_check_format(fmt)得到的余數(shù)不是0逾柿，那么這個碼被損壞或者是讀取錯誤了（譯注：即使是0也不能100%保證就對了）。下一步是要找出哪一個格式碼最可能是原數(shù)據(jù)宅此。雖然對于BCH碼有許多復雜的解碼算法机错，但是在這里殺雞用不上牛刀。因為總共只有32個格式串（譯注：15 bits 中前5個是原信息父腕，后10個是根據(jù)原信息生成的）弱匪，因此遍歷找出所有碼字中與fmt不同位數(shù)最小的那個會更簡單（這個被稱為漢明距離, hamming distance）

def hamming_weight(x): #不同bit的數(shù)量
    weight = 0
    while x > 0:
        weight += x & 1
        x >>= 1
    return weight
def qr_decode_format(fmt):
    best_fmt = -1
    best_dist = 15
    for test_fmt in range(0,32):
        test_code = (test_fmt<<10) ^ qr_check_format(test_fmt<<10)
        test_dist = hamming_weight(fmt ^ test_code)
        if test_dist < best_dist:
            best_dist = test_dist
            best_fmt = test_fmt
        elif test_dist == best_dist: #如多個碼字與fmt距離相同，則都不選
            best_fmt = -1
    return best_fmt

如果發(fā)現(xiàn)fmt不能被唯一地解碼侣诵，則函數(shù)qr_decode_format返回-1痢法。這種情況是由于有多個碼字與fmt具有相同的距離狱窘。

>>> print(qr_decode_format(int("000111101011001",2))) # no errors
3
>>> print(qr_decode_format(int("111111101011001",2))) # 3 bit-errors
3
>>> print(qr_decode_format(int("111011101011001",2))) # 4 bit-errors
-1

3. 有限域理論

在講解用于編碼數(shù)據(jù)的RS碼之前杜顺，還需要介紹一點數(shù)學。類似于乘法和除法蘸炸，我們定義兩個對應的操作躬络，應用于 8-bit 字節(jié)（譯注：這里應該指的是整數(shù)，下同）搭儒，且其結(jié)果也是 8-bit 字節(jié)穷当。許多類似的算術(shù)規(guī)則對于新的定義也仍然有效，例如淹禾，任意元素乘以1（幺元）結(jié)果不變馁菜，任意元素乘以0（零元）得0，不允許除以0铃岔。其他有用的數(shù)學屬性（例如分布率）也仍然有效汪疮。（基于這些操作的元素集合）被稱為一個有限域(finite field)，或者叫伽羅華域(Galois field)毁习。更多關(guān)于GF(2^w)的意義

3.1 乘法

乘法是基于上面定義的多項式智嚷。以多項式的形式把輸入寫下來，然后像我們熟悉的那樣纺且，用乘法分配率來計算盏道。我們用 10001001 x 00101010 來舉例：
(x^7 + x^3 + 1) * (x^5 + x^3 + x^1)= x^7 * (x^5 + x^3 + x^) + x^3 * (x^5 + x^3 + x^1) + 1 * (x^5 + x^3 + x^1)= x^12 + x^10 + x^6 + x^5 + x^4 + x^3 + x^1

譯注：由于最后的加法是異或（模2加），因此系數(shù)為偶數(shù)的項都消去了载碌。

相同的結(jié)果也可以通過豎式計算的一個修改版得到猜嘱，只要把其中的加改成異或即可：

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

下面這個Python函數(shù)實現(xiàn)了這個多項式乘法衅枫。

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

當然，由于結(jié)果不能用8-bit來表示了（在這個例子里結(jié)果是13 bits ）泉坐，所以我們在得到最終結(jié)果之前還需要一個步驟：模100011101为鳄，使用前面提到的除法過程。在本例中腕让，最終的結(jié)果是 11000011孤钦。
上述的乘法運算也可以直接用100011101進行異或運算:

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

下面的python代碼用來實現(xiàn)上述等效乘法的異或運算:

def gf_mult_noLUT(x, y, prim=0):
    #GF域乘法，無預算查詢表（precomputed look-up table）輔助纯丸，速度較慢
    #通過使用標準的無進位乘法和不可約多項式規(guī)約實現(xiàn)

    def cl_mult(x,y):
        '''二進制無進位整數(shù)乘法'''
        z = 0
        i = 0
        while (y>>i) > 0:
            if y & (1<<i):
                z ^= x<<i
            i += 1
        return z

    def bit_length(n):
        # 計算整數(shù)最高有效位（整數(shù)二進制格式第一位）. 等效于 int.bit_length()
        bits = 0
        while n >> bits: bits += 1
        return bits

    def cl_div(dividend, divisor=None):
        '''二進制無進位整數(shù)除法偏形，返回余數(shù)'''
        dl1 = bit_length(dividend)
        dl2 = bit_length(divisor)
        # 余數(shù)<除數(shù)
        if dl1 < dl2:
            return dividend
        # 余數(shù)>=除數(shù), 通過移位對齊除數(shù)與余數(shù)的最高有效位
        for i in range(dl1-dl2,-1,-1):
            # 檢查余數(shù)是否可除（為下一次循環(huán)鋪墊）
            if dividend & (1 << i+dl2-1):
                # 可除，則對齊最高有效位并執(zhí)行異或運算(無進位減法)
                dividend ^= divisor << i
        return dividend
    
    ### 調(diào)用GF域乘法 ###
    result = cl_mult(x,y)
    # 用不可約多項式規(guī)約觉鼻，保證結(jié)果仍然在GF域內(nèi)
    if prim > 0:
        result = cl_div(result, prim)
    return result

測試結(jié)果:

>>> a = 0b10001001
>>> b = 0b00101010
>>> print(bin(gf_mult_noLUT(a, b, 0))) # multiplication only
0b1010001111010
>>> print(bin(gf_mult_noLUT(a, b, 0x11d))) # multiplication + modular reduction
0b11000011

至于為什么要用100011101作為約數(shù)進行規(guī)約(reduction)俊扭，其數(shù)學原理較為復雜，簡單來說坠陈，100011101所代表的GF(2^8)多項式是一個“不可約(irreducible)”多項式萨惑，也叫“本原多項式”。100011101(0x11d)是Reed-Solomon (RS)碼中常用的本原多項式仇矾。

譯注：從概念上講庸蔼，本原多項式和不可約多項式并不等價，但在本文適用的范疇內(nèi)是可以等價理解的贮匕，而且完全可以和整數(shù)中的素數(shù)做類比來理解姐仅，只不過這個是在多項式的范疇內(nèi)。

上面的乘法代碼效率較差刻盐，下面是不同算法實現(xiàn)的更高效率的版本：

def gf_mult_noLUT(x, y, prim=0, field_charac_full=256, carryless=True):
    '''采用 Russian Peasant 算法實現(xiàn)GF域整數(shù)乘法 (主要使用位運算, 比上面的方法快).
    當設(shè)定參數(shù)prim = 0 且 carryless=False 時, 返回普通整數(shù)乘法(進位乘法)計算結(jié)果.'''
    r = 0
    while y: # y > 0
        if y & 1: r = r ^ x if carryless else r + x #原文有注釋, 似懂非懂, 譯不了, 同樣的情況后面用*****表示
        y = y >> 1 # 等價于 y // 2
        x = x << 1 # 等價于 x*2
        if prim > 0 and x & field_charac_full: x = x ^ prim #*****
    return r

3.2 基于對數(shù)的乘法

上述過程并不是最適合用來計算伽羅華域乘法的方法掏膏。對兩個數(shù)進行乘法的過程需要8次循環(huán)，除法的過程也需要8次循環(huán)敦锌，然而我們實際上可以通過使用查表的方式來避免循環(huán)馒疹。一個科學的結(jié)果是在內(nèi)存中建立完整的乘法表，但是那需要創(chuàng)建64k大小的表格（譯注：256×256）乙墙。下面給出的結(jié)果要簡潔許多颖变。
首先，注意到使用00000010(=2)來進行乘法是相當簡單的（通常把它記為α或 generator number ）伶丐，只需要向左移1位悼做，然后與100011101進行異或操作（如果有必要）。下面是α的前幾次冪：

α^0 = 00000001 α^4 = 00010000 α^8 = 00011101 α^12 = 11001101
α^1 = 00000010 α^5 = 00100000 α^9 = 00111010 α^13 = 10000111
α^2 = 00000100 α^6 = 01000000 α^10 = 01110100 α^14 = 00010011
α^3 = 00001000 α^7 = 10000000 α^11 = 11101000 α^15 = 00100110

譯注： α^7 * α超過了8-bit哗魂，需要與100011101異或得到α^8肛走，依此類推。

如果這個表格依次類推录别，α的i次冪不會重復朽色，直到α^255=00000001邻吞。由此，這個域中的每個元素葫男，除了0之外抱冷，都是α的某次冪。α被稱為伽羅華域的本原元素（primitive element）或者生成子(generator)梢褐。
仔細觀察旺遮，可以發(fā)現(xiàn)另一種實現(xiàn)乘法的方法：把α的冪加起來。

10001001 * 00101010 = α^74 * α^142 = α^74 + 142 = α^216 = 11000011

問題是盈咳，我們怎么計算10001001是α的幾次冪呢耿眉？這個問題被稱為離散對數(shù)（discrete logarithm），目前沒有已知的高效解法鱼响。然而鸣剪，因為域中只有256個元素，我們可以很容易地建立起一個指數(shù)表丈积；而同時筐骇，對應的對數(shù)表也也有用。

gf_exp = [0] * 512 #512個元素的列表. Python 2.6+可以考慮使用bytearray類型
gf_log = [0] * 256

def init_tables(prim=0x11d):
    #使用參數(shù)prim給的本原多項式計算指數(shù)表和對數(shù)表備用
    global gf_exp, gf_log
    gf_exp = [0] * 512 # anti-log (exponential指數(shù)) table
    gf_log = [0] * 256 # log(對數(shù)) table
    # 計算每一個GF(2^8)域內(nèi)正整數(shù)的指數(shù)和對數(shù)
    x = 1
    for i in range(0, 255):
        gf_exp[i] = x # 存儲指數(shù)表
        gf_log[x] = i # 存儲對數(shù)表
        # 更一般的情況用下面這行江滨，不過速度較慢
        #x = gf_mult_noLUT(x, 2, prim)

        # 只用到 generator==2 或指數(shù)底為 2的情況下铛纬，用下面的代碼速度快過上面的 gf_mult_noLUT():
        x <<= 1 
        if x & 0x100: # 等效于 x >= 256, 但要更快些 (because 0x100 == 256，位運算速度優(yōu)勢)
            x ^= prim # substract the primary polynomial to the current value (instead of 255, so that we get a unique set made of coprime numbers), this is the core of the tables generation

    # Optimization: 雙倍指數(shù)表大小可以省去為了不出界而取模255的運算 (因為主要用這個表來計算GF域乘法牙寞，僅此而已).
    for i in range(255, 512):
        gf_exp[i] = gf_exp[i - 255]
    return [gf_log, gf_exp]

表gf_exp的冗余（實際只需要256饺鹃，但是重復了一遍）是為了簡化乘法莫秆。這樣我們就不用保證gf_log[x]和gf_log[y]會超出表的范圍间雀。

def gf_mul(x,y):
    if x==0 or y==0:
        return 0
    return gf_exp[gf_log[x] + gf_log[y]] # should be gf_exp[(gf_log[x]+gf_log[y])%255] if gf_exp wasn't oversized

3.3 除法

實現(xiàn)對數(shù)表的另一個好處是可以使用冪的差來定義除法。在下面的代碼里镊屎，加上255是為了保證差不是負數(shù)惹挟。

def gf_div(x,y):
    if y==0:
        raise ZeroDivisionError()
    if x==0:
        return 0
    return gf_exp[(gf_log[x] + 255 - gf_log[y]) % 255]

Python注記：raise語句拋出一個異常從而終止gf_div函數(shù)的執(zhí)行。

使用這種方式來實現(xiàn)除法缝驳，對于任意x和非0元素y连锯，gf_div(gf_mul(x, y), y) == x。

3.4 指數(shù)運算和逆運算

對數(shù)表的方法同時還可以簡化并加速指數(shù)運算和逆運算：

def gf_pow(x, power):
    return gf_exp[(gf_log[x] * power) % 255]

def gf_inverse(x):
    return gf_exp[255 - gf_log[x]] # gf_inverse(x) == gf_div(1, x)

3.5 多項式

繼續(xù)介紹RS碼之前用狱，我們需要定義一些伽羅華域多項式（系數(shù)屬于伽羅華域）的操作运怖。這可能會帶來一點混淆，因為伽羅華域的元素本身也是多項式（譯注：是系數(shù)僅為0夏伊、1的多項式）摇展。我的建議是不要想太多。更混亂的是溺忧，x也被用來當作占位符（多項式里的未知數(shù)）咏连。這個x和前面提到的多項式的x沒有任何關(guān)系盯孙，因此別把它們混起來。
前面給出的伽羅華域元素的二進制表示發(fā)在這里用起來會很羅嗦祟滴，因此我換成十六進制振惰。

00000001 x^4 + 00001111 x^3 + 00110110 x^2 + 01111000 x^1 + 01000000= 01 x^4 + 0f x^3 + 36 x^2 + 78 x^1 + 40

在Python里，多項式可以按x的冪遞減的list來表示垄懂，list中的元素是對于項的系數(shù)骑晶，從而上述多項式變成了[0x01, 0x0f, 0x36, 0x78, 0x40]。（實現(xiàn)中也可以用反過來的順序草慧，各有優(yōu)劣透罢。）
第一個函數(shù)將一個多項式和一個標量相乘：

def gf_poly_scale(p,x):
    r = [0] * len(p)
    for i in range(0, len(p)):
        r[i] = gf_mul(p[i], x)
    return r

Python程序員注意：這個函數(shù)不是用Pythonic的方式實現(xiàn)的，它可以用列表推導式(list comprehension)更優(yōu)雅地實現(xiàn)（譯注：指的是可以用return [gf_mul(i, x) for i in p] 來實現(xiàn)）冠蒋。但是我限制使用這些語言特性以便它可以更方便地被轉(zhuǎn)換成其他語言羽圃。

這個函數(shù)將兩個多項式"加"起來（照舊使用異或操作）

def gf_poly_add(p,q):
    r = [0] * max(len(p),len(q))
    for i in range(0,len(p)):
        r[i+len(r)-len(p)] = p[i]
    for i in range(0,len(q)):
        r[i+len(r)-len(q)] ^= q[i]
    return r

下一個函數(shù)將兩個多項式乘起來：

def gf_poly_mul(p,q):
    r = [0] * (len(p)+len(q)-1)
    for j in range(0, len(q)):
        for i in range(0, len(p)):
            r[i+j] ^= gf_mul(p[i], q[j])
    return r

最后，我們需要對一個多項式求值抖剿，給定某個x值朽寞，求出其標量結(jié)果。我們引入秦九韶算法(Horner Scheme, 霍納算法斩郎，中國人很牛逼有沒有!)來避免計算x的n次冪：

01 x^4 + 0f x^3 + 36 x^2 + 78 x^1 + 40= (((01 x^ + 0f) x^ + 36) x^ + 78) x^ + 40

def gf_poly_eval(p,x):
    y = p[0]
    for i in range(1, len(p)):
        y = gf_mul(y,x) ^ p[i]
    return y

還剩一個相對復雜的多項式除法沒有講到脑融，留待下一節(jié)與RS碼一塊講。

4 RS碼

好了缩宜，背景知識都介紹完了肘迎，我們可以開始看看RS碼了。 *譯注吐槽：終于進入重點了锻煌。

4.1 洞察編碼原理

前面費了這么大的篇幅介紹有限域和多項式妓布，是因為它們是糾錯編碼技術(shù)的基礎(chǔ)。通過有限域和多項式的算法宋梧，我們給數(shù)據(jù)添加了一種結(jié)構(gòu)匣沼。盡管信息不變，但是數(shù)據(jù)的結(jié)構(gòu)化允許我們通過定義其上的規(guī)則操作數(shù)據(jù)捂龄，這種獨立于數(shù)據(jù)的結(jié)構(gòu)讓我們可以通過它來恢復損壞的數(shù)據(jù)释涛。

4.2 編碼

4.2.1 簡介

類似于BCH編碼，RS碼通過生成(irreducible generator)和分解(dividing)多項式來表達信息倦沧，分解余下的多項式(remainder)就是RS碼唇撬，最后我們將RS碼附加在原信息后。
BCH編碼乃至其它大多數(shù)的糾錯碼編碼原則是：采用有限的字典(limited dictionary)存儲差異最大的信息元展融，越長的信息元差異往往越大窖认。在上述RS編碼方式中，我們通過將RS碼后綴在原信息碼后面的方式，加長了原碼長耀态；而RS碼基本是獨一無二的轮傍，通過巧妙的算法，可以用RS碼來推斷原信息首装。
通俗的講创夜，類比加密的過程，生成多項式(generator polynomial)就是我們編碼用的字典仙逻，而用它除原信息多項式的運算就是我們加密（RS編碼）的過程驰吓。

4.2.2 錯誤處理(同python，略)

4.2.3 RS生成多項式

RS碼使用類似BCH碼的生成多項式系奉，將 (x - α^n) 乘起來檬贰。例如：

g(x) = (x - α^0) (x - α^1) (x - α^2) (x - α^3) = 01 x^4 + 0f x^3 + 36 x^2 + 78 x^1 + 40

下面這個函數(shù)計算對于給定數(shù)量nsym個校驗符號的RS碼需要的生成多項式：

def rs_generator_poly(nsym):
    g = [1]
    for i in range(0,nsym):
       # g = gf_poly_mul(g, [1, gf_pow(2, i)]) # 指數(shù)運算，跟下面等效
       g = gf_poly_mul(g, [1, gf_exp[i]])
    return g

這個函數(shù)不是很高效缺亮，因為它在每次循環(huán)中陸續(xù)分配了更大的數(shù)組給g翁涤；但是在實際應用中它并不是瓶頸，優(yōu)化癖讀者有興趣的話可以重寫它萌踱，讓它只給g分配一次葵礼。

4.2.4 多項式除法

現(xiàn)成的多項式除法算法中，最簡單的是類似整數(shù)除法的算法并鸵。下面的例子展示了數(shù)據(jù) 12 34 56 （注記：這是16進制表示的）是如何被編碼的：

                             12 da df
                ---------------------
01 0f 36 78 40 ) 12 34 56 00 00 00 00
            ^     12 ee 2b 23 f4
               ----------------
                    da 7d 23 f4 00
                ^   da a2 85 79 84
                  -----------------
                       df a6 8d 84 00
                  ^    df 91 6b fc d9
                     ----------------
                          37 e6 78 d9

余數(shù)和原數(shù)據(jù)連起來得到編碼后的數(shù)據(jù) 12 34 56 37 e6 78 d9鸳粉。這種整數(shù)除法由于算法實現(xiàn)中需要用的很多的規(guī)約循環(huán)導致算法效率低下。綜合除法(synthetic division)是一種更高效的算法园担。下面的代碼是該算法對GF(2^p)多項式的一個擴展實現(xiàn)：

def gf_poly_div(dividend, divisor):
    '''適用于GF(2^p)域的快速多項式除法.'''
    # 注意: 多項式系數(shù)需要按冪次由高到低排序. 例如: 1 + 2x + 5x^2 = [5, 2, 1], 而非 [1, 2, 5]

    msg_out = list(dividend) # 復制被除數(shù)(尾部后綴ecc字節(jié), 用0填充)
    #normalizer = divisor[0] # precomputing for performance
    for i in range(0, len(dividend) - (len(divisor)-1)):
        #msg_out[i] /= normalizer 
        coef = msg_out[i]
        if coef != 0: # 避免log(0)未定義錯誤.
            for j in range(1, len(divisor)): # 因為多項式的首位都是1, (1^1==0)所以可以跳過
                if divisor[j] != 0: # log(0) is undefined
                    msg_out[i + j] ^= gf_mul(divisor[j], coef) # 等價于數(shù)學表達式:msg_out[i + j] += -divisor[j] * coef ,但異或運高效

    # msg_out 包含商和余數(shù), 余數(shù)的最高冪次( == length-1)和除數(shù)一樣, 下面計算分斷點.
    separator = -(len(divisor)-1)
    return msg_out[:separator], msg_out[separator:] # 返回商, 余數(shù).

4.2.5 編碼函數(shù)

下面的這個函數(shù)實現(xiàn)了編碼過程（譯注：可以參考2.1的BCH編碼過程 (fmt<<10) ^ qr_check_format(fmt << 10)届谈，那里有個譯注解釋為什么這么編碼）：

def rs_encode_msg(msg_in, nsym):
    '''Reed-Solomon main encoding function'''
    gen = rs_generator_poly(nsym)

    # 后綴ecc字節(jié)位用0填充, 之后用生成子(irreducible generator polynomial)除
    _, remainder = gf_poly_div(msg_in + [0] * (len(gen)-1), gen)
    # 余數(shù)就是 RS 碼! 后綴到原信息之后形成全部編碼
    msg_out = msg_in + remainder
    # Return the codeword
    return msg_out

簡單吧！事實上弯汰，RS編碼流程中艰山，最后的編碼是最簡單的步驟，它只是一個多項式除法運算蝙泼。而解碼才是RS碼中最難的步驟程剥，因為根據(jù)你需要的不同劝枣，你能找到很多種不同的算法來實現(xiàn)汤踏。（譯注：選擇恐懼癥怕怕！Ｌ蛱凇Ｏ骸）解碼的部分我們以后再講。將之前的多項式綜合除法的代碼和上面的編碼函數(shù)內(nèi)聯(lián)在一起稳诚，得到的下面的代碼就是在很多RS編碼軟件庫中你將看到的樣子：

def rs_encode_msg(msg_in, nsym):
    '''Reed-Solomon main encoding function, using polynomial division (algorithm Extended Synthetic Division)'''
    if (len(msg_in) + nsym) > 255: raise ValueError("Message is too long (%i when max is 255)" % (len(msg_in)+nsym))
    gen = rs_generator_poly(nsym)
    # 用msg_in值初始化 msg_out 并后綴ecc字節(jié), 用0填充.
    msg_out = [0] * (len(msg_in) + len(gen)-1)
    msg_out[:len(msg_in)] = msg_in

    # Synthetic division main loop
    for i in range(len(msg_in)):
        coef = msg_out[i]

        if coef != 0:  # 避免log(0) 未定義錯誤
            for j in range(1, len(gen)): # 因為多項式的首位都是1, (1^1==0)所以可以跳過
                msg_out[i+j] ^= gf_mul(gen[j], coef) # equivalent to msg_out[i+j] += gf_mul(gen[j], coef)

    # 除完之后, msg_out 包含商 msg_out[:len(msg_in)]哗脖，余數(shù) msg_out[len(msg_in):]. 
    #RS碼只用到余數(shù), 所以我們用原信息覆蓋商得到完整編碼.
    msg_out[:len(msg_in)] = msg_in

    return msg_out

下面這個例子演示了編碼函數(shù)如何對第一節(jié)的樣例QR碼中的數(shù)據(jù)進行編碼。它計算出的結(jié)果中的第二行（譯注：輸出的第二行）與樣例QR碼中解碼出來的糾錯碼相符

>>> gf_exp = [0] * 512; gf_log = [0] * 256;
>>> init_tables();
>>> msg_in = [ 0x40, 0xd2, 0x75, 0x47, 0x76, 0x17, 0x32, 0x06,
...      0x27, 0x26, 0x96, 0xc6, 0xc6, 0x96, 0x70, 0xec ]
>>> msg = rs_encode_msg(msg_in, 10)
>>> for i in range(0,len(msg)):
...  print(hex(msg[i]), end=' ')
...
0x40 0xd2 0x75 0x47 0x76 0x17 0x32 0x6 0x27 0x26 0x96 0xc6 0xc6
 0x96 0x70 0xec0xbc 0x2a 0x90 0x13 0x6b 0xaf 0xef 0xfd 0x4b 0xe0

Python版本注記：這是Python 3.0+的例子，因為print的語法變了才避。早版本的Python可以用這句來代替：“print hex(msg[i]), ” （包含末尾的逗號）

4.3 解碼

4.3.1 簡介

RS解碼是從可能損毀的經(jīng)過RS編碼處理的信息中還原信息的過程橱夭。
盡管RS編碼的過程是不變的，但是卻有多種解碼方法桑逝，相應的就有多種解碼算法棘劣。
RS解碼流程大致可以分為5個步驟：

1. 計算伴隨式(syndromes polynomial). 這個方法需要我們通過Berlekamp-Massey等算法分析哪些字符是錯誤的，或者快速判斷輸入的信息是否完全損壞楞遏。
2. 由伴隨式計算錯誤定位式(erasure/error locator polynomial). 用Berlekamp-Massey算法探測哪些字符損毀茬暇。
3. 由前兩式計算錯誤判別式(evaluator polynomial). 判斷有多少字符被篡改的必要步驟。
4. 由前三式計算錯誤等級式(magnitude polynomial). 這個式子保存了有哪些值需要從原信息中減去寡喝，所以也叫損毀式糙俗。換個角度講，在這一步我們提取錯誤的信息保存到這個式子中预鬓。
5. 修復輸入信息. 通過將原輸入信息中的錯誤信息減去來完成修復巧骚。
   下面我們分步驟講解。

4.3.2 伴隨式(Syndrome)計算

（譯注：這個要怎么翻譯啊格二。网缝。Syndrome是癥狀的意思，這里的確也是在計算收到的RS碼碼字的癥狀蟋定，從而判斷是否接受到了錯誤的碼字粉臊。）（譯注：后來看了各種文檔以后，得知一般都翻譯為“伴隨式”）
RS碼的譯碼操作需要多個步驟驶兜。第一個步驟就是計算數(shù)據(jù)的伴隨式扼仲。把數(shù)據(jù)當成一個多項式，使用x = α^0, α^1, α^2, ..., α^n對其求值（譯注：得到n個伴隨式的值）抄淑。因為這些x值使得生成多項式的值為0屠凶，因此如果讀取到的數(shù)據(jù)沒有損壞，結(jié)果應當是0肆资。如果不是矗愧，這些伴隨式里就包含了完成糾錯所必需的信息。計算伴隨式的實現(xiàn)很簡單：

def rs_calc_syndromes(msg, nsym):
    '''給定原信息和糾錯標記數(shù)郑原，計算伴隨式
    從數(shù)學角度看唉韭，這個過程就是一個傅里葉變換 (Chien搜索剛好相反).
    '''    
    synd = [0] * nsym
    for i in range(0, nsym):
        synd[i] = gf_poly_eval(msg, gf_pow(2,i))
    return [0] + synd # pad with one 0 for mathematical precision (else we can end up with weird calculations sometimes)

繼續(xù)上面的例子，我們可以看到（編碼后的）數(shù)據(jù)的伴隨式的確都是0犯犁；而引入一個錯誤以后就會得到非零伴隨式属愤。

>>> synd = rs_calc_syndromes(msg, 10)
>>> print(synd)[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
>>> msg[0] = 0 # deliberately damage the message
>>> synd = rs_calc_syndromes(msg, 10)
>>> print(synd)
[64, 192, 93, 231, 52, 92, 228, 49, 83, 245]

下面的代碼用來自動檢查信息是否有誤：

def rs_check(msg, nsym):
    '''Returns true if the message + ecc has no error, false otherwise '''
    return ( max(rs_calc_syndromes(msg, nsym)) == 0 )

=== 譯注開始 ===
收到的數(shù)據(jù) r(x) 是由編碼后的數(shù)據(jù)（RS碼的碼字） c(x) 和錯誤 e(x) 組合得到的：
r(x) = c(x) + e(x)
記v=nsym由于c(x)在α^0, α^1, α^2, ..., α^(v-1)都是生成多項式的根（去看生成多項式的定義就知道了，注意+實際上也是-）酸役，而c(x)能夠表示為 q(x) * g(x) （不明白的話可以參看2.1節(jié)的譯注）住诸，所以這里對c(x)求值的結(jié)果必然都是0驾胆。
如果對r(x)的求值都是0，說明e(x)也都是0贱呐，說明要么是沒問題丧诺，要么就是錯到根本發(fā)現(xiàn)不出來（正好是另一個碼字）。如果r(x)不等于0奄薇，那么伴隨式S正好就是所有e(x)的值：
S[0] = e(1), S[1] = e(α^1), S[2] = e(α^2), ..., S[v-1] = e(α^(v-1))
如果錯誤的碼字個數(shù)不超過v/2锅必，通過S是可以還原信息的麸粮。
=== 譯注結(jié)束 ===

4.3.3 消除(erasure)糾正

如果錯誤的位置是已知的（譯注：某些信道可以預知秉犹，比如BEC信道，wiki這么說的译秦，我也不知道是啥远搪；不過QR碼是另一個例子）劣纲，糾正它是最簡單的。這被稱為消除糾正谁鳍。對于每一個添加的糾錯碼癞季，都可以糾正一個消除（譯注：這里應該是說，添加的n個糾錯碼能夠保證糾正n個消除碼字）倘潜。如果錯誤位置未知绷柒，那么對于一個錯誤，需要2個錯誤校驗符號涮因。這在實際應用中非常有用废睦，比如QR碼的某些位置被覆蓋或被剪掉什么的。掃描器很難知道發(fā)生了什么养泡，所以不是所有的QR碼掃描器都能糾正消除嗜湃。
有了伴隨式，接下來計算定位式：

def rs_find_errata_locator(e_pos):
    '''誤碼定位
       eg："_"為誤碼澜掩，"h_ll_ worldxxxxxxxxx"中誤碼位置應為： n-1 - [1, 4] = [18, 15] = erasures_loc.'''

    e_loc = [1] # 初始化為1而非0是因為要計算乘法
    # erasures_loc = product(1 - x*alpha**i) for i in erasures_pos and where alpha is the alpha chosen to evaluate polynomials.
    for i in e_pos:
        e_loc = gf_poly_mul( e_loc, gf_poly_add([1], [gf_pow(2, i), 0]) )
    return e_loc

接下來計算錯誤判別式(erasure/error evaluator polynomial)购披，通過一個多項式乘法而后一個多項式除法來實現(xiàn)：

def rs_find_error_evaluator(synd, err_loc, nsym):
    '''Compute the error (or erasures if you supply sigma=erasures locator polynomial, or errata) evaluator polynomial Omega
       from the syndrome and the error/erasures/errata locator Sigma.'''

    # Omega(x) = [ Synd(x) * Error_loc(x) ] mod x^(n-k+1)
    _, remainder = gf_poly_div( gf_poly_mul(synd, err_loc), ([1] + [0]*(nsym+1)) ) # 除法運算只是為了截短

    # Faster way that is equivalent
    #remainder = gf_poly_mul(synd, err_loc) # 乘法
    #remainder = remainder[len(remainder)-(nsym+1):] # 截短

    return remainder

最后，F(xiàn)orney算法實現(xiàn)第4步肩榕，下面的函數(shù)包含了該算法刚陡，實現(xiàn)糾正消除：

def rs_correct_errata(msg_in, synd, err_pos): # err_pos is a list of the positions of the errors/erasures/errata
    '''Forney algorithm, computes the values (error magnitude) to correct the input message.'''
    # 結(jié)合前面的定位結(jié)果進行錯誤判別(errata locator polynomial)
    coef_pos = [len(msg_in) - 1 - p for p in err_pos] # 位置換算為多項式系數(shù) (eg: 由 [0, 1, 2] 換算為 [len(msg)-1, len(msg)-2, len(msg) -3])
    err_loc = rs_find_errata_locator(coef_pos)
    # 計算 errata evaluator polynomial (術(shù)語叫 Omega 或 Gamma)
    err_eval = rs_find_error_evaluator(synd[::-1], err_loc, len(err_loc)-1)[::-1]

    # Chien 算法搜索 err_pos 得到誤碼位置 X  (the roots of the error locator polynomial, ie, where it evaluates to 0)
    X = [] # will store the position of the errors
    for i in range(0, len(coef_pos)):
        l = 255 - coef_pos[i]
        X.append( gf_pow(2, -l) )

    # Forney algorithm: compute the magnitudes
    E = [0] * (len(msg_in)) # 用來保存誤碼數(shù)據(jù). This is sometimes called the error magnitude polynomial.
    Xlength = len(X)
    for i, Xi in enumerate(X):

        Xi_inv = gf_inverse(Xi)

        # 用errata locator的形式導數(shù)(formal derivative)作為Forney算法的分母(denominator), 則第i個error值為： error_evaluator(gf_inverse(Xi)) / error_locator_derivative(gf_inverse(Xi)). 
        err_loc_prime_tmp = []
        for j in range(0, Xlength):
            if j != i:
                err_loc_prime_tmp.append( gf_sub(1, gf_mul(Xi_inv, X[j])) )
        # 下面的乘法計算Forney算法所需的分母 (errata locator derivative)
        err_loc_prime = 1
        for coef in err_loc_prime_tmp:
            err_loc_prime = gf_mul(err_loc_prime, coef)
        # 等價于： err_loc_prime = functools.reduce(gf_mul, err_loc_prime_tmp, 1)

        # Compute y (evaluation of the errata evaluator polynomial)
        # 更直白的Forney算法實現(xiàn)： 
        # Yl = omega(Xl.inverse()) / prod(1 - Xj*Xl.inverse()) for j in len(X)
        y = gf_poly_eval(err_eval[::-1], Xi_inv) # Forney 算法的分子
        y = gf_mul(gf_pow(Xi, 1), y)
        
        # Compute the magnitude
        magnitude = gf_div(y, err_loc_prime) # magnitude value of the error, calculated by the Forney algorithm 
        E[err_pos[i]] = magnitude # store the magnitude for this error into the magnitude polynomial

    # 糾正誤碼! (糾錯碼也會被修正!)
    # (這里并不是Forney 算法, 只是直接套用了解碼結(jié)果)
    msg_in = gf_poly_add(msg_in, E) # 等價于Ci = Ri - Ei ； Ci 為正確信息, Ri 為包含錯誤的信息, Ei 為修正量表(errata magnitudes) .
    return msg_in

Python注記：這個函數(shù)使用[::-1]來反轉(zhuǎn)列表元素的順序株汉。

譯注：這個算法我沒去看筐乳，有興趣的讀者可以看本文提到的另一篇講解AN2407.pdf。

代碼測試如下：

>>> msg[0] = 0
>>> synd = rs_calc_syndromes(msg, 10)
>>> rs_correct_errata(msg, synd, [0]) # [0] is the list of the erasures locations, here it's the first character, at position 0
>>> print(hex(msg[0]))
0x40

4.3.4 錯誤(error)糾正

更可能的情況是未知位置的錯誤郎逃，所以第一步是找出它們哥童。Berlekamp-Massey算法（通常簡稱為BM算法）用來計算錯誤定位多項式。然后我們只需要計算這個（錯誤定位）多項式的零點（譯注：應該就是指多項式的根）褒翰，這標志了錯誤的位置。

def rs_find_error_locator(synd, nsym, erase_loc=None, erase_count=0):
    '''Find error/errata locator and evaluator polynomials with Berlekamp-Massey algorithm'''
    # 使用BM算法迭代計算error locator polynomial.計算差異項 Delta, 由此判斷是否更新 error locator polynomial.

    # Init the polynomials
    if erase_loc: 
        err_loc = list(erase_loc)
        old_loc = list(erase_loc)
    else:
        err_loc = [1] # Sigma(errors/errata locator polynomial) .
        old_loc = [1] # BM 是一個迭代算法，需要對比新舊值(Sigma)來判斷哪些變量需要更新优训。

    synd_shift = 0
    if len(synd) > nsym: synd_shift = len(synd) - nsym

    for i in range(0, nsym-erase_count): 
        if erase_loc: # if an erasures locator polynomial was provided to init the errors locator polynomial, then we must skip the FIRST erase_count iterations (not the last iterations, this is very important!)
            K = erase_count+i+synd_shift
        else: # if erasures locator is not provided, then either there's no erasures to account or we use the Forney syndromes, so we don't need to use erase_count nor erase_loc (the erasures have been trimmed out of the Forney syndromes).
            K = i+synd_shift

        # Compute the discrepancy Delta
        # 計算 Delta: error locator和 the syndromes進行多項式乘法運算朵你。
        #delta = gf_poly_mul(err_loc[::-1], synd)[K] # 嚴格的講應該是 gf_poly_add(synd[::-1], [1])[::-1] , 但對于正確解碼而言并不必要.
        # 而且同時還能優(yōu)化效率: 我們只計算第K項的多項式乘法. 避免套嵌循環(huán).
        delta = synd[K]
        for j in range(1, len(err_loc)):
            delta ^= gf_mul(err_loc[-(j+1)], synd[K - j]) .
        #print("delta", K, delta, list(gf_poly_mul(err_loc[::-1], synd))) # debugline

        # Shift polynomials to compute the next degree
        old_loc = old_loc + [0]

        # 計算 the errata locator and evaluator polynomials
        if delta != 0: # delta非0才更新
            if len(old_loc) > len(err_loc): 
                # 計算 Sigma(errata locator polynomial)
                new_loc = gf_poly_scale(old_loc, delta)
                old_loc = gf_poly_scale(err_loc, gf_inverse(delta)) # err_loc * 1/delta = err_loc // delta
                err_loc = new_loc

            # 更新Delta
            err_loc = gf_poly_add(err_loc, gf_poly_scale(old_loc, delta))

    # 查看結(jié)果是否正確, 或者錯誤太多無法修復
    while len(err_loc) and err_loc[0] == 0: del err_loc[0] # 排除前面的0
    errs = len(err_loc) - 1
    if (errs-erase_count) * 2 + erase_count > nsym:
        raise ReedSolomonError("Too many errors to correct")    # 錯誤太多

    return err_loc

接下來，我們用一種簡單直接的試驗方法揣非，通過這個定位式來定位多項式中的0抡医，進而定位錯誤位置。代碼如下：

def rs_find_errors(err_loc, nmess): # nmess is len(msg_in)
    '''Find the roots (ie, where evaluation = zero) of error polynomial by brute-force trial, this is a sort of Chien's search
    (but less efficient, Chien's search is a way to evaluate the polynomial such that each evaluation only takes constant time).'''
    errs = len(err_loc) - 1
    err_pos = []
    for i in range(nmess): # normally we should try all 2^8 possible values, but here we optimize to just check the interesting symbols
        if gf_poly_eval(err_loc, gf_pow(2, i)) == 0: # It's a 0? Bingo, it's a root of the error locator polynomial,
                                                     # in other terms this is the location of an error
            err_pos.append(nmess - 1 - i)
    # Sanity check: the number of errors/errata positions found should be exactly the same as the length of the errata locator polynomial
    if len(err_pos) != errs:
        # couldn't find error locations
        raise ReedSolomonError("Too many (or few) errors found by Chien Search for the errata locator polynomial!")
    return err_pos

譯注：這個算法我也沒細看早敬，大體上忌傻，它用到了牛頓恒等式（根據(jù)4.3.3計算的伴隨式）來生成錯誤定位多項式，然后求值得到錯誤的位置搞监。更詳細一點的信息可參考AN2407.pdf

這樣的方法效率不高水孩，Chien搜索是一種更高效的算法。<div color="white"（譯注：英文原作者在此留了個坑琐驴，這不是我的風格俘种，以后找到了補上！）</div>
下面是一個糾正了編碼后數(shù)據(jù)中3個錯誤的例子：

>>> print(hex(msg[10]))
0x96
>>> msg[0] = 6
>>> msg[10] = 7
>>> msg[20] = 8
>>> synd = rs_calc_syndromes(msg, 10)
>>> pos = rs_find_errors(synd, len(msg))
>>> print(pos)
[20, 10, 0]
>>> rs_correct_errata(msg, synd, pos)
>>> print(hex(msg[10]))
0x96

4.3.5 消除和錯誤糾正

為了能夠糾正錯誤和消除绝淡，我們需要讓已知的消除不影響查找錯誤位置宙刘。這可以通過計算Forney syndrome來實現(xiàn)，如下所示：

def rs_forney_syndromes(synd, pos, nmess):
    # Compute Forney syndromes, which computes a modified syndromes to compute only errors (erasures are trimmed out). 別混淆Forney syndromes 和 Forney algorithm(根據(jù)錯誤位置糾正信息的算法).
    erase_pos_reversed = [nmess-1-p for p in pos] # prepare the coefficient degree positions (instead of the erasures positions)

    # Optimized method, all operations are inlined
    fsynd = list(synd[1:])      # make a copy and trim the first coefficient which is always 0 by definition
    for i in range(0, len(pos)):
        x = gf_pow(2, erase_pos_reversed[i])
        for j in range(0, len(fsynd) - 1):
            fsynd[j] = gf_mul(fsynd[j], x) ^ fsynd[j + 1]

    # 上述算法是進過優(yōu)化之后的, 等價的理論算法應該是: fsynd = (erase_loc * synd) % x^(n-k)
    #erase_loc = rs_find_errata_locator(erase_pos_reversed, generator=generator) # computing the erasures locator polynomial
    #fsynd = gf_poly_mul(erase_loc[::-1], synd[1:]) # then multiply with the syndrome to get the untrimmed forney syndrome
    #fsynd = fsynd[len(pos):] # 刪除(trim)沒用的前 erase_pos 個系數(shù). 似乎沒必要, 不過能降低后面BM算法的計算量.

    return fsynd

Forney syndrome可以用來替換常規(guī)錯誤位置查找中的syndrome牢酵。
下面的這個rs_correct_errata函數(shù)給出了完整的過程悬包，在msg_in被消除的位置中使用-1來表示。

def rs_correct_msg(msg_in, nsym, erase_pos=None):
    '''Reed-Solomon main decoding function'''
    if len(msg_in) > 255: # can't decode, message is too big
        raise ValueError("Message is too long (%i when max is 255)" % len(msg_in))

    msg_out = list(msg_in)     # copy of message
    # erasures: 為了解碼是debug方便, 設(shè)置為 null (不必要)
    if erase_pos is None:
        erase_pos = []
    else:
        for e_pos in erase_pos:
            msg_out[e_pos] = 0
    # check if there are too many erasures to correct (beyond the Singleton bound)
    if len(erase_pos) > nsym: raise ReedSolomonError("Too many erasures to correct")
    # prepare the syndrome polynomial using only errors.
    synd = rs_calc_syndromes(msg_out, nsym)
    # 查驗 codeword中是否有誤碼. 如果沒有 (伴隨式全部系數(shù)為 0), codeword原樣返回.
    if max(synd) == 0:
        return msg_out[:-nsym], msg_out[-nsym:]  # no errors

    # compute the Forney syndromes, which hide the erasures from the original syndrome (so that BM will just have to deal with errors, not erasures)
    fsynd = rs_forney_syndromes(synd, erase_pos, len(msg_out))
    # compute the error locator polynomial using Berlekamp-Massey
    err_loc = rs_find_error_locator(fsynd, nsym, erase_count=len(erase_pos))
    # locate the message errors using Chien search (or brute-force search)
    err_pos = rs_find_errors(err_loc[::-1] , len(msg_out))
    if err_pos is None:
        raise ReedSolomonError("Could not locate error")    # error location failed

    # Find errors values and apply them to correct the message
    # compute errata evaluator and errata magnitude polynomials, then correct errors and erasures
    msg_out = rs_correct_errata(msg_out, synd, (erase_pos + err_pos)) # note that we here use the original syndrome, not the forney syndrome
                                                                                                                                  # (because we will correct both errors and erasures, so we need the full syndrome)
    # check if the final message is fully repaired
    synd = rs_calc_syndromes(msg_out, nsym)
    if max(synd) > 0:
        raise ReedSolomonError("Could not correct message")     # message could not be repaired
    # return the successfully decoded message
    return msg_out[:-nsym], msg_out[-nsym:] # also return the corrected ecc block so that the user can check()

Python注記：erase_pos和err_pos這兩個數(shù)組用 + 運算連接起來馍乙。

這是實現(xiàn)一個全功能的錯誤/消除糾正RS解碼器所需的最后一個片段玉罐。如果想進一步深究，原作者推薦了一本書潘拨，《Algebraic Codes for Data Transmission》

5 打包一個實例

下面的例子展示了如何使用上文中的代碼進行編碼和解碼：

# Configuration of the parameters and input message
prim = 0x11d
n = 20 # set the size you want, it must be > k, the remaining n-k symbols will be the ECC code (more is better)
k = 11 # k = len(message)
message = "hello world" # input message

# Initializing the log/antilog tables
init_tables(prim)

# Encoding the input message
mesecc = rs_encode_msg([ord(x) for x in message], n-k)
print("Original: %s" % mesecc)

# Tampering 6 characters of the message (over 9 ecc symbols, so we are above the Singleton Bound)
mesecc[0] = 0
mesecc[1] = 2
mesecc[2] = 2
mesecc[3] = 2
mesecc[4] = 2
mesecc[5] = 2
print("Corrupted: %s" % mesecc)

# Decoding/repairing the corrupted message, by providing the locations of a few erasures, we get below the Singleton Bound
# Remember that the Singleton Bound is: 2*e+v <= (n-k)
corrected_message, corrected_ecc = rs_correct_msg(mesecc, n-k, erase_pos=[0, 1, 2])
print("Repaired: %s" % (corrected_message+corrected_ecc))
print(''.join([chr(x) for x in corrected_message]))

應該得到如下的輸出：

Original:  [104, 101, 108, 108, 111,  32, 119, 111, 114, 108, 100, 145, 124, 96, 105, 94, 31, 179, 149, 163]
Corrupted: [  0,   2,   2,   2,   2,   2, 119, 111, 114, 108, 100, 145, 124, 96, 105, 94, 31, 179, 149, 163]
Repaired:  [104, 101, 108, 108, 111,  32, 119, 111, 114, 108, 100, 145, 124, 96, 105, 94, 31, 179, 149, 163]
hello world

6 總結(jié)（略）

譯注："Bobmath"將上述代碼整合以后提交到了pypi吊输，包名叫reedsolo，可以在https://pypi.python.org/pypi/reedsolo找到铁追。稍微要注意的一點是季蚂，上述的所有代碼實際上是針對GF(256)上面的RS碼實現(xiàn)的，且使用了固定的生成多項式琅束，所以并不是完全通用扭屁。