# Error-Correcting Codes (ECC) > 貢獻者: RinHizakura 本文先探討比較容易理解的 Hamming codes,並主要針對 Reed–Solomon codes 及其在 Linux 核心中的應用作探討。 ## 什麼是 ECC? Error-Correcting Codes(ECC) 廣泛用於通訊傳輸,補正因為資料損壞、雜訊所導致的錯誤。 過往當人們談及通訊系統的改進時,幾乎只著眼於如何讓訊號接收得更清楚,像是藉由加大天線或增強訊號的傳輸功率,這些方法都以物理學的角度來改進。 1948 年,[Claude Shannon](https://en.wikipedia.org/wiki/Claude_Shannon) 發表一篇劃時代的論文《通訊的數學理論》([A Mathematical Theory of Communications](https://en.wikipedia.org/wiki/A_Mathematical_Theory_of_Communication)),他挑戰當時的主流思維,他指出,通訊系統的改進不必依賴於實體設備的提升,他認為只要傳送資料的速度低於該頻道的容量上限 (亦即頻道容量 [channel capacity],指通訊頻道能傳輸資訊的速度,以每秒位元數 [bps] 為測量單位),就能大幅降低資料傳輸錯誤的機率,甚至可將錯誤率控制到幾乎為零。 Shannon 的理論不僅顛覆了人們對於通訊系統改進的傳統觀點,更是以一己之力將通訊乃至資料的儲存、壓縮、傳輸,帶進全新的紀元,從而奠定資訊理論 (information theory) 這個學科的基礎。 ```graphviz digraph{ rank = LR data[shape = box, color = white] Encoder[shape = box] Channel[shape = box] Decoder[shape = box] gooddata[label="good data" shape = box color = white] subgraph{ rank = same data->Encoder Encoder->Channel [label = "redundancy" ] Encoder->Channel [label = "data" color=white] Channel->Decoder [label = "noisy redundancy" ] Channel->Decoder [label = "noisy data" color=white] Decoder->gooddata } } ``` 通常伺服器等級主機配置 BMC (一種 out-of-band 系統,裡面往往是小型的 Linux 系統),後者便於監控硬體的錯誤,且有 EDAC 機制,觸發 ECC 修復記憶體錯誤,從而提升伺服器整體的穩定性。於是乎,主流伺服器品牌會改進其 BMC 裡頭的 ECC 機制。 科普短片: [計算機怎麼知道自己錯了?](https://youtu.be/ey0Yk3I2ARY) 短片 [What are Reed-Solomon Codes? How computers recover lost data](https://youtu.be/1pQJkt7-R4Q) 提供 Reed-Solomon (RS) 在內演算法的視覺化解說。 ## Single–Error Correcting (SEC) Codes > 此節內容摘錄並重新整理自 [錯誤更正碼簡介](https://w3.math.sinica.edu.tw/math_media/d184/18404.pdf) 其中的一類被稱 Single–Error Correcting (SEC) Codes,顧名思義,這個方法可以修正編碼中的單位元的錯誤,其中一個經典的例子為 [(7, 4) hamming code](https://en.wikipedia.org/wiki/Hamming(7,4)),意指資料以 7 個位元編碼的資料中,有 4 個位元是真實的資料,其做法大致為: 對於一個資料 `d1 d2 d3 d4`,加上額外的 parity `p1 p2 p3`,此 parity 定義為使下圖中的三個大圓中的 1 數量為偶數,則在僅有至多一碼錯誤的前提下: * 僅有其中一個大圓內 1 的個數不符合,這表示是 `p1` 或 `p2` 或 `p3` 產生錯誤,因此反轉 `p1` 或 `p2` 或 `p3` 即可 * 其中兩個大圓內的 1 的個數不符合,表示是 `d1` 或 `d2` 或 `d3` 產生錯誤,因此反轉 `d1` 或 `d2` 或 `d3` * 若是 `d4` 產生錯誤,三個大圓內的 1 個數皆不符合,反轉 `d4` ![image](https://hackmd.io/_uploads/S16K6Taaa.png) ### Hamming codes 的數學結構 以數學的角度來探討,為什麼 (7, 4) hamming code 可修正一個位元的錯誤呢? 考慮集合$S=\{0,1\}$ 以及其元素 $x_i\in S$ 對 $i=1,2,\dots5$。假設 $d_1,d_2,d_3,d_4$ 對應到 $x_4,x_3,x_2,x_1$,而 $p_3,p_2,p_1$ 對應到 $x_5,x_6,x_7$,則表示其需滿足(注意這裡的加法指的是二進位單一位元的加法,例如: $1+1=0$ ): $$ \left\{ \begin{array}{c} x_1 + x_2 + x_3 + x_5 = 0 \\ x_1 + x_2 + x_4 + x_6 = 0 \\ x_1 + x_3 + x_4 + x_7 = 0 \\ \end{array} \right. $$ 第一個方程對應上述文式圖紅圖部分,第二個方程對應藍圖部分,第三個方程對應綠圖部分。 將上述線性方程式寫成矩陣則是: $$ \begin{pmatrix} 1 & 1 & 1 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 1\\ \end{pmatrix} $$ 令 $C$ 為所有編碼資料的集合即是 H 的 null space,也就是 $\forall x \in C$ $Hx^T = 0$,又 rank($H$) = 3,根據 [Rank–nullity theorem](https://en.wikipedia.org/wiki/Rank%E2%80%93nullity_theorem) 則 rank($C$) = 4 ,表示對於上述式子的解,事實上只要給定 x1~x7 中任意 4 個,其他 3 個就會因應產生,又這任意 4 個可以為任意的 0 或 1,因此共有 $2^4$ = 16 個可能,對應可以編碼的 4 個位元。 接著,我們可以從 Hamming distance 和 Hamming weight 來思考。 * Hamming distance: 指兩個 bitstring 0 和 1 相異的個數,例如 x = 0100, y = 1100,則 Hamming distance $d_H(x,y)$ = 1 * Hamming weight: 一個 bitsting 中非 0 的數量,例如 x = 0100,其 Hamming weight $w_{H}(x)$ = 1 * 我們可以從這兩個定義推廣,發現到 x 和 y 的 Hamming distance 事實上就是 x + y 的 Hamming weight 則我們可以定義一段 codeword 的 minimum distance 為此 codeword 集合中任兩個編碼的最小 Hamming distance,也就是 $$ d_{min}{C} = \min \left\{ d_{H}(x,y) \mid x,y \in C, x \neq y \right\} $$ 而 minimum weight,則定義為 $$ w_{min}{C} = \min \left\{ w_{H}(x)\mid x \in C \right\} $$ 由於 null space $C$ 是同時也是一個 [linear subspace](https://en.wikipedia.org/wiki/Linear_subspace),根據加法封閉性,集合中的任兩個 x、y 相加仍會存在集合 C 中,又前面提到 x、y 相加的 Hamming weight 會等同 x、y 的 Hamming distance,因此可以歸納出 $d_{min}{C} = w_{min}{C}$。 那麼回過頭看,$C$ 的 minimum weight 會是多少呢?,我們知道 $C$ 是 $H$ 的 nullspace,由此下去推理的話: * 1 是不可能的,因為 H 的任一個 column 中都不存在 [0 0 0] * 2 是不可能的,因為 H 的任兩個 column 各異,任取兩個相加後不會得到 [0 0 0] * 如果取第 1、2、7 個 column 相加,可得 [0 0 0],因此 $C$ 的 minimum weight / minumum distance 為 3 由於對於 d 個 bits 的錯誤,其 minumum distance 必須至少要是 2d + 1,因此我們知道這個演算法僅能修正 1 個位元的錯誤 (詳情請參考[錯誤更正碼簡介](https://w3.math.sinica.edu.tw/math_media/d184/18404.pdf)第 7 頁) ### Hamming codes 的矩陣運算 Hamming codes 的編碼與解碼操作可以被轉換成矩陣的運算,首先要先定義編碼的生成矩陣 $G$ 和檢驗的矩陣 $H$ $$ \begin{split} (d4, \space d3, \space d2, \space d1) \\ G = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 1 & 1 & 0 \\ 1 & 1 & 0 & 1 \\ 1 & 0 & 1 & 1 \\ \end{pmatrix}\end{split} $$ 可以看到前四列形成一個單位矩陣,後三列則分別是 p1、p2、p3 所在的大圓中,根據包含的 $d_{n}$ 為何,其對應的位元就設為 1 形成的。 而 $H$ 就是我們前面的線性方程式所形成的矩陣: $$ H = \begin{pmatrix} 1 & 1 & 1 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 1\\ \end{pmatrix} $$ 則對於一段資料,舉例來說 p = $[1 \space 1 \space 0\space 1]^T$,其編碼的計算就是。再次注意到,這裡的矩陣運算是對單一位元的加乘法: $$ x = Gp = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 1 & 1 & 0 \\ 1 & 1 & 0 & 1 \\ 1 & 0 & 1 & 1 \\ \end{pmatrix}\ \begin{pmatrix} 1 \\ 1 \\ 0 \\ 1 \end{pmatrix}\ = \begin{pmatrix} 1 \\ 1 \\ 0 \\ 1 \\ 0 \\ 1 \\ 0\end{pmatrix}\ $$ 假設經過傳輸後得到的向量為 r,如果在傳輸過程中沒有錯誤的位元產生,則會得到 0 向量: $$ z = Hr = \begin{pmatrix} 1 & 1 & 1 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 1\\ \end{pmatrix}\begin{pmatrix} 1 \\ 1 \\ 0 \\ 1 \\ 0 \\ 1 \\ 0\end{pmatrix}\ = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}\ $$ 如果 r 產生 1位元的錯誤,舉例來說,第 6 個位元有錯,則輸入第 6 個位元有錯的 codeword 之後,得到的結果為: $$ z = Hr = \begin{pmatrix} 1 & 1 & 1 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 1\\ \end{pmatrix}\begin{pmatrix} 1 \\ 1 \\ 0 \\ 1 \\ 0 \\ 0 \\ 0\end{pmatrix}\ = \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}\ $$ 因為我們知道對第 6 個位元為 1 的向量運算會得到: $$ z = Hr = \begin{pmatrix} 1 & 1 & 1 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 1\\ \end{pmatrix}\begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 1 \\ 0\end{pmatrix}\ = \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}\ $$ 則表示輸入的 codeword 第 6 個位元出錯,因此需要對其取補數調整,就可得原始的 codeword。(如果我們把 parity位元擺在 2^n 的位置的話,事實上這個矩陣運算可得更美妙(?)的結果,詳情請參考[(7, 4) hamming code](https://en.wikipedia.org/wiki/Hamming(7,4))) ## Reed–Solomon error correction > 此節內容引用譯自 [Reed-Solomon error-correcting code decoder](https://www.nayuki.io/page/reed-solomon-error-correcting-code-decoder#preliminaries) Reed–Solomon error correction 也是 ECC 的一種,透過在原始的資料中加上額外的資訊,則在定義內的錯誤數量以下,在資料傳輸後產生的錯仍可以被修正,這個 ECC 算法被廣泛使用於諸多領域。 ### 先備知識 1. Reed-solomon codes 的操作是對於在 [finite field](https://en.wikipedia.org/wiki/Finite_field) 集合中的 symbol。對於使用者定義的 field $F$ (以電腦領域來說,通常是 $GF(2^8)$),存在 primitive element / generator $\alpha$,這表示 field 中的所有數值可以被表示成非零且不重複的 {$\alpha^0、\alpha^1、...\alpha^{\lvert F\rvert -2}$} :::info 對於 primitive element 的敘述可以參考 [Primitive element (finite field)](https://en.wikipedia.org/wiki/Primitive_element_(finite_field)) 和 [Primitive root modulo n](https://en.wikipedia.org/wiki/Primitive_root_modulo_n)。 ::: 2. 定義資料的長度為 $k$,$k$ 需滿足 $1\leq k < \lvert F\rvert$,意及每個要編碼的資料都是來自 field $F$ 的 $k$ 值數量的 block 3. 定義 $m$ 為延伸的糾正碼長度,$m$ 為滿足 $1 \leq m < \lvert F\rvert - k$ 的整數,這可以糾正上至 $\lfloor m/2 \rfloor$ 的錯誤數量 4. 因此,最終我們會定義 $n = k + m$ 為最後的 block size \ codeword 長度,$n$ 需要滿足 $2 \leq n < \lvert F\rvert$,這表示如果我們想要糾正大量的錯誤,則需要愈大的 field $F$ ### Systematic encoder 1. Reed-Solomon error-correcting codes 具有不同的形式與處理方式,這裡所說明的是比較普遍的 [BCH code](https://en.wikipedia.org/wiki/BCH_code)。這意味著原始資料會視為是一個多項式的係數,而接續的糾正碼也會接續著成為新的多項式項。 2. 根據 $\alpha$ 和 $m$ 定義一個生成多項式 $$ g(x) = \prod_{i=0}^{m-1} (x-\alpha^i) $$ 則此多項式為 m degree,且最大冪次項 $x^m$ 的係數為 1 3. 假設原始的資料是由 $k$ 個值組成的 ($M_0$, $M_1$, ... ,$M_k-1$),$M_i$ 是 field $F$ 中的一個數值,則將以這 $k$ 個值形成的多項式定義為 $$ M(x) = \sum\limits_{i= 0}^{k-1}{M_ix^i} $$ 4. 則 Reed-solomon codeword 的編碼被定義為 $$ s(x) = M(x)x^m - [(M(x)x^m) \space mod \space g(x)] $$ 這裡我們可以看到, $M(x)x^m$ 會形成 $x^m$, $x^{m+1}$ 到 $x^{m+k-1}$,總共 $k$ 項,而這 $k$ 項的係數其實也是資料的內容。 而因為 $[(M(x)x^m) \space mod \space g(x)]$ 會是 degree $m - 1$ 更低,因此不會有和 $M(x)x^m$ 的計算互動,因此整個 $s(x)$ 的 degree 會為 $m + k - 1$ 也就是 $n - 1$ 4. 因此,最後的編碼即為 $$ s(x) = \sum\limits_{i= 0}^{n-1}{s_ix^i} $$ 多項式的係數,也就是 $(s_0,s_1,...,s_{n-1})$ ### Peterson-Gorenstein-Zierler decoder #### Calculating syndromes 1. 假收我們收到的編碼為 $(r_0,r_1,...r_{n-1})$,這此編碼形成的多項式為: $$ r(x) = \sum\limits_{i= 0}^{n-1}{r_ix^i} $$ 2. 對於接收者而言,我們並不知道 $s_0, s_1, ... , s_{n-1}$ 是多少,但我們可以定義誤差 $e_0 = r_0 -s_0\text{, } e_1 = r_1 - s_1,...,e_{n-1} = r_{n-1}-s_{n-1}$,則誤差多項式為 $$ e(x) = r(x) - s(x ) = \sum\limits_{i= 0}^{n-1}{e_ix^i} $$ 3. 定義 $m$ 個 syndrome value $S_i$ 中,$0 \leq i < m$ ,則 $$ S_i = r(\alpha^i) = s(\alpha^i) + e(\alpha^i) = 0 + e(\alpha^i) = e(\alpha^i) $$ * 這是因為 $s(x)$ 可以被 $g(x)=(x-\alpha^0)(x-\alpha^1)...(x-\alpha^{m-1})$ 整除(參考 Systematic encoder 的第四點),代表 syndrome 的值只會根據 error 而改變 * 所以,如果所有的 syndrome value 皆為 $0$ 的話,這表示得到的編碼是正確的,不需要做任何的修正。 4. 如果把所有的 syndrome 方程式都列出來,那就是 $$ \left\{ \begin{array}{r} S_0 =& e_0\alpha^{0 \times 0} + e_1\alpha^{1 \times 0} + ... + e_{n-1}\alpha^{(n-1) \times 0}\\ S_1 =& e_0\alpha^{0 \times 1} + e_1\alpha^{1 \times 1} + ... + e_{n-1}\alpha^{(n-1) \times 1}\\ ... \\ S_{m-1} =& e_0\alpha^{0 \times (m-1)} + e_1\alpha^{1 \times (m-1)} + ... + e_{n-1}\alpha^{(n-1) \times (m-1)} \end{array} \right. $$ 5. 並重寫成矩陣形式: $$ \begin{bmatrix} e_0\alpha^{0 \times 0} + e_1\alpha^{1 \times 0} + ... + e_{n-1}\alpha^{(n-1) \times 0} \\ e_0\alpha^{0 \times 1} + e_1\alpha^{1 \times 1} + ... + e_{n-1}\alpha^{(n-1) \times 1} \\ ...\\ e_0\alpha^{0 \times (m-1)} + e_1\alpha^{1 \times (m-1)} + ... + e_{n-1}\alpha^{(n-1) \times (m-1)} \end{bmatrix} = \begin{bmatrix} S_0 \\ S_1 \\ ... \\ S_{m-1} \end{bmatrix} $$ 再把誤差項單獨取出 $$ \begin{bmatrix} \alpha^{0 \times 0} & \alpha^{1 \times 0} & ... & \alpha^{(n-1) \times 0} \\ \alpha^{0 \times 1} & \alpha^{1 \times 1} & ... & \alpha^{(n-1) \times 1} \\ ... \\ \alpha^{0 \times (m-1)} & \alpha^{1 \times (m-1)} & ... & \alpha^{(n-1) \times (m-1)} \end{bmatrix} \begin{bmatrix} e_0 \\ e_1 \\ ... \\ e_{n-1} \end{bmatrix} = \begin{bmatrix} S_0 \\ S_1 \\ ... \\ S_{m-1} \end{bmatrix} $$ 由於 $\alpha$ 是已知的,且 $S_i$ 也為已知,理論上我們可以對 $e_i$ 求解,不過可以看到對於 $n - 1$ 個未知數,僅有 $m - 1$ 個方程式,顯然該如何得到誤差項(或者比較準確的說,預測)就需要一些技巧了 #### Finding error locations 1. 我們需要先找到錯誤的位置。假設訊息中有 $\nu$ 個錯誤,如我們在前提中所說,$\nu$ 需滿足 $1 \leq \nu \leq \lfloor\frac{m}{2}\rfloor$。所以說,除非有時間或空間的限制,我們會期望 $\nu$ 可以大一些,可以最大化能糾正的錯誤 2. 假設這 $\nu$ 錯誤的位置在 $I_0,I_1,...I_{\nu-1}$,$I_i$ 之間是沒有順序的且滿足 $0 \leq I_i < n$,這表示在這位置上的誤差 $e_{I_0},e_{I_1},...,e_{I_{\nu-1}}$ 會是任意值(可能是 $0$),而其他的 $e_i$ 則必須等於 $0$ 3. 接著,再根據錯誤的地方 $I_i$ 定義,對於 $0 \leq i < \nu$ $$ X_i = \alpha^{I_i} \\ Y_i = e_{I_i} $$ 4. 因為我們知道其他的 $e_i$ 都為 $0$,因此可以重新定義 syndrome,把前面的矩陣方程式改寫成 $$ \begin{bmatrix} X^{0}_0 & X^{0}_1 & ... & X^{0}_{\nu-1} \\ X^{1}_0 & X^{1}_1 & ... & X^{1}_{\nu-1}\\ ... \\ X^{m-1}_0 & X^{m-1}_1 & ... & X^{m-1}_{\nu-1} \end{bmatrix} \begin{bmatrix} Y_0 \\ Y_1 \\ ... \\ Y_{\nu-1}\end{bmatrix} = \begin{bmatrix} S_0 \\ S_1 \\ ... \\ S_{m-1}\end{bmatrix} $$ 5. 接著,定義一個 error locator polynomial $$ \Lambda(x) = \prod_{i=0}^{\nu-1} (1-X_ix) \\ = 1+\Lambda_{1} x+\Lambda_{2} x^2+...\Lambda_{\nu} x^\nu $$ 6. 因為 $(1-X_iX_i^{-1}) = 0$,因此我們知道 $$ 0 = \Lambda(X_{i}^{-1}) = 1+\Lambda_{1} X_i^{-1}+\Lambda_{2} X_i^{-2}+...\Lambda_{\nu} X_i^{-\nu} $$ 7. 現在,對於 $0 \leq i < \nu$ 且 j $\in \mathbb{Z}$,將兩側同乘上 $Y_iX_i^{j+\nu}$ $$ 0 = Y_iX_i^{j+\nu}\Lambda(X_{i}^{-1}) = Y_iX_i^{j+\nu} + \Lambda_{1}Y_iX_i^{j+\nu} + ... \Lambda_{\nu}Y_iX_i^{j+\nu} $$ 8. 然後把範圍中的所有 $i$ 得到的結果相加 $$ \begin{aligned} 0 &= \sum\limits_{i= 0}^{\nu-1}{Y_iX_i^{j+\nu}\Lambda(X_{i}^{-1})} \\ &= \sum\limits_{i= 0}^{\nu-1}{(Y_iX_i^{j+\nu} + \Lambda_1 Y_iX_i^{j+\nu - 1} + \Lambda_2 Y_iX_i^{j+\nu - 2} + ... +\Lambda_\nu Y_iX_i^{j})}\\ &= \sum\limits_{i= 0}^{\nu-1}{Y_iX_i^{j+\nu}} + \Lambda_1 \sum\limits_{i= 0}^{\nu-1}{Y_iX_i^{j+\nu - 1}} + ... +\Lambda_\nu \sum\limits_{i= 0}^{\nu-1}{Y_iX_i^{j}} \\ &= S_{j+\nu} + \Lambda_1 S_{j+\nu - 1} + ... + \Lambda_\nu S_{j} \end{aligned} $$ 然後再將式子調整成 $$ -S_{j+\nu} = \Lambda_\nu S_{j} +\Lambda_{\nu-1} S_{j+1}+ ... +\Lambda_1 S_{j+\nu - 1} $$ 其中 $0 \leq j < \nu$ 9. 則如果將每個 $j$ 都列出等式: $$ \left\{ \begin{array}{r} \Lambda_\nu S_0 + \Lambda_{\nu-1} S_1 + ... + \Lambda_{1} S_{\nu-1} = -S_\nu \\ \Lambda_\nu S_1 + \Lambda_{\nu-1} S_2 + ... + \Lambda_{1} S_{\nu} = -S_{\nu+1}\\ ...\\ \Lambda_\nu S_{\nu-1} + \Lambda_{\nu} S_{\nu} + ... + \Lambda_{1} S_{2\nu-2} = -S_{2\nu-1} \end{array} \right. $$ 一樣可以將此轉換成矩陣運算的形式: $$ \begin{bmatrix} S_0 & S_1 & ... & S_{\nu-1} \\ S_1 & S_2 & ... & S_{\nu} \\ ...\\ S_{\nu-1} & S_\nu & ... & S_{2\nu-2} \\ \end{bmatrix} \begin{bmatrix} \Lambda_\nu \\ \Lambda_{\nu-1} \\ ... \\ \Lambda_{1}\end{bmatrix} = \begin{bmatrix} -S_\nu \\ -S_{\nu+1} \\ ... \\ -S_{2\nu-1}\end{bmatrix} $$ 10. 現在,我們終於得到了可以求解的方程式了,我們可以將其寫成擴增矩陣: $$ \begin{bmatrix} S_0 & S_1 & ... & S_{\nu-1} & -S_\nu\\ S_1 & S_2 & ... & S_{\nu} & -S_{\nu+1}\\ ...\\ S_{\nu-1} & S_\nu & ... & S_{2\nu-2} & -S_{2\nu-1}\\ \end{bmatrix} $$ 並且透過 [Gauss-Jordan elimination](https://online.stat.psu.edu/statprogram/reviews/matrix-algebra/gauss-jordan-elimination) 求解。 如果對此線性矩陣的解為 [inconsistent](https://en.wikipedia.org/wiki/Consistent_and_inconsistent_equations)(無解),那表示 codeword 中有超過 $\nu$ 數量的錯誤,因此無法得到解。如果我們的 $\nu$ 可以更大,則應該將其調大並且重新解碼,否則,就是無法完成 error correction。 如果是 consistent 但是 under-determined(無限多解),這可能發生於錯誤少於 $\nu$ 個錯誤時,因為任何一個 $e_i = 0$ 的位置都可以被視為是一個 virtual error,此時就必須把 $\nu$ 降低並重新解碼。 11. 一旦我們得到係數 $\Lambda_1$、$\Lambda_2$、...$\Lambda_\nu$,則此時可以把 x = $\alpha^{-i}$($0 \leq i < n$) 代入 $\Lambda(x)$ 中進行運算看看 $\Lambda(\alpha^{-i})$ 是否為 0,表示此位置有誤差,即 $e^i$ 不為 0,把這個 i 加入 $I_i$ 如果找到不為 $0$ 的 $e^i$ 多餘 $\nu$ 個,表示傳輸後的資料錯誤太多,無法修正。如果少於 $\nu$ 個,表示誤差比預計中的少,我們只需將 $\nu$ 重新定義為較低的數字,並刪除編號較高的 $I_i$ 變量即可。 對於大於等於 $n$ 的 $i$ 毋需被代入,這是因為錯誤不能出現在長度為 $n$ 的資料以外。 #### Calculating error values 1. 現在,我們已經得知了錯誤的位置 $I_0,I_1,...I_{\nu-1}$,因此我們就知道 $X_i$ = $\alpha^{I_i}$, $0 \leq i < \nu$ 2. 由於 $Y_i = e_{I_i}$ ,且我們已知 $S_i$,所以 $$ \begin{bmatrix} X^{0}_0 & X^{0}_1 & ... & X^{0}_{\nu-1} \\ X^{1}_0 & X^{1}_1 & ... & X^{1}_{\nu-1}\\ ... \\ X^{m-1}_0 & X^{m-1}_1 & ... & X^{m-1}_{\nu-1} \end{bmatrix} \begin{bmatrix} Y_0 \\ Y_1 \\ ... \\ Y_{\nu-1}\end{bmatrix} = \begin{bmatrix} S_0 \\ S_1 \\ ... \\ S_{m-1}\end{bmatrix} $$ 現在就可以求解了。 同樣可以透過 [Gauss-Jordan elimination](https://online.stat.psu.edu/statprogram/reviews/matrix-algebra/gauss-jordan-elimination) 來求解。 如果 linear system 是 consistent 且 independent 的,則存在唯一解,可以成功修正錯誤。如果是 inconsistent,則表示錯誤無法被修正,或者結果是 consistent 但 dependent,那表示沒有唯一的解(即使唯一解是必須的)。 3. 則對於修正後的訊息 $r'(x) = r'_0x^0+r'_1x^1+...r'_{n-1}x^{n-1}$,對於在錯誤位置 $I_i$ $0 \leq i < \nu$ 的部份,根據 $r'_{I_i} = r_{I_i} \!-\! Y_i = r_{I_i} \!-\! e_{I_i}$ 去修正,其他則維持 $r_i' = r_i$ 即可。 ### 註記 decode 的部份,Peterson-Gorenstein-Zierler decoder 是$Θ(m^3+mk)$ 的解,但在速度的考量上,則應該透過 [Berlekamp-Massey algorithm](https://en.wikipedia.org/wiki/Reed%E2%80%93Solomon_error_correction#Berlekamp%E2%80%93Massey%20decoder) 求出位置,加上 [Forney algorithm](https://en.wikipedia.org/wiki/Forney_algorithm) 還原錯誤,其中兩者的時間複雜度皆為 $Θ(m^2)$ ## 如何實作 Reed–Solomon codes? > 此節內容譯自: [Reed–Solomon codes for coders](https://en.wikiversity.org/wiki/Reed%E2%80%93Solomon_codes_for_coders) 這裡引用另一篇文章,來理解相關的操作可以如何被實作。 ### 基本想法 想像一下,資料在傳輸的過程中有資料遺失了,最終我們收到的資料叫作 `a p p * e`。不過我們可以猜到遺失的字母是什麼,因為當我們把字典翻開,我們找到字典中最接近的字就是 `apple`。 那麼如果我們收到的資料是 `b * d` 呢? 這時候它應該要是 `bad` 還是 `bed`? 為了解決這個問題,我們可以在原始資料後面加上一下額外的字母然後傳輸,例如: ``` bad -> badxyz bed -> bedqwe ``` 這時候,假如我們收到的資料有遺失,變成 `b * d x * z` ,那麼對照字典,最接近的字應該是 `badxyz`,所以我們就有足夠的證據猜測原始的資料應該是 `bad` 了。 ### BCH code [BCH code](https://en.wikipedia.org/wiki/BCH_code) 是 Reed–Solomon codes 經常使用的編碼形式,在前面的章節也有提到,可以簡單理解為是透過多項式的係數來表達資料。 #### BCH error detection BCH 的錯誤檢查類似整數長除法,不過減法的運算會改用 XOR。舉例來說,對於 000111101011001 這個 formating infomation 以及 10100110111 這個 generator,如果可以整除(餘數為 0),那就是表示沒有錯誤發生,此時得到的商就是原始的資料,否則的話就表示有錯誤需要被修正。 :::info 因為原文中是以 QR 為例來說明的,關於 [formating infomation](https://en.wikiversity.org/wiki/Reed%E2%80%93Solomon_codes_for_coders#Formatting_information) 請參考原文,還蠻有趣的XD,有空我就把完整內容整理上來。我們可以先簡單理解成 000111101011001 的前 5 個位元是真正的資料,後面的 10 個位元則是用來修正錯誤的編碼。 ::: ``` 00011 ---------------- 10100110111 ) 000111101011001 ^ 10100110111 -------------- 010100110111 ^ 10100110111 ------------- 00000000000 ``` 對應的 python code 實作如下 ```python def qr_check_format(fmt): g = 0x537 # = 0b10100110111 in python 2.6+ for i in range(4,-1,-1): if fmt & (1 << (i+10)): fmt ^= g << i return fmt ``` 所以我們也可以反過來用此來編碼 5 位元的 formating infomation ```python format = 0x3 encoded_format = (format<<10) ^ qr_check_format(format<<10) print(hex(encoded_format)) # 反過來可得 formating infomation 000111101011001 ``` #### BCH error correction 如果得到的餘數不為 0,表示資料有誤,那麼我們就要「猜出」原本的資料為何。對於 BCH 的解碼有許多演算法,不過這裡我們可以用一個簡單的方式來理解: 因為我們的原始資料只有 32 種可能,因此我們可以一一比較,並以 hamming distance 最小的那個作為猜出的結果 ```python def hamming_weight(x): 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: best_fmt = -1 return best_fmt ``` 事實上 Reed–Solomon 就是依據類似的概念來解碼的,不過顯然資料很大的時候,我們就不能採取這種搜尋方法了,一些比較複雜的作法,如 Berlekamp-Massey,就必須被使用。 ### Finite field arithmetic 要了解 Reed–Solomon,就不得不對 [finite field](https://en.wikipedia.org/wiki/Finite_field) 有一定的了解。 我們需要定義對於一個 8 位元數字加法、減法、乘法、除法,且經過這些運算後仍為 8 bit,一個比較 naive 的作法是,我們可以定義這些運算和自然數的作法相同,並 mod 256 以避免 overflow,對於這些範圍的數字,我們就將其叫作 Galois Field 2^8,這個作法事實上是可以成立的(你可能會質疑在除法的可行性?下面會解釋之)。 要簡單的說明什麼是 Galois Fields 的話: 一個 finite field 就是在有限集合中的一群數字,在這個 field 裡,我們需要定義加法、減法、乘法、除法滿足: [封閉性](https://en.wikipedia.org/wiki/Closure_(mathematics))、[結合律](https://en.wikipedia.org/wiki/Associative_property)、[交換律](https://en.wikipedia.org/wiki/Commutative_property)、[分配律](https://en.wikipedia.org/wiki/Distributive_property)、[Identity](https://simple.wikipedia.org/wiki/Identity_Property) 和 inverse 六大 property。 例如,實數域 ℝ 是一個 field,但整數域 ℤ 就不是,因為不滿足乘法的 inverse property(找不到整數域中的 x 可以滿足 x * 4 = 5)。為了處理此問題,一個可行的作法是用質數來當作 modulo,例如 257,或者是任何質數的正整數數次方。也就是說,整數域 ℤ modulo 任何的質數都可以被稱為 Galois Field。 而 modulo 2 則是有趣的一類,例如這裡的例子中,因為 256 = $2^8$,這是質數 2 的 8 次方,這使得 GF(2^8) 的計算定義有一些特性。 :::warning 細節可以參考[原文](https://en.wikiversity.org/wiki/Reed%E2%80%93Solomon_codes_for_coders#Introduction_to_mathematical_fields)以及 [Non-prime fields](https://en.wikipedia.org/wiki/Finite_field#Non-prime%20fields) ::: #### 加減法 對於 GF(2^n) ,加法和減法可以同時被 XOR 運算所替代,這很符合邏輯,因為相加 modulo 2 就像 XOR,而相減 modulo 2 和相加 modulo 2 並沒有不同。 ``` 0101 + 0110 = 0101 - 0110 = 0101 XOR 0110 = 0011 ``` 或者從多項式的角度去思考 $(x^2 + 1) + (x^2 + x) = 2 x^2 + x + 1 = 0 x^2 + x + 1 = x + 1$ 注意到這是因為我們是在 GF(2^n) 中,才可以用 XOR 取代加減法,在其他的 GF 中則不一定如此。 ```python= def gf_add(x, y): return x ^ y def gf_sub(x, y): return x ^ y ``` #### 乘法 乘法的運算是基於多項式的乘法,例如 10001001 乘以 00101010: $$ \begin{align} &(x^7 + x^3 + 1) * (x^5 + x^3 + x^1)\\ &= x^7 * (x^5 + x^3 + x^1) + 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 \end{align} $$ 這等同於 ``` 10001001 * 00101010 ---------- 10001001 ^ 10001001 ^ 10001001 ------------- 1010001111010 ``` 注意要把加法替換成 XOR。 顯然,此時得到的位元已經超過原本的 8 bits,因此我們需要 modulo 一個神祕的數字 100011101,同等的操作也可以寫成: ``` 1010001111010 ^ 100011101 ---------------- 0010110101010 ^ 100011101 ---------------- 00111011110 ^ 100011101 --------------- 011000011 ``` 對應的程式碼為: ```python def gf_mult_noLUT(x, y, prim=0): '''Multiplication in Galois Fields without using a precomputed look-up table (and thus it's slower) by using the standard carry-less multiplication + modular reduction using an irreducible prime polynomial''' ### Define bitwise carry-less operations as inner functions ### def cl_mult(x,y): '''Bitwise carry-less multiplication on integers''' z = 0 i = 0 while (y>>i) > 0: if y & (1<<i): z ^= x<<i i += 1 return z def bit_length(n): '''Compute the position of the most significant位元(1) of an integer. Equivalent to int.bit_length()''' bits = 0 while n >> bits: bits += 1 return bits def cl_div(dividend, divisor=None): '''Bitwise carry-less long division on integers and returns the remainder''' # Compute the position of the most significant位元for each integers dl1 = bit_length(dividend) dl2 = bit_length(divisor) # If the dividend is smaller than the divisor, just exit if dl1 < dl2: return dividend # Else, align the most significant 1 of the divisor to the most significant 1 of the dividend (by shifting the divisor) for i in range(dl1-dl2,-1,-1): # Check that the dividend is divisible (useless for the first iteration but important for the next ones) if dividend & (1 << i+dl2-1): # If divisible, then shift the divisor to align the most significant bits and XOR (carry-less subtraction) dividend ^= divisor << i return dividend ### Main GF multiplication routine ### # Multiply the gf numbers result = cl_mult(x,y) # Then do a modular reduction (ie, remainder from the division) with an irreducible primitive polynomial so that it stays inside GF bounds if prim > 0: result = cl_div(result, prim) return result ``` 為什麼需要 mod 100011101 呢,其數學證明可能會有些複雜,不過簡而言之,這是為了確保 $GF(2^8)$ 多項式是 "irreducible" (可以想成是一個質數的多項式,所以不能被任兩個 degree 小的多項式相乘而得到),這個數字被稱做 primitive polynomial \ irreducible polynomial \ prime polynomial, 100011101 (0x11d) 是其中常用的一種。至於該如何生成這個數字呢,請參考 [Universal Reed-Solomon Codec](https://en.wikiversity.org/wiki/Reed%E2%80%93Solomon_codes_for_coders/Additional_information#Universal_Reed-Solomon_Codec)。 #### Multiplication with logarithms 雖然前面的作法可以達成 GF 中數字的乘法,然而並不是一個高效率的作法,一個更好的作法是透過查表,不過8位元數字會有 256 種結果,我們也沒必要把 256 x 256 的乘法結果完全的紀錄起來。下面介紹一種更精簡的作法: 首先,定義 $\alpha$ 為一個 generator number = 2,則對於 $\alpha^i$ 的計算就只是對 $\alpha$ 做位元的位移,然後和 100011101 做 XOR (如果有 overflow) 運算而已。我們可得如下: $$ \begin{split} & \alpha^0 = 00000001, \alpha^1 = 00000010, \alpha^2 = 00000100, \alpha^3 = 00001000 \\ & \alpha^4 = 00010000, \alpha^5 = 00100000, \alpha^6 = 01000000, \alpha^7 = 10000000 \\ & \alpha^8 = 00011101, \alpha^9 = 00111010, \alpha^{10} = 01110100, \alpha^{11} = 11101000 \\ & \alpha^{12} = 11001101, \alpha^{13} = 10000111, \alpha^{14} = 00010011, \alpha^{15} = 00100110\\ \end{split} $$ 如果依此類推建立表格,從 $\alpha^0$ 直到 $\alpha^{255}$ 前的數值都不會有重複(這是 GF 的特性),也就是說,除了 0 之外的數都會是某個 $\alpha^i$。 那麼,我們知道 10001001 乘以 00101010 就會是某兩個 $\alpha$ 相乘 $$ 10001001 * 00101010 = \alpha^{74} \times \alpha^{142} = \alpha^{216} = 11000011 $$ 因此我們只要建立起對於 $\alpha^i$ 查到的值,以及反過來對於一個值對應的 $\alpha^i$ 就可以有效的進行乘法運算了! ```python gf_exp = [0] * 512 # Create list of 512 elements. In Python 2.6+, consider using bytearray gf_log = [0] * 256 def init_tables(prim=0x11d): '''Precompute the logarithm and anti-log tables for faster computation later, using the provided primitive polynomial.''' # prim is the primitive (binary) polynomial. Since it's a polynomial in the binary sense, # it's only in fact a single galois field value between 0 and 255, and not a list of gf values. global gf_exp, gf_log gf_exp = [0] * 512 # anti-log (exponential) table gf_log = [0] * 256 # log table # For each possible value in the galois field 2^8, we will pre-compute the logarithm and anti-logarithm (exponential) of this value x = 1 for i in range(0, 255): gf_exp[i] = x # compute anti-log for this value and store it in a table gf_log[x] = i # compute log at the same time x = gf_mult_noLUT(x, 2, prim) # If you use only generator==2 or a power of 2, you can use the following which is faster than gf_mult_noLUT(): #x <<= 1 # multiply by 2 (change 1 by another number y to multiply by a power of 2^y) #if x & 0x100: # similar to x >= 256, but a lot faster (because 0x100 == 256) #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: double the size of the anti-log table so that we don't need to mod 255 to # stay inside the bounds (because we will mainly use this table for the multiplication of two GF numbers, no more). for i in range(255, 512): gf_exp[i] = gf_exp[i - 255] return [gf_log, gf_exp] ``` 稍微瀏覽一下程式碼,我們可以看出其中的 `gf_exp` 是用來計算從 $\alpha^1$ 到 $\alpha^{256}$ 的數值的,而 `gf_log` 則是反過來的 `index`。也可以看到程式碼中紀錄的事實上是 $\alpha^1$ 到 $\alpha^{512}$,這是使我們在查表時 gf_log[x] + gf_log[y] 超過範圍也可以找到對應的 i。 ```python def gf_mul(x,y): if x==0 or y==0: return 0 return gf_exp[gf_log[x] + gf_log[y]] ``` #### 除法 透過對 $\alpha$ 建立的表格,則除法的運算也就可以一併加速了 ```python 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] ``` 注意這裡 `+ 255` 可以保證 index > 0 #### Power and Inverse 同樣的,冪次和 inverse 也可以一併加速 ```python 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) ``` #### Polynomials 接下來,我們還需要再定義一些對於多項式係數的操作。因為 GF 的元素本身,做為多項式的係數也是個多項式,所以勢必不要搞混了。以下面的多項式為例: $$ 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 $$ 在程式中,我們可以把這個係數在矩陣中寫成 [0x01, 0x0f, 0x36, 0x78, 0x40],根據實作的不同考量,也可以把順序反過來。 則第一個 function 用來與一個 scalar 相乘 ```python 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 ``` 將兩個多項式相加 (記得用 XOR) ```python 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 ``` 將兩個多項式相乘 ```python= def gf_poly_mul(p,q): '''Multiply two polynomials, inside Galois Field''' # Pre-allocate the result array r = [0] * (len(p)+len(q)-2) # Compute the polynomial multiplication (just like the outer product of two vectors, # we multiply each coefficients of p with all coefficients of q) for j in range(0, len(q)): for i in range(0, len(p)): r[i+j] ^= gf_mul(p[i], q[j]) # equivalent to: r[i + j] = gf_add(r[i+j], gf_mul(p[i], q[j])) # -- you can see it's your usual polynomial multiplication return r ``` 最後,我們需要一個可以對特定的 x 求出多項式值的 function。這裡使用了 [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method) 使得 $x^i$ 的計算可以迭代累積,不需要每次都重複計算 ```python def gf_poly_eval(poly, x): '''Evaluates a polynomial in GF(2^p) given the value for x. This is based on Horner's scheme for maximum efficiency.''' y = poly[0] for i in range(1, len(poly)): y = gf_mul(y, x) ^ poly[i] return y ``` 最後,我們還需要一個較為複雜的多項式除法,這個留待下一節來講解。 --- ## Linux 核心中的 Reed–Solomon codes 終於來到 Reed–Solomon codes 的關鍵實作了。 在前一節的文章中,使用 python 實作 Reed–Solomon codes,作為認識該如何實作 RS 是個很好的出發點。接著我們從另一個應用來嘗試了解 Reed-solomon codes 的實作,那就是 linux 中的 [linux/lib/reed_solomon/](https://github.com/torvalds/linux/tree/master/lib/reed_solomon)。沒錯! linux 中是有相關的 API 可以使用的。 下面,我們就以 [Reed-Solomon Library Programming Interface](https://www.kernel.org/doc/html/latest/core-api/librs.html#reed-solomon-library-programming-interface) 為出發點,嘗試從 API 的使用深入到程式內部吧! ### `rs_control` 我們可以在 [linux/include/linux/rslib.h](https://github.com/torvalds/linux/blob/0fe5f9ca223573167c4c4156903d751d2c8e160e/include/linux/rslib.h) 找到相關的資料結構: ```c struct rs_control { struct rs_codec *codec; uint16_t buffers[]; }; ``` * `codec`: 整個 Reed-Solomon 的編碼相關參數 * `buffer`: 在 `decode_rs()` 內部會使用的 buffer ```c struct rs_codec { int mm; int nn; uint16_t *alpha_to; uint16_t *index_of; uint16_t *genpoly; int nroots; int fcr; int prim; int iprim; int gfpoly; int (*gffunc)(int); int users; struct list_head list; }; ``` rs_codec 則是整個 Reed-Solomon 中最關鍵的資料結構: * `mm`: 每個 symbol 總共的位元數量 * `nn`: 每個 block 中共包含幾個 symbol (就是 `(1<<mm)-1`) * `alpha_to`: 為了增加在 [Galois field](https://en.wikipedia.org/wiki/Finite_field) 中乘法運算的效率,會建立每個 $\alpha^i$ 的表 * `index_of`: 同樣考量效率,對於每個值,也要可以反回去找到對應的 $alpha^i$ * `genpoly`: 生成式的係數矩陣 * `nroots`: 生成多項式的 degree,也就是額外的 parity 數量 * `fcr`: first consecutive root(index form)。這可以使計算出生成多項式時的 consecutive root 改變,通常都設為 0 * 舉例來說,對於 degree 3 的生成多項式,如果指定 fcr = 3,計算出的生成多項式就是 $g(x) = \prod_{i=3}^{5} (x-\alpha^i)$ * 也因此有範圍限制( 舉例來說 255 for GF(2^8))) * `prim`: 表示用 Galois field 中的 primitive element / generator $\alpha^{prim}$ 來當 generator,參照 [Primitive element (finite field)](https://en.wikipedia.org/wiki/Primitive_element_(finite_field)) :::info 在 linux 的實作中,如果不考慮 `gffunc` 的使用,則 α 固定 = 2。輸入的 `prim` 是對此的 index,因此實際的 genearator $\beta$ = $\alpha^{prim}$,於是 $\beta^2 = \alpha^{2 * prim}$... 以此類推。 為了方便解說,後面我們也會把這裡設定的 primitive element 稱為 $\beta$ ::: * `iprim`: ?? * `gfpoly`: [Primitive polynomial](https://en.wikipedia.org/wiki/Primitive_polynomial_(field_theory)) 的多項式,這是為了把多項式制約在 Galois field 中。以 `0x409 = 0b0100 0000 1001` 為例,表示 primitive polynomial 為 $x^{10}+x^3+1$ ### `init_rs` 想要使用 Reed–Solomon codes,首先我們需要先初始化相關的資料結構。 ```c static struct rs_control *rs_decoder; rs_decoder = init_rs (10, 0x409, 0, 1, 6); ``` 我們可以在 [linux/include/linux/rslib.h](https://github.com/torvalds/linux/blob/0fe5f9ca223573167c4c4156903d751d2c8e160e/include/linux/rslib.h) 找到其定義,首先先來看相關的參數設定: ```c static inline struct rs_control *init_rs(int symsize, int gfpoly, int fcr, int prim, int nroots) { return init_rs_gfp(symsize, gfpoly, fcr, prim, nroots, GFP_KERNEL); } struct rs_control *init_rs_gfp(int symsize, int gfpoly, int fcr, int prim, int nroots, gfp_t gfp) { return init_rs_internal(symsize, gfpoly, NULL, fcr, prim, nroots, gfp); } ``` `init_rs(int symsize, int gfpoly, int fcr, int prim, int nroots)` * `symsize`: 每個 symbol 的位元數量 :::warning 根據 [Reed-Solomon Library Programming Interface](https://www.kernel.org/doc/html/latest/core-api/librs.html) > The databytes are expanded to the given symbol size on the fly. There is no support for decoding continuous bitstreams with a symbolsize != 8 at the moment. If it is necessary it should be not a big deal to implement such functionality. ::: * 其他的參數基本上就對應資料結構 `rs_codec` 中的同名稱 `init_rs` 會呼叫 `init_rs_gfp`,`init_rs_gfp` 再呼叫 `init_rs_internal`,針對額外的參數說明: ```c init_rs_internal(int symsize, int gfpoly, int (*gffunc)(int), int fcr, int prim, int nroots, gfp_t gfp) ``` * `gffunc`: 在 [non-canonical](https://en.wikipedia.org/wiki/Canonical_form) 模式下生成不同的 $\alpha$,不過在 `init_rs` 下會被設為 NULL(對此的使用還沒概念,暫時先忽略qq) * `gfp`: 使用 [kmalloc](https://www.kernel.org/doc/htmldocs/kernel-api/API-kmalloc.html) 時的 flag,會影響記憶體配置的區域和行為不同 ### `init_rs_internal` 初始 `rs_control` 結構。 ```c if (symsize < 1) return NULL; if (fcr < 0 || fcr >= (1<<symsize)) return NULL; if (prim <= 0 || prim >= (1<<symsize)) return NULL; if (nroots < 0 || nroots >= (1<<symsize)) return NULL; ``` 幾個基本的檢查被進行: * `symsize` 是每個 symbol 的位元數量,因此不能小於 1 * `fcr` 的範圍,參見前面的說明 * `prim` 的範圍,因為 primitive element 不能是 $\alpha^0 = 1$ 也不能是 GF 外的 $\alpha^i$ * 每個 block 只會有 `(1<<symsize)-1` 個 symbol,`nroots` 大於這個範圍是不合理的,小於 0 顯然也是 ```c /* * The decoder needs buffers in each control struct instance to * avoid variable size or large fixed size allocations on * stack. Size the buffers to arrays of [nroots + 1]. */ bsize = sizeof(uint16_t) * RS_DECODE_NUM_BUFFERS * (nroots + 1); rs = kzalloc(sizeof(*rs) + bsize, gfp); if (!rs) return NULL; ``` * decode 的時候,會需要一個固定的大小的 heap buffer 空間,可以看到 rs_control 中使用了 [Arrays of Length Zero](https://gcc.gnu.org/onlinedocs/gcc/Zero-Length.html) 技巧,因此可以在執行時期才決定分配的空間。`RS_DECODE_NUM_BUFFERS` = 8。 ```c mutex_lock(&rslistlock); /* Walk through the list and look for a matching entry */ list_for_each(tmp, &codec_list) { struct rs_codec *cd = list_entry(tmp, struct rs_codec, list); if (symsize != cd->mm) continue; if (gfpoly != cd->gfpoly) continue; if (gffunc != cd->gffunc) continue; if (fcr != cd->fcr) continue; if (prim != cd->prim) continue; if (nroots != cd->nroots) continue; /* We have a matching one already */ cd->users++; rs->codec = cd; goto out; } ``` * `list_for_each(tmp, &codec_list)` 逐一走訪整個鏈結串列,檢查是否已經存在相同的編碼結構,如果有,則直接 reference 一個相同的即可,不需要另外 init 一份 ```c /* Create a new one */ rs->codec = codec_init(symsize, gfpoly, gffunc, fcr, prim, nroots, gfp); if (!rs->codec) { kfree(rs); rs = NULL; } out: mutex_unlock(&rslistlock); return rs; ``` * 如果此結構不存在,則呼叫 `codec_init` 去建立一個新的結構 ### `codec_init` ```c int i, j, sr, root, iprim; struct rs_codec *rs; rs = kzalloc(sizeof(*rs), gfp); if (!rs) return NULL; INIT_LIST_HEAD(&rs->list); rs->mm = symsize; rs->nn = (1 << symsize) - 1; rs->fcr = fcr; rs->prim = prim; rs->nroots = nroots; rs->gfpoly = gfpoly; rs->gffunc = gffunc; /* Allocate the arrays */ rs->alpha_to = kmalloc_array(rs->nn + 1, sizeof(uint16_t), gfp); if (rs->alpha_to == NULL) goto err; rs->index_of = kmalloc_array(rs->nn + 1, sizeof(uint16_t), gfp); if (rs->index_of == NULL) goto err; rs->genpoly = kmalloc_array(rs->nroots + 1, sizeof(uint16_t), gfp); if(rs->genpoly == NULL) goto err; ``` * 初始化一個 `rs_codec` 所需的空間,並透過 `INIT_LIST_HEAD` 初始化這個鏈結串列節點 * 初始化相關成員、矩陣,關於這些成員的作用請參考前面 ```c /* Generate Galois field lookup tables */ rs->index_of[0] = rs->nn; /* log(zero) = -inf */ rs->alpha_to[rs->nn] = 0; /* alpha**-inf = 0 */ if (gfpoly) { sr = 1; for (i = 0; i < rs->nn; i++) { rs->index_of[sr] = i; rs->alpha_to[i] = sr; sr <<= 1; if (sr & (1 << symsize)) sr ^= gfpoly; sr &= rs->nn; } } else { /* else case 在 init_rs_non_canonical 才會用到, * 尚未研究,請先讓我省略這部份 qq */ ``` * 定義 `rs->alpha_to[rs->nn] = 0` * 接著,我們需要初始化 $\alpha^0$、$\alpha^1$、...、$\alpha^{2^{symsize}-1}$ 對應的數值,以及反過來的查找,這個加速我們之後編碼與解碼的運算 * 因為 $\alpha$ = 2,因此每次對 $\alpha^i$ 的運算就是把 $\alpha^{i-1}$ 左移 1 bit * `if (sr & (1 << symsize))` 判斷是否有 overflow,如果有,則需要額外對 `gfpoly`,也就是 primitive polynomial 做 XOR 運算,最後則 `sr &= rs->nn` 以確保得到的數字仍在 GF 中 ```c /* If it's not primitive, exit */ if(sr != rs->alpha_to[0]) goto err; /* Find prim-th root of 1, used in decoding */ for(iprim = 1; (iprim % prim) != 0; iprim += rs->nn); /* prim-th root of 1, index form */ rs->iprim = iprim / prim; ``` * 如果 sr 不等於 $\alpha^0$,表示 $\alpha^i$ 建立起的集合並非 GF,有錯誤 * 這裡要計算出 decode 用到的 `iprim` (??) ```c /* Form RS code generator polynomial from its roots */ rs->genpoly[0] = 1; for (i = 0, root = fcr * prim; i < nroots; i++, root += prim) { rs->genpoly[i + 1] = 1; /* Multiply rs->genpoly[] by @**(root + x) */ for (j = i; j > 0; j--) { if (rs->genpoly[j] != 0) { rs->genpoly[j] = rs->genpoly[j -1] ^ rs->alpha_to[rs_modnn(rs, rs->index_of[rs->genpoly[j]] + root)]; } else rs->genpoly[j] = rs->genpoly[j - 1]; } /* rs->genpoly[0] can never be zero */ rs->genpoly[0] = rs->alpha_to[rs_modnn(rs, rs->index_of[rs->genpoly[0]] + root)]; } ``` * genpoly 的 index 從 0 到 rs->nroots,就對應多項式從 1 到 $x^{nroots-1}$ * 建立生成多項式,我們知道生成多項式長成 $g(x) = \prod_{i=fcr}^{fcr+nroots} (x-\beta^i)$,這裡顯然是在對其進行計算 * 對於 `for (j = i; j > 0; j--)` 迴圈中的運算,可以理解為例如當下的 genpoly 對應的多項式可能是 $a_j+ a_{j-1}x^1...+ a_2x^{j-2} + a_1x^{j-1}$,然後要乘上下一個 $(x-\beta^{root})$,因此新的多項式係數,假設是 $x$ 項好了,就會變成 $a_{j-1} \times \beta^{root} + a_j \times 1$ * 呈上,$a_{j-1}$ 其實就是某一個 $\alpha^i$,而 $\beta^{root}$ 也是某個 $\alpha^j$,則因為 $\alpha^{(i+j)}$ 可能會超出能 index 到的範圍,因此就要去使用 `rs_modnn(rs,rs->index_of[rs->genpoly[j]] + root)` 去調整 ```c for (i = 0; i <= nroots; i++) rs->genpoly[i] = rs->index_of[rs->genpoly[i]]; rs->users = 1; list_add(&rs->list, &codec_list); return rs; err: kfree(rs->genpoly); kfree(rs->index_of); kfree(rs->alpha_to); kfree(rs); return NULL; ``` * 前面計算出來的 `genpoly` 是係數的真實值,為了方便計算,透過 `index_of` 改成用 index 來代表 * `rs->users = 1` 計算在 reference 這個結構的數量 ### `rs_modnn` ```c static inline int rs_modnn(struct rs_codec *rs, int x) { while (x >= rs->nn) { x -= rs->nn; x = (x >> rs->mm) + (x & rs->nn); } return x; } ``` * 如果 x 大於 `rs->nn = (1<<mm)-1`,表示 overflow 所以不存在定義的 GF 中,此時需調整成 `x % rs->nn` * 這裡對取餘運算進行了優化,可以參考[Mersenne Numbers: Mod 3, Mod 7, Mod ($2^i$ – 1)](http://homepage.divms.uiowa.edu/~jones/bcd/mod.shtml#exmod3) (沒有找到比較明確的數學證明,待補) ### `encode_rs8` 初始化 RS 相關的資料結構後,接著我們就可以來進行編碼了,下面以 `encode_rs8` 為範例 ```c /* Parity buffer. Size = number of roots */ uint16_t par[6]; /* Initialize the parity buffer */ memset(par, 0, sizeof(par)); /* Encode 512 byte in data8. Store parity in buffer par */ encode_rs8 (rs_decoder, data8, 512, par, 0); ``` ```c int encode_rs8(struct rs_control *rsc, uint8_t *data, int len, uint16_t *par, uint16_t invmsk) ``` * `rsc` 就是我們初始化過的結構 * `data` 是要被 encode 的資料,其長度為 `len` * `par` 則是呈裝 parity data 的地方,注意需要由使用者自行初始化 * `invmsk` 可以讓我們在編碼時加上一個 mask(關於為什麼需要 mask,可以稍微參考 [Reed–Solomon codes for coders: Masking](https://en.wikiversity.org/wiki/Reed%E2%80%93Solomon_codes_for_coders#Masking)) ```c pad = nn - nroots - len; if (pad < 0 || pad >= nn) return -ERANGE; ``` * 首先檢查 `nn - nroots - len`,也就是檢查資料數量加上 parity 數量是否超過總 symbol 數量,這是因為編碼數量必須小於 GF 的元素總量 $\lvert F\rvert$ ```c= for (i = 0; i < len; i++) { fb = index_of[((((uint16_t) data[i])^invmsk) & msk) ^ par[0]]; /* feedback term is non-zero */ if (fb != nn) { for (j = 1; j < nroots; j++) { par[j] ^= alpha_to[rs_modnn(rs, fb + genpoly[nroots - j])]; } } /* Shift */ memmove(&par[0], &par[1], sizeof(uint16_t) * (nroots - 1)); if (fb != nn) { par[nroots - 1] = alpha_to[rs_modnn(rs, fb + genpoly[0])]; } else { par[nroots - 1] = 0; } } ``` 以編碼而言,事實上就是把資料對生成多項式做長除法,並把餘數紀錄到 `par` 中,為了儘量說清楚,這裡我借用 [Polynomial_division](https://en.wikiversity.org/wiki/Reed%E2%80%93Solomon_codes_for_coders#Polynomial_division) 為例來解說,以這個長除法為例,假設我們的資料是 `12 34 56`,生成多項式為 `01 0f 36 78 40`,然後其 `par` 的大小為 4 個 uint8: ``` 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 ``` * 先假設沒有設 invmsk * 那麼最外層迴圈的第 1 層,i = 0,一開始的 `par` 成員就是 `00 00 00 00` ,`fb` 先取出 `12` * 在第 5 行的 for 迴圈要計算出後面的 `ee 2b 23`,這會是透過 `par[j] ^= alpha_to[rs_modnn(rs, fb + genpoly[nroots - j])]` 來得到 * 我們要計算 `12` * `0f 36 78` 為何。在這裡 `12` 就是 $\alpha^{fb}$,而 generator 的 `0f 36 78` 會對應某個 $\alpha^{genpoly[nroots - j]}$ * 因此我們要計算的 `12` * `0f 36 78` 對應每個 $\alpha^{fb+{genpoly[nroots - j]}}$ ,也就是程式中的 `alpha_to[rs_modnn(rs, fb + genpoly[nroots - j])];` * 至此的 par 是 `12 ee 2b 23`,接著我們把它向左位移 (也就是這裡的 `memmove`),變成 `ee 2b 23 23` (雖然這裡我說位移,不過最後一個 byte 沒有變成 0 ,這是因為事實上是使用 memmove,不過稍候它就會被覆蓋掉) * 第 13 行的概念如前,因此 par 變成 `ee 2b 23 f4` * 來到最外層迴圈的第二層,i = 1,fb 會是 `34 ^ ee = da`,然後就依尋前面的步驟,一直做到處理完整個 data * 最後得到的餘數會在 parity 中,那就是編碼的 parity ### `decode_rs8` 最複雜的部份應該就是這裡了。對於 RS 的 decode,可以有許多不同的演算法被應用。不過我們可以大致總結成 5 點: 1. 計算出 syndromes polynomial,透過 Berlekamp-Massey 等演算法檢查資料是否有錯誤產生 2. 根據 syndromes 計算出 erasure/error locator polynomial,用 Berlekamp-Massey 找到確切是哪些位置產生了錯誤 3. 根據前兩者算出 erasure/error evaluator polynomial 有算出有多少 symbol 被竄改 4. 由前3個多項式計算 erasure/error magnitude polynomial, 5. 你可想像成這個 polynomial 就是揪出來的 noise,則我們把 noise 從接收到的資料後去除後,就可得到原始的資料 ```c int decode_rs8(struct rs_control *rsc, uint8_t *data, uint16_t *par, int len, uint16_t *s, int no_eras, int *eras_pos, uint16_t invmsk, uint16_t *corr) ``` * `rsc` 就是我們初始化好的 `rs_control` 結構 * 輸入長度為 `len` 的 `data` 跟 `par` * `s`: 解碼的過程中會需要 syndromes polynomial,找出發生錯誤的位置,因為 syndromes 的計算可能可以用硬體來計算,因此 API 提供此參數可以直接依據這個 syndromes 來修正錯誤 * 下面的範例先假設這個欄位為 NULL,也就是需要軟體負責計算 syndromes * `no_eras`、`eras_pos`: 有些 channel 中,錯誤會發生的情形是可以預知的,這種錯誤被稱做 erasure。要修正一個 erasure,只要添加一個 symbol 就可以做 [erasure correction](https://en.wikiversity.org/wiki/Reed%E2%80%93Solomon_codes_for_coders#Erasure_correction) 了 * 因此,我們可以告知 decoder 相關的訊息來做 erasure 修正 * 如果不知道錯誤的位置,那就稱為 error,這種情況下,就需要添加 2 個 symbol 才可以做修正 * `invmsk` 就等於我們在 encode 時所使用的 * `corr`: 根據 `eras_pos` 可以要做 correction 相關的 bitmsk * 所以這個 API 可以不輸入 data 和 parity,僅輸入 `s`、`eras_pos` 等得到 `corr` 的 mask,再透過自定義的程式去處理問題 ```c enum { RS_DECODE_LAMBDA, RS_DECODE_SYN, RS_DECODE_B, RS_DECODE_T, RS_DECODE_OMEGA, RS_DECODE_ROOT, RS_DECODE_REG, RS_DECODE_LOC, RS_DECODE_NUM_BUFFERS }; ... /* * The decoder buffers are in the rs control struct. They are * arrays sized [nroots + 1] */ uint16_t *lambda = rsc->buffers + RS_DECODE_LAMBDA * (nroots + 1); uint16_t *syn = rsc->buffers + RS_DECODE_SYN * (nroots + 1); uint16_t *b = rsc->buffers + RS_DECODE_B * (nroots + 1); uint16_t *t = rsc->buffers + RS_DECODE_T * (nroots + 1); uint16_t *omega = rsc->buffers + RS_DECODE_OMEGA * (nroots + 1); uint16_t *root = rsc->buffers + RS_DECODE_ROOT * (nroots + 1); uint16_t *reg = rsc->buffers + RS_DECODE_REG * (nroots + 1); uint16_t *loc = rsc->buffers + RS_DECODE_LOC * (nroots + 1); /* Check length parameter for validity */ pad = nn - nroots - len; BUG_ON(pad < 0 || pad >= nn - nroots); ``` * 首先,還記得我們預先拿到的 `rs->buffers` 嗎? 當時我們 allocate 了 `8 * (nroots + 1)` uint16_t 大小的空間,現在我們把它切在 8 等份給不同的需求使用 * encode 時相同的理由,做對應的檢查 ```c /* Does the caller provide the syndrome ? */ if (s != NULL) { for (i = 0; i < nroots; i++) { /* The syndrome is in index form, * so nn represents zero */ if (s[i] != nn) goto decode; } /* syndrome is zero, no errors to correct */ return 0; } ``` * 如果 caller 有輸入 syndrome 參數的話,檢查 syndrome 是否為 0 * 因為 syndrome 是 index 的,而我們在初始化時也定義了 `rs->alpha_to[rs->nn] = 0`,所以這邊只要檢查 syndrome 是否等於 `nn` 即可 ```c /* form the syndromes; i.e., evaluate data(x) at roots of * g(x) */ for (i = 0; i < nroots; i++) syn[i] = (((uint16_t) data[0]) ^ invmsk) & msk; for (j = 1; j < len; j++) { for (i = 0; i < nroots; i++) { if (syn[i] == 0) { syn[i] = (((uint16_t) data[j]) ^ invmsk) & msk; } else { syn[i] = ((((uint16_t) data[j]) ^ invmsk) & msk) ^ alpha_to[rs_modnn(rs, index_of[syn[i]] + (fcr + i) * prim)]; } } } for (j = 0; j < nroots; j++) { for (i = 0; i < nroots; i++) { if (syn[i] == 0) { syn[i] = ((uint16_t) par[j]) & msk; } else { syn[i] = (((uint16_t) par[j]) & msk) ^ alpha_to[rs_modnn(rs, index_of[syn[i]] + (fcr+i)*prim)]; } } } s = syn; /* Convert syndromes to index form, checking for nonzero condition */ syn_error = 0; for (i = 0; i < nroots; i++) { syn_error |= s[i]; s[i] = index_of[s[i]]; } if (!syn_error) { /* if syndrome is zero, data[] is a codeword and there are no * errors to correct. So return data[] unmodified */ return 0; } ``` * 如果沒有提供 `s` 參數,我們就要自行計算 syndromes polynomial * syndromes polynomial 的計算就是把 data(包含資料部份和 parity) 當成多項式,然後代入 $\beta^{fcr}$ 到 $\beta^{fcr+nroots}$ 看看是否為 0 * 這邊有使用 [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method) 避免直接計算 $x^n$ 做優化 * 把計算完的轉換為 index 形式方便後續的處理 * `syn_error` 計算是不是所有 syndrome 都是 0,如果不是就要去修正 ```c decode: memset(&lambda[1], 0, nroots * sizeof(lambda[0])); lambda[0] = 1; if (no_eras > 0) { /* Init lambda to be the erasure locator polynomial */ lambda[1] = alpha_to[rs_modnn(rs, prim * (nn - 1 - (eras_pos[0] + pad)))]; for (i = 1; i < no_eras; i++) { u = rs_modnn(rs, prim * (nn - 1 - (eras_pos[i] + pad))); for (j = i + 1; j > 0; j--) { tmp = index_of[lambda[j - 1]]; if (tmp != nn) { lambda[j] ^= alpha_to[rs_modnn(rs, u + tmp)]; } } } } for (i = 0; i < nroots + 1; i++) b[i] = index_of[lambda[i]]; ``` * 如果已知哪些位置出錯,如前所述要做 erasure correction * 首先初始 $\Lambda(x)$ = 1 * 如果有 erasure,則 要計算 erasure locator polynomial $\Lambda(x) = \prod (1-\beta^ix)$ * 其中 i 的值是根據錯誤位置計算 * 要注意到例如我們的資料是 `01 02 03 04`,對應的多項式是 $01x^3 + 02x^2 + 03x + 04$,index 0 的地方是多項式的最高次方,但是如果錯誤是發生在 `03`(index = 2) 的地方,要代入的是 $\beta^1$ ,`nn - 1 - (eras_pos[0] + pad)` 就是在做此計算 * 則迴圈的計算邏輯就和前面的 g(x) 計算大致相同 * `b` 對應 Berlekamp-Massey algorithm 的多項式 B(x),用 lamda 的 $\alpha$ index 初始化 ```c /* * Begin Berlekamp-Massey algorithm to determine error+erasure * locator polynomial */ r = no_eras; el = no_eras; while (++r <= nroots) { /* r is the step number */ /* Compute discrepancy at the r-th step in poly-form */ discr_r = 0; for (i = 0; i < r; i++) { if ((lambda[i] != 0) && (s[r - i - 1] != nn)) { discr_r ^= alpha_to[rs_modnn(rs, index_of[lambda[i]] + s[r - i - 1])]; } } discr_r = index_of[discr_r]; /* Index form */ ... ``` * 我們需要透過 [Berlekamp-Massey algorithm](https://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm) 來找出未知位置的錯誤(error locator polynomial) * 迴圈要進行 `nroots` 減去已知的錯誤次(換句話說,如果沒有 erasure 的話就是進行 nroots 次) * 首先要計算每一步的 discrepancy,這個數字等於 $\sum_{i=0}^{r}\Lambda_iS_{r-i-1}$ ```c ... if (discr_r == nn) { /* 2 lines below: B(x) <-- x*B(x) */ memmove (&b[1], b, nroots * sizeof (b[0])); b[0] = nn; } else { /* 7 lines below: T(x) <-- lambda(x)-discr_r*x*b(x) */ t[0] = lambda[0]; for (i = 0; i < nroots; i++) { if (b[i] != nn) { t[i + 1] = lambda[i + 1] ^ alpha_to[rs_modnn(rs, discr_r + b[i])]; } else t[i + 1] = lambda[i + 1]; } if (2 * el <= r + no_eras - 1) { el = r + no_eras - el; /* * 2 lines below: B(x) <-- inv(discr_r) * * lambda(x) */ for (i = 0; i <= nroots; i++) { b[i] = (lambda[i] == 0) ? nn : rs_modnn(rs, index_of[lambda[i]] - discr_r + nn); } } else { /* 2 lines below: B(x) <-- x*B(x) */ memmove(&b[1], b, nroots * sizeof(b[0])); b[0] = nn; } memcpy(lambda, t, (nroots + 1) * sizeof(t[0])); } } ``` * `discr_r == nn` 表示原求的 discrepancy `d` 為 0 * 此時 $B(x) \Longleftarrow x*B(x)$ * 然後前進到下個 iteration * 如果 discrepancy 不為 0: * $\Lambda'(x) \Longleftarrow \Lambda(x) - d \times x \times B(x)$ * 可以看到程式中的 `t` 就是 $\Lambda(x) - d \times x \times B(x)$ 的暫存空間 $\Lambda'(x)$ * 每個 iteration 中的 `r + no_eras` 表示目前是第幾個 iteration (從 1 開始數) * `el` 的用途是用來判斷某些 iteration 中需要額外的調整 `B(x)`,在條件滿足的 `el` 下: * $el \Longleftarrow$ 現在的 iteration 數量 $- \space el$ * 此時 $B(x) \Longleftarrow d^{-1}\Lambda(x)$ * 除此之外 $B(x) \Longleftarrow x*B(x)$ * 最後,$\Lambda(x) \Longleftarrow \Lambda'(x)$,然後進入下一個 iteration ```c /* Convert lambda to index form and compute deg(lambda(x)) */ deg_lambda = 0; for (i = 0; i < nroots + 1; i++) { lambda[i] = index_of[lambda[i]]; if (lambda[i] != nn) deg_lambda = i; } if (deg_lambda == 0) { /* * deg(lambda) is zero even though the syndrome is non-zero * => uncorrectable error detected */ return -EBADMSG; } ``` * 最終,我們會得到一個 erasure /error locator polynomial $\Lambda(x)$ * 把 lambda 的值轉換成 index 形式,並且計算是幾階的多項式 `deg_lambda` * 如果結果為 0 階但是 syndrome 卻是非零值,那表示此錯誤不能被修正 ```c /* Find roots of error+erasure locator polynomial by Chien search */ memcpy(&reg[1], &lambda[1], nroots * sizeof(reg[0])); count = 0; /* Number of roots of lambda(x) */ for (i = 1, k = iprim - 1; i <= nn; i++, k = rs_modnn(rs, k + iprim)) { q = 1; /* lambda[0] is always 0 */ for (j = deg_lambda; j > 0; j--) { if (reg[j] != nn) { reg[j] = rs_modnn(rs, reg[j] + j); q ^= alpha_to[reg[j]]; } } if (q != 0) continue; /* Not a root */ if (k < pad) { /* Impossible error location. Uncorrectable error. */ return -EBADMSG; } /* store root (index-form) and error location number */ root[count] = i; loc[count] = k; /* If we've already found max possible roots, * abort the search to save time */ if (++count == deg_lambda) break; } ``` * 則對於 erasure /error locator polynomial $\Lambda(x)$,其根($\beta^i$ 使得 $\Lambda(x)$ = 1) 就表示錯誤發生的地方是 len(msg) - 1 - i * 這裡使用了 [Chien search](https://en.wikipedia.org/wiki/Chien_search) 加速求根 * ??? ```c if (deg_lambda != count) { /* * deg(lambda) unequal to number of roots => uncorrectable * error detected */ return -EBADMSG; } /* * Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo * x**nroots). in index form. Also find deg(omega). */ deg_omega = deg_lambda - 1; for (i = 0; i <= deg_omega; i++) { tmp = 0; for (j = i; j >= 0; j--) { if ((s[i - j] != nn) && (lambda[j] != nn)) tmp ^= alpha_to[rs_modnn(rs, s[i - j] + lambda[j])]; } omega[i] = index_of[tmp]; } ``` * 如果找到的根數量和 $\Lambda$ 的階數不一致,表示錯誤無法被修正 * 接下來要求 error / evaluator polynomial $\Omega(x) = S(x) * \Lambda(x) \space mod \space x^{nroots}$,且存成 index form ```c /* * Compute error values in poly-form. num1 = omega(inv(X(l))), num2 = * inv(X(l))**(fcr-1) and den = lambda_pr(inv(X(l))) all in poly-form * Note: we reuse the buffer for b to store the correction pattern */ num_corrected = 0; for (j = count - 1; j >= 0; j--) { num1 = 0; for (i = deg_omega; i >= 0; i--) { if (omega[i] != nn) num1 ^= alpha_to[rs_modnn(rs, omega[i] + i * root[j])]; } if (num1 == 0) { /* Nothing to correct at this position */ b[j] = 0; continue; } num2 = alpha_to[rs_modnn(rs, root[j] * (fcr - 1) + nn)]; den = 0; /* lambda[i+1] for i even is the formal derivative * lambda_pr of lambda[i] */ for (i = min(deg_lambda, nroots - 1) & ~1; i >= 0; i -= 2) { if (lambda[i + 1] != nn) { den ^= alpha_to[rs_modnn(rs, lambda[i + 1] + i * root[j])]; } } b[j] = alpha_to[rs_modnn(rs, index_of[num1] + index_of[num2] + nn - index_of[den])]; num_corrected++; } ``` * [Forney algorithm](https://en.wikipedia.org/wiki/Forney_algorithm) ```c /* * We compute the syndrome of the 'error' and check that it matches * the syndrome of the received word */ for (i = 0; i < nroots; i++) { tmp = 0; for (j = 0; j < count; j++) { if (b[j] == 0) continue; k = (fcr + i) * prim * (nn-loc[j]-1); tmp ^= alpha_to[rs_modnn(rs, index_of[b[j]] + k)]; } if (tmp != alpha_to[s[i]]) return -EBADMSG; } ``` * 經過前面的步驟後,這裡我們會得到誤差多項式 e(x)(重用 `b` 來保存係數) $$ e(x) = e_0 + e_1x + ... + e_{n-1}x^{(n-1)} $$ * 理論上,我們前面得到的 syndromes $S_i$ 應該要等於 $e(\beta^i)$,這裡再做額外的檢查以確保正確性 ```c /* * Store the error correction pattern, if a * correction buffer is available */ if (corr && eras_pos) { j = 0; for (i = 0; i < count; i++) { if (b[i]) { corr[j] = b[i]; eras_pos[j++] = loc[i] - pad; } } } else if (data && par) { /* Apply error to data and parity */ for (i = 0; i < count; i++) { if (loc[i] < (nn - nroots)) data[loc[i] - pad] ^= b[i]; else par[loc[i] - pad - len] ^= b[i]; } } return num_corrected; ``` * 如果 input 是 `corr` 和 `eras_pos` 的話,就把該誤差多項式填入 `corr` 並把有錯的位置填入 `eras_pos`, 當成回傳 * 如果是 input 是 `data` 和 `par` 的話,那就是看誤差項的位置在 `data` 還是 `par`,並且將結果修正回去 ## 參考資料 * [錯誤更正碼簡介](https://w3.math.sinica.edu.tw/math_media/d184/18404.pdf) * [Hamming(7,4)](https://en.wikipedia.org/wiki/Hamming(7,4)) * [Reed-Solomon error-correcting code decoder](https://www.nayuki.io/page/reed-solomon-error-correcting-code-decoder#preliminaries) * [Reed–Solomon error correction(wikipedia)](https://en.wikipedia.org/wiki/Reed%E2%80%93Solomon_error_correction) * [Reed–Solomon codes for coders](https://en.wikiversity.org/wiki/Reed%E2%80%93Solomon_codes_for_coders)