# Error-Correcting Codes (ECC)
> 貢獻者: RinHizakura, [jouae](https://github.com/jouae)
本文先探討比較容易理解的 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) 在內演算法的視覺化解說。
:::success
Under working: 閱讀第一手資料。整理 [Claude Shannon](https://en.wikipedia.org/wiki/Claude_Shannon) 的 [A Mathematical Theory of Communications](https://en.wikipedia.org/wiki/A_Mathematical_Theory_of_Communication) 以中文敘述跟更易懂的例子解釋論文中的數學原理。並解釋當初基於什麼樣的公理與假設建立起
$$
H(p_1,p_2,\dots,p_n) = -\sum_{i=1}^n p_i\log (p_i).
$$
參閱 [閱讀筆記：通訊領域中的數學理論 by jouae](https://hackmd.io/@jouae/Shannon_entropy)
:::
## 為什麼會有錯誤
> Refs:
1. Error Correction Coding Mathematical Methods and Algorithms by Todd K. Moon
2. Error Control Coding: Fundamentals and Applications by Shu Lin
在 Todd 所著一書中描述：
一序列的數個位元可以藉由一對應關係寫成波形(waveform)的形式透過通道(channel)傳遞。二元相移鍵控(Binary phase-shift keying, BPSK)為其中一種調幅(amplitude modulation)形式，其中波形的正負號用於表示位元為 `0` 或是 `1`。
以 $b_i \in \lbrace 0, 1 \rbrace$ 構成的序列
$$
\lbrace \dots, b_{-2}, b_{-1}, b_0, b_1, b_2, \dots \rbrace
$$
表示傳遞的訊息。其中每個位元到達的時間需要 $T$ 單位時間，換言之，週期為 $T$ ，頻率則為 $\frac{1}{T}$。
要以波形的正負號來表示位元為 `0` 或 `1` ，需要先做一個變數變換。以 $\tilde{b_i}$ 表示 $\lbrace 0, 1 \rbrace$ 經過映射至 $\lbrace -1, 1 \rbrace$ 的結果，具體的表示為：
$$
\tilde{b_i}=(2b_i - 1) \text{ or } \tilde{b_i}=-(2b_i - 1)
$$
我們要確立該波形的數學表示式，在有該波形的週期 $T$ 時，再來要知道該波形的振福(amplitude)，所以藉由找出一位元對應振福的映射關係，來確定該波形的振福。
令振福 $a_i$ 表示為
$$
a_i= \sqrt{E_b}(2b_i-1) = -\sqrt{E_b}(-1)^{b_i} = \sqrt{E_b}\tilde{b_i}
$$
得到該振幅後，我們想要藉由一個坐標系來描述在一個訊號中什麼樣的位置會表示為位元為 `0` 或 `1`。先從我們熟知的一維坐標系來思考，一維坐標系就是數線。數線的構成要素為原點、方向，及單位度，二維坐標系就是在一維的基礎上，找出兩同原點且互相正交的單位度。所以在以訊號為坐標軸該怎麼表示單位度？所以要先定義單位訊號是什麼，更具體的說法是單位能量(unit-energy)訊號。
對一連續時間的訊號，可以藉由對時間積分的方式來描述在該段時間內的能量為多少，在《Signal and System》一書 1.1.2 節提及從伏特 $v(t)$ 與電流 $i(t)$ 流經一電阻 $R$ 的物理角度思考該電阻的順時功率為：
$$
p(t) = v(t)i(i) = \dfrac{1}{R} v(t)^2
$$
第二個等式藉由歐姆定律(Ohm's law) $v=iR$ 得來。
在自 $t_0$ 至 $t_1$ 時間內的功耗則可以表示為：
$$
\int_{t_0}^{t_1} p(t) dt= \dfrac{1}{R}\int_{t_0}^{t_1} v(t)^2 dt
$$
藉由該物理概念，將一訊號 $s(t)$ 的能量計算方式定義成自身對自身的 $L_2$ 內積
$$
E = \left< s(t), s(t) \right> = \int^{\infty}_{-\infty} s(t)^2 dt.
$$
其中，一單位能量訊號表示為：
$$
\psi_1(t) = \left< \psi_1(t), \psi_1(t) \right> =1.
$$
要注意的是該訊號 $\psi_1(t)$ 在區間 $[0,T)$ 有[支撐集(support)](https://en.wikipedia.org/wiki/Support_(mathematics))。簡單的想在該區間 $t\in[0,T)$ 中 $\psi_1(t)$ 有非零值， $t\in\mathbb{R}- [0,T)$ 中皆為零。
所以藉由該單位能量訊號 $\psi_1(t)$，我們可以假想一個一維訊號坐標系以 $\psi_1(t)$ 表示為一個軸上的單位訊號，該軸上的點表示不同的相位(phase)。這樣的空間稱之為訊號空間(signal space)。
所以該位元 $b_i$ 在時間點 $iT$ 到達的波形，能以一訊號表示：
$$
a_i\psi_1(t-iT)=\sqrt{E_b}\tilde{b_i}\psi_1(t-iT)
$$
其中， $b_i=0$ 對應到的訊號為 $-\sqrt{E_b}\psi_1(t)$，$b_i=1$ 對應到的訊號為 $\sqrt{E_b}\psi_1(t)$。在此以 $\psi_1(t)$ 為一維訊號空間維基底的空間，會有兩個座標點 $\sqrt{E_b}$ 和 $-\sqrt{E_b}$。
將此訊號空間與相位繪製成圖，該圖稱之為調變號分佈圖(constellation diagram)。
<!-- 需要重新繪製 1d_signal_space -->
顯而易見的傳送一位元 $b_i$ 需要的能量為：
$$
\left< a_i\psi_1(t), a_i\psi_1(t) \right> = E_b
$$
以此方式，可以將有限個位元以依連續訊號表示：
$$
s(t) = \sum_i a_i\psi_i(t-iT)
$$
這個手法讓我能夠將離散的資訊以連續的訊號傳送。再來討論更一般化的情況，在單一個訊號波形僅能傳遞一個位元的資訊，即 `0` 和 `1`。假設想要傳遞 `00`、`01`、 `10` 和 `11` 該怎麼以連續訊號表示？
## 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`
### 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.
$$
第一個方程對應上述文式圖紅圖部分，第二個方程對應藍圖部分，第三個方程對應綠圖部分。
將上述線性方程式寫成矩陣則是：
$$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}
$$
令 $C$ 為所有編碼資料的集合即是 H 的 null space，也就是 $\forall x \in C$ $Hx^T = 0$，又 $\text{rank}(H) = 3$，根據 [Rank–nullity theorem](https://en.wikipedia.org/wiki/Rank%E2%80%93nullity_theorem) 則 $\text{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
則我們可以定義一段 code word 的 minimum distance 為此 code word 集合中任兩個編碼的最小 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 個位元有錯的 code word 之後，得到的結果為:
$$
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)))
### 漢明邊界 hamming bound
> 參考自 Introduction to Abstract Algebra by W. Keith Nicholson - 2.11 An Application to Binary Linear Code
<!-- 此小節使用點集拓樸與集合來描述，試圖不引入任何矩陣解釋偵測到錯誤碼與可糾正錯誤碼的上界 -->
令 $B^n=\lbrace 0,1\rbrace \times \dots \times \lbrace 0,1\rbrace$ 為一 $n$ 摺[笛卡爾積](https://en.wikipedia.org/wiki/Cartesian_product)(n-fold cartesian product)，以及集合
$$
B_r(w) = \lbrace v \in B^n \vert d(w,v) < r\rbrace \\
$$
為以 $w$ 為中心，漢明距離 $d(\cdot,\cdot)$ 為測度(metric)，蒐集與 $w$ 距離 $r$ 以內的所有元素。$B_r(w)$ 也可以叫做以 $d(\cdot,\cdot)$ 為距離函數，圓心為 $w$ ，半徑為 $r$ 在 $B^n$ 中的 t-ball，在點集拓樸中該球又稱作一**開球(open ball)**。與之相對的為**閉球(closed ball)** 差別在於明漢距離可等於 $r$
$$
\bar{B_r}(w) = \lbrace v \in B^n \vert d(w,v) \leq r\rbrace
$$
我們可以用從點集拓樸的術語，來描述錯誤碼的偵測與糾正。以下節錄自《Introduction to Abstract Algebra by W. Keith Nicholson》第 147 頁：
假設 $c\in C$ 且 $C$ 最小距離為 $d$ ，換言之 $C$ 集合中任意兩不相同元素的集合距離不小於 $d$。更進一步的具象來講，對任意的 $c,c'\in C$ 且 $c\neq c'$ 為圓心，以半徑 $d$ 的開球(open ball)皆不會相交，用集合符號敘述則為 $B_d(c) \cap B_d(c') = \emptyset$。
假設一個 code word $c$ 被傳送後，接收者得到的 code word $w$ 帶有 $s$ 個錯誤，其中 $1\leq s \leq t$，藉由漢明距離可以得到 $s=d(c,w)$。因此，開球 $B_t(c)$ 在以漢明距離 $t$ 為半徑所圈起來的圓，可以收到最多 $t$ 個錯誤。
那在此我們要以集合論的方式說明兩個漢明碼的功能**偵測錯誤**與**校正錯誤**：
先假設 $C$ 該集合沒有任何 code word 落在以漢明距離 $t$ 為半徑，另一個 code word 的開球中，則
1. 偵測錯誤：
$w\in B_t(c)$ 且 $w\neq c$ 時， $w$ 為一個錯誤碼。
2. 校正錯誤：
偵測錯誤的條件下，如果再加上每個 code words 形成的半徑 $t$ 開球皆式兩兩不相交。則 $w$ 可以被校正成正確的 code word。
依據這個概念，我們**可以藉由漢明距離和開球等集合論方式寫出漢明碼偵測錯誤與校正錯誤的能力**。
* 下圖將集合 $C$ 與其中元素視覺化的結果，該圖顯示當 $C$ 集合中不存在任何錯誤碼的情況的長相。對所有 $c_i,c_j\in C$，$B_d(c_i)\cap B_d(c_j)\neq \emptyset$ 其中 $i\neq j$。
:warning: 值得注意的是，一般常用用語圓形與半徑距離是在歐式空間 $\mathbb{R}^2$ 底下以 $l_2$ 範數為距離定義而成的，漢明碼中使用的是漢明距離，實際的「圓形」與「半徑」與圖中不一定相符。此處為便於解釋以常用圖示表示。
> 不同距離函數下的圓形，可參考 [Unit circles with $l_p$-norm](https://en.wikipedia.org/wiki/Lp_space#/media/File:Vector-p-Norms_qtl1.svg)
- [ ] 假如 $C$ 為 $B^n$ 一非空子集合，同時 $C$ 的最小距離(minimum distance) 為 $d$，則：
1. 如果 $t+1 \leq d$ ， $C$ 可以偵測到 $t$ 個錯誤。
2. 如果 $2t+1 \leq d$ ， $C$ 可以校正 $t$ 個錯誤。
* 證明：
1. 對任意 $c\in C$ 的開球 $B_t(c)$ ，因為 $t<d$ ，僅有一元素圓心 $c$，意即 $B_t(c)=\lbrace c\rbrace$，該球不包含其他元素，換言之沒有任何錯誤碼出現在 $B_t(c)$ 中。
2. 藉由反證法。在 $2t+1\leq d$ 條件下，假設所有在 $C$ 中的元素以 $t$ 為半徑所構成的開球皆兩兩不相交，換言之 $C$ 無法校正錯誤。但集合 $C$ 中的元素 $c\neq c'$ 和 $w\in B_t(c)\cap B_t(c')$ ，藉由三角不等式：
$$
d(c,c') \leq d(c,w) + d(w,c') \leq t+t = 2t < d
$$
則得 $d(c,c') < d$ 與 $C$ 的最短距離為 $d$ 矛盾。所以與假設「所有在 $C$ 中的元素以 $t$ 為半徑所構成的開球皆兩兩不相交」不成立，故以 $t$ 為半徑所構成的開球，至少會有一對開球會相交。
再來我們有興趣的是，開球 ${B_r}(w)$ 究竟有多少個元素，換言之該集合的[基數(cardinal number)](https://en.wikipedia.org/wiki/Cardinal_number)為多少？藉由漢明距離的定義， $d(w,v)=r$ 成立時，若且為若 $w,v$ 兩者有 $r$ 個不同的位元。換句話說對於一個 $n$ 位元長的元素，選取 $r$ 個不同位置為 $w,v$ 位元差異處，寫成組合數可以表示成 ${n \choose r}$。藉此，該集合的基數可以表示成：
$$
\vert B_r(w) \vert =
{n \choose 0} + \dots + {n \choose r}
$$
- [ ] 漢明邊界：假如 $C$ 為 $B^n$ 一非空子集合，可以校正 $t$ 個錯誤，則
$$
\vert C \vert \leq \dfrac{2^n}{{n \choose 0} + \dots + {n \choose t}}
$$
* 證明：
令 $N=\vert B_t(w) \vert =
{n \choose 0} + \dots + {n \choose t}$。對不同 $c\in C$ 每個開球有 $N$ 個 code words，又開球數量總共有 $\vert C \vert$ 個(因為有 $\vert C \vert$ 個圓心)。如果兩兩開球皆不相交，則有 $N\vert C \vert$ 獨特的 code words，且 $\vert B^n \vert$，因此
$$
N\vert C \vert \leq 2^n
$$
移項後得到 $\vert C \vert \leq \dfrac{2^n}{{n \choose 0} + \dots + {n \choose t}}$
:::success
TODO: 如何從機率的觀點解釋 hamming code
TODO: ECC 校正錯誤的能力要如何量化？
可以參考：
* shu lin and daniel j. costello, "error control coding: fundamentals and applications," pearson/ prentice hall, second edition, 2004.
:::
## 基礎代數
### 二元運算子 binary operation
在一個非空集合 $G$ 定義的二元運算 $*$ 為一函數 $f$ （或稱為映射），其映射關係為該集合的[笛卡兒積](https://en.wikipedia.org/wiki/Cartesian_product)至該集合本身 $f: G\times G \rightarrow G$ ，另一種更常見的寫法為 $x*y$ ，其中 $x,y$ 屬於 $G$。
同時該二元運算滿足：
- [ ] 封閉性：對任意元素 $x,y\in G$， $x*y\in G$
* 假如任意 $x,y,z\in G$ 其二元運算 $*$ 滿足：
$$
(x*y)*z=x*(y*z)
$$
則稱該二元運算子具有**結合律**。
* 假如任意 $x,y\in G$ 其二元運算 $*$ 滿足：
$$
x*y=y*x
$$
則稱該二元運算具有**交換律**。
### 群 group
當一非空集合 $G$ 中的二元運算 $*$ 滿足以下性質時稱之為群：
1. 該二元運算子具有結合律。
2. $G$ 存在一單位(identity)元素 $e$ ，使得對任意 $x\in G$ 滿足$a*e=e*a=a$
3. 對任意 $x\in G$ 存在一反元素 $x'$ ，使得 $x*x'=e$
例子， $G=\lbrace 0, 1 \rbrace$ 與一運算 $*$ ，有以下映射關係：
$$
1*1 = 0, \quad 0*0 = 0, \quad 1*0 = 0
$$
該集合的運算子是否為二元運算？假如是，該集合是否為有一二元運算 $*$ 的群？
顯而易見的該運算為二元運算，對任意元素的運算結果皆是集合中的元素。同時， $0$ 為一單位元素； $1$ 的反元素為自己本身； $0$ 的反元素為自己本身。結合律可以藉由列舉的方式說明，在此不贅述。所以 $(G,*)$ 為一群。
假如對任意 $x,y\in G$ 該二元運算 $*$ 滿足交換律 $x*y=y*x$ 則稱為**阿貝爾群(abelian group)**。顯而易見的上述例子 $(G,*)$ 為一阿貝爾群。
## 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(®[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)
Sign in with Wallet
Connect another wallet