有限体のプログラミング後編: リードソロモン符号のしくみ

[前編](https://qiita.com/bellbind/items/5636c99be7ed8d7eaaf1)内容

1. 素数べき要素数の有限体のための基礎知識
2. 素数べき要素数の有限体とその作り方
3. 素数べき要素数の有限体の例
4. 素数べき要素数の有限体のガロア群

[中編](https://qiita.com/bellbind/items/5c977cc945a6ff4a8853)内容

5. JavaScriptで実装する素数べき要素数有限体
6. 2べき要素数の有限体に特化した実装

[後編](https://qiita.com/bellbind/items/b85fe5f5055fc9d42472)内容

7. 有限体の応用: リードソロモン符号とBCH符号

---

# 注: 記事内のJavaScriptコードについて

この記事内のコードで用いるJavaScriptの仕様は、`Array.flatMap()`のある[ECMAScript2019仕様](https://www.ecma-international.org/ecma-262/10.0/index.html)までを用います。JSファイルは、すべて拡張子は`.js`のES moduleです。

非ECMAScriptな要素として、`console.log`と`console.assert`のみを用います。その他のWeb APIやnode.jsビルトイン等は一切使用していません。

この記事内のコードは、以下のgistに置いています。

- https://gist.github.com/bellbind/ba4b307fea8251f7f40aeb7e9349814b

多項式や有限体等の要素は、数値と`Array`にマッピングし、それらのデータ処理を`function`式の中で記述します。体や多項式等は、この`function`実装をメンバーに持つ`Object`として提供します。たとえば、有理数体Qでの加算は以下のようになります。

- $c = \frac{1}{2} \times -\frac{3}{4}$ :  `const c = Q.mul([1, 2], [-3, 4]);` 


# リードソロモン符号とは

有限体の応用として、エラー検出だけでなく**エラー訂正**も可能な[リードソロモン符号](https://ja.wikipedia.org/wiki/%E3%83%AA%E3%83%BC%E3%83%89%E3%83%BB%E3%82%BD%E3%83%AD%E3%83%A2%E3%83%B3%E7%AC%A6%E5%8F%B7)があります。

このリードソロモン符号は[CD](https://en.wikipedia.org/wiki/Cross-interleaved_Reed%E2%80%93Solomon_coding)や[MPEG2](https://www.etsi.org/deliver/etsi_ts/102400_102499/102427/01.01.01_60/ts_102427v010101p.pdf)、[ORコード](https://www.thonky.com/qr-code-tutorial/format-version-information)などでも使われている機能です。

- リードソロモン符号のエンコード: バイト列データを入力とし、パリティデータを付与する符号化バイト列を出力
- リードソロモン符号のデコード: (エラー混入した)符号化バイト列を入力とし、エラー補正したバイト列データを出力

リードソロモン符号では、有限体GFおよび、符号化バイト列長を`N`、データバイト列長を`K`、非負整数値bというパラメータがあります(bの意味は後述)。**エンコード側とデコード側双方で同じパラメータを共有**することが前提です。
また、パラメータの有限体では、同じ原始元を用いることが必須となります。つまり、PF(p)の場合は原始元となる整数値aを、GF($p^n$)の場合は多項式aを解とする既約多項式fを、pやnとともに指定する必要があります。

ただし、リードソロモン符号で用いる有限体GFは、多くのケースで、1バイト8ビットのビット列に対応するGF($2^8$)を用いています。

パリティデータのバイト長dは、N-Kとなります。パリティデータは素データの直後似付与されます。
リードソロモン符号で**訂正可能なエラー数は、dの1/2以下**までです。これは素データとパリティデータをあわせた、全符号化データ内でのエラー出現数であり、どこにエラーが出てもデコードでの扱いの差はありません。

# リードソロモン符号の実装

リードソロモン符号の内部では、**入出力のバイト列は、有限体GF係数の多項式**である、として扱います。

リードソロモン符号実装`rscode.js`では、中編で実装したコードを使用します。
多項式として`Polynomial`を用い、パラメータとして`GF`や`PF`、`GF2n`を受け付けます。また、デコード実装では、連立方程式を解くための`FieldUtils`を使用します。

<details>
<summary>`rscode.js`コード全体</summary>
<div>

```js:rscode.js
// Reed-Solomon code
import {range, FieldUtils} from "./field-utils.js";
import {Polynomial, PolynomialUtils} from "./polynomial.js";

export const RSCode = (FF, N, K, b, debug = false) => {
  const d = N - K, t = d >>> 1;
  const poly = Polynomial(FF), futils = FieldUtils(FF);
  // explict conversion between FF-element array and FF-polynomial
  const a2p = poly.fromArray, rev2p = a => a2p([...a].reverse());
  const p2a = poly.toArray, p2rev = (e, len) => [...p2a(e, len)].reverse();

  // gen = (x-a^b)*(x-a^(b+1)*...*(x-a^(b+d-1)) as extracted poly
  const rsGenerator = () => poly.prod(range(d).map(k => poly.add(
    poly.monomial(FF.one(), 1),
    poly.monomial(FF.neg(FF.alpha(b + k)), 0))));
  const gen = rsGenerator();
  if (debug) console.log("gen", gen);

  const encode = msg => {
    const base = poly.carry(rev2p(msg.map(FF.fromNum)), d);
    const coded = poly.sub(base, poly.mod(base, gen));
    if (debug) console.assert(poly.mod(coded, gen).every(FF.isZero));
    return p2rev(coded, N).map(FF.toNum);
  };

  const rsSyndromes = cx => range(d).map(
    k => poly.apply(cx, FF.alpha(b + k)));
  const hasError = synds => !synds.every(si => FF.isZero(si));

  const rsLx = (synds) => {
    const leqs0 = range(t).map(i => range(t + 1).map(j => synds[i + j]));
    const v = futils.rank(leqs0);
    const leqs = v === t ? leqs0 :
          range(v).map(i => range(v + 1).map(j => synds[i + j]));
    const rlx = futils.solve(leqs);
    return a2p([FF.one(), ...rlx.reverse()]);
  };
  const rsPositions = lx => range(N).filter(
    k => FF.isZero(poly.apply(lx, FF.alpha(-k))));

  const rsOx = (sx, lx) =>
        poly.mod(poly.mul(sx, lx), poly.carry(poly.one(), d));
  const rsDLx = lx => poly.diff(lx);
  const rsErrors = (positions, ox, dlx) => positions.map(k => {
    const akinv = FF.alpha(-k);
    const oAkinv = poly.apply(ox, akinv);
    const dlAkinv = poly.apply(dlx, akinv);
    const ak1b = FF.alpha(k * (1 - b));
    return FF.neg(FF.mul(ak1b, FF.div(oAkinv, dlAkinv)));
  });

  const rsEx = (positions, errors) =>
        poly.sum(positions.map((k, i) => poly.monomial(errors[i], k)));

  const decode = code => {
    const rx = rev2p(code.map(FF.fromNum));
    const synds = rsSyndromes(rx);
    if (debug) console.log("sx", synds);
    if (!hasError(synds)) return code.slice(0, K);

    const lx = rsLx(synds);
    //const lx = PolynomialUtils(poly).findLinearRecurrence(synds);
    if (debug) console.log("lx", lx);
    const positions = rsPositions(lx);
    if (debug) console.log("positions", positions);
    if (positions.length === 0) {
      throw Error(`Cannot recover: error count > ${t}`);
    }

    const sx = a2p(synds);
    const ox = rsOx(sx, lx);
    const dlx = rsDLx(lx);
    const errors = rsErrors(positions, ox, dlx);
    if (debug) console.log("errors", errors);
    const ex = rsEx(positions, errors);
    const cx = poly.sub(rx, ex);
    return p2rev(cx, N).slice(0, K).map(FF.toNum);
  };
  return {encode, decode};
};
```

</div>
</details>

<details>
<summary>`rscode.js`利用例と結果</summary>
<div>

QRコードでのリードソロモン符号の例として、

- https://kouyama.sci.u-toyama.ac.jp/main/education/2007/infomath/pdf/text/text09.pdf

の73ページの例を利用しました。

コード

```js:rscode-example.js
import {GF2n} from "./gf2n.js";
import {RSCode} from "./rscode.js";

{
  console.log("[Reed-Solomon error detective code]: N=26,K=7,b=0 on GF(2^8)");
  // GF(2^8) System for QR Code
  const n = 8, f = 0b100011101; //a^8+a^4+a^3+a^2+1=0
  // example from https://kouyama.sci.u-toyama.ac.jp/main/education/2007/infomath/pdf/text/text09.pdf
  const N = 26, K = 19, b = 0;
  const {encode, decode} = RSCode(GF2n(n, f), N, K, b);
  const msg = [
    0b10000000, 0b01000100, 0b10000101, 0b10100111, 0b01001001, 0b10100111,
    0b10001011, 0b01101100, 0b00000000, 0b11101100, 0b00010001, 0b11101100,
    0b00010001, 0b11101100, 0b00010001, 0b11101100, 0b00010001, 0b11101100,
    0b00010001,
  ];
  console.log("msg", msg);

  const coded = encode(msg);
  console.log("coded", coded);

  // error injection
  coded[24] ^= 0xff;
  //coded[16] ^= 0x02;
  coded[11] ^= 0xac;
  console.log("error coded", coded);

  const fixed = decode(coded);
  console.log("fixed", fixed);

  console.log("Same as original", fixed.every((c, i) => c === msg[i]));
}
```

結果

```
[Reed-Solomon error detective code]: N=26,K=7,b=0 on GF(2^8)
msg [
  128,  68, 133, 167, 73, 167,
  139, 108,   0, 236, 17, 236,
   17, 236,  17, 236, 17, 236,
   17
]
coded [
  128,  68, 133, 167,  73, 167, 139,
  108,   0, 236,  17, 236,  17, 236,
   17, 236,  17, 236,  17, 249, 187,
   11, 161,  75,  69, 244
]
error coded [
  128,  68, 133, 167,  73, 167, 139,
  108,   0, 236,  17,  64,  17, 236,
   17, 236,  17, 236,  17, 249, 187,
   11, 161,  75, 186, 244
]
fixed [
  128,  68, 133, 167, 73, 167,
  139, 108,   0, 236, 17, 236,
   17, 236,  17, 236, 17, 236,
   17
]
Same as original true
```

この`RSCode()`は、`GF2n(n, f)`だけでなく、PF(p)や`GF(p, n, f)`でも同様に機能します。`GF.toNum()`、`GF.fromNum()`により、`GF`であっても数値配列が入出力となります(GFの既約多項式fは配列表現です)。

</div>
</details>

以下、この`rscode.js`のコードの順をおって説明していきます。

## リードソロモン符号のパラメータ

```js
export const RSCode = (FF, N, K, b, debug = false) => {
  const d = N - K, t = d >>> 1;
  const poly = Polynomial(FF), futils = FieldUtils(FF);
```

`RSCode(FF, N, K, b)`のパラメータFFは多項式の係数にする有限体です。
この`FF`としては、先に実装した`GF2n(n, f)`、`PF(p)`、`GF(p, n, f)`が使えます。

リードソロモン符号では、入力データとパリティデータをあわせた符号データのバイト数を`N`、入力データのバイト数を`K`として用いています。そして、符号データバイト数Nは、有限体FFの要素数-1(つまり$p^n-1$)以下である必要があります。

パリティデータ配列の長さ`d`は、N-Kです。
そして、最大訂正可能エラー数`t`は、d/2を超えない整数値となります。

最後のパラメータ`b`は整数値で、0や1などを使います。bを変えるとパリティデータ値は変化します。このbの意味は後の説明の中で行います。

実装によっては、後述の生成多項式g(x)もパラメータ扱いになります。
ただし、この多項式値g(x)は、他のパラメータFF、N、K、bで算出できる値です。

### 入出力データとしての有限体係数の多項式

```js
  const a2p = poly.fromArray, rev2p = a => a2p([...a].reverse());
  const p2a = poly.toArray, p2rev = (e, len) => [...p2a(e, len)].reverse();
```

リードソロモン符号の入出力バイト列は、有限体GF($2^8$)を係数とする多項式として扱います。

ただし、これまのインデックスを指数値とする多項式の表現とは逆で、リードソロモン符号のデータでは、**バイト列の先頭バイトが多項式での最上位桁**になります。

この実装では、多項式計算での簡便性のためバイト列の配列を入出力時に`reverse()`することで、`Polynomial(FF)`として使う想定にしています。

この逆順配列aと多項式eの相互変換のために`rev2p(a)`と`p2rev(e)`を用意します。

また、コード上でデータ配列と有限体係数多項式を使い分けるために、ArrayとPolynomialの相互変換を示す`a2p`や`p2a`も用意しています。

## リードソロモン符号の生成多項式

```js
  const rsGenerator = () => poly.prod(range(d).map(k => poly.add(
    poly.monomial(FF.one(), 1),
    poly.monomial(FF.neg(FF.alpha(b + k)), 0))));
  const gen = rsGenerator();
```

リードソロモン符号でのエンコードでは、パリティデータ長`d`に応じた「**生成多項式**(generator polynomial)」という多項式を用意します。

リードソロモン符号での生成多項式g(x)は、有限体の**原始元aの指数表現で連続する値**$a^0$、$a^1$、$a^2$、$a^3$、...$a^{d-1}$`を解に持つ多項式です。つまり、

- $g(x) = \prod_{i=0}^{d-1}(x-a^{i})$

を展開した多項式 $x^d-(a^0+a^1+a^2+...+a^{d-1})x^{d-1}+...+a^{0+1+2+...+d-1}$ を`rsGenerator()`で作成しています。

前出の`RSCode(FF, N, K, b)`の最後のパラメータ`b`は、この生成多項式の始点を$a^0=1$以外から始めるためのものであり、生成多項式は

- $g(x) = \prod_{i=0}^{d-1}(x-a^{b+i})$

を展開したものになります。上記のコードはこの右辺に相当します。

このパラメータ`b`はエラー値算出のときに影響し、b=1のとき計算がすこし少なくなる効果がありますが、CD, MPEG2, QR codeなどの仕様ではb=0であるようです。

そして、リードソロモン符号での生成多項式g(x)が持つ意味とは、i=b, ..., b+d-1のどれかのときは、**$g(a^i)=0$となる多項式**である、ということです。

## リードソロモン符号のエンコード

```js
  const encode = msg => {
    const base = poly.carry(rev2p(msg.map(FF.fromNum)), d);
    const coded = poly.sub(base, poly.mod(base, gen));
    if (debug) console.assert(poly.mod(coded, gen).every(FF.isZero));
    return p2rev(coded, N).map(FF.toNum);
  };
```

リードソロモン符号でのエンコードは、

- 入力メッセージを、多項式として、パリティデータの分だけ桁上げする
- 桁上げした多項式を、生成多項式でモジュロをとり、その余りの多項式を出す
- 桁上げしたメッセージから、余りの多項式を**引く**

を行うだけです。

最後に生成多項式g(x)での余りの多項式を引いているため、エンコードした結果の符号化データ多項式c(x)は、生成多項式g(x)で割り切れる多項式となっています。

- c(x) = d(x) * g(x)

つまり、エンコード後の**符号化データc(x)はg(x)で割り切れる**ことから、代入するとg(x)が0になる値**$a^b$,...$a^{b+d-1}$を代入すると0になる**多項式となっています。

つまり、多項式にこれらを代入して**すべてが(有限体での)0になればエラーはない**、ということになります。

逆に、**エラーが一つでもあれば、0にならない値が出てくる**、ということになります。

## エンコードしたデータと混入したエラーの関係

上記のエンコードした結果の、生成多項式g(x)で割り切れる、(エラーなし)符号化データの多項式をc(x)とします。

c(x)にエラーが混入したエラー混入後データの多項式をr(x)とし、c(x)との差分であるエラーも多項式e(x)で表現します。

- r(x) = c(x) + e(x)

この和は、多項式としての各項の係数ごとの和です。

そして、e(x)の0でない各係数を「**エラー値**(error value)」とよび、各エラー値ごとのe(x)の指数を「**エラー位置**(error position)」と呼びます。
このエラー位置は、多項式としての指数値であり、バイト列でのインデックス値ではありません(逆順であるため)。

リードソロモン符号で訂正可能なエラー数はt個まであるため、

- $e(x) = e_0x^{i_0} + e_1x^{i_1} + ... + e_{t-1}x^{i_{t-1}}$

となります。たとえば混入したエラーの個数が3個なら

- $e(x) = e_0x^{i_0} + e_1x^{i_1} + e_{2}x^{i_{2}}$

となります。この$e_0$、$e_1$、$e_2$が各エラー値、$i_0$、$i_1$、$i_2$が各エラー位置になります。

エラーなしデータ多項式c(x)は生成多項式g(x)で割り切れるため、g(x)=0となるxを代入すれば、c(x)も0になります。
たとえば、$c(a^b)=0$です。

よって、このc(x)が0となる$a^b$等を、エラー混入データ多項式r(x)に代入すると、

- $r(a^b) = c(a^b) + e(a^b) = e(a^b) = e_0a^{i_0 b} + e_1a^{i_1 b} + e_2a^{i_2 b}$

となって、**エラー多項式e(x)だけに関する(有限体の)値**になる性質があります。

そして、**符号データ長Nが$p^n-1$以下である**必要があるのは、この代入結果の値で、エラー位置$i_k$ごとの$a^{i_k}$が区別できる必要があるためです。もし$i_k$として$p^n-1$より大きいものがあるばあい、$i_k$と$i_{k_2} = i_k - (p^n-1)$との区別ができず、エラー位置の特定が不可能になります。

## デコード

リードソロモン符号のデコードは少し複雑で、以下の手順に分けられます:

- シンドローム値の算出
- エラー位置の導出
- エラー値の導出
- エラー値の補正

## シンドローム値の算出

```js
  const rsSyndromes = cx => range(d).map(
    k => poly.apply(cx, FF.alpha(b + k)));
  const hasError = synds => !synds.every(si => FF.isZero(si));
```

リードソロモン符号の「**シンドローム**」とは、デコード入力r(x)にg(x)が0になる$a^b$、$a^{b+1}$、$a^{b+d-1}$それぞれを代入した、`d`個の有限体の値(の配列)のことです。

- $s_k = r(a^{b+k}) = e(a^{b+k}) = \sum_{j=0}^{v-1}e_ja^{(b+k)i_j}$

先に示したように、この**シンドロームの値がすべて(有限体での)0であればエラーはない**、ということになります。

具体的には、b=０でエラーが3個あるとき、$s_k$は以下のようになっています。

- $s_0 = e_0 + e_1 + e_2$ 
- $s_1 = e_0a^{i_0} + e_1a^{i_1} + e_2a^{i_2}$ 
- $s_2 = e_0a^{2i_0} + e_1a^{2i_1} + e_2a^{2i_2}$ 
- $s_3 = e_0a^{3i_0} + e_1a^{3i_1} + e_2a^{3i_2}$ 
- ...
- $s_{d-1} = e_0a^{(d-1)i_0} + e_1a^{(d-1)i_1} + e_2a^{(d-1)i_2}$ 

エラー値算出のときは、このシンドロームの配列を、d-1次多項式s(x)として利用します。

- $s(x) = s_0 + s_1x + ... + s_{d-1}x^{d-1}$

## エラー位置の算出

```js
  const rsLx = (synds) => {
    const leqs0 = range(t).map(i => range(t + 1).map(j => synds[i + j]));
    const v = futils.rank(leqs0);
    const leqs = v === t ? leqs0 :
          range(v).map(i => range(v + 1).map(j => synds[i + j]));
    const rlx = futils.solve(leqs);
    return a2p([FF.one(), ...rlx.reverse()]);
  };
  const rsPositions = lx => range(N).filter(
    k => FF.isZero(poly.apply(lx, FF.alpha(-k))));
```

まず、検出可能エラー個数`t`以下となる、**実際のエラーの個数**を`v`とします。

k=0 ... v-1でのエラー位置$i_k$を求めるための、**エラー位置多項式**l(x)として、以下の式について考えます。

- $l(x) = \prod_{k=0}^{v-1}(1-a^{i_k}x) $

この多項式は、$l(a^{-i_k})=0$となるv次多項式になります(代入値は、指数が負である、つまり逆数であることに注意)。

v=3のときは、

```math
\begin{align}
l(x) &= (1-a^{i_0}x)(1-a^{i_1}x)(1-a^{i_2}x) \\
  &= 1 - (a^{i_0}+a^{i_1}+a^{i_2})x + (a^{i_0+i_1}+a^{i_0+i_2}+a^{i_1+i_2})x^2-a^{i_0+i_1+i_2}x^3
\end{align}
```

です。

このl(x)が得られたなら、$a^{-0}$、$a^{-1}$、...、$a^{-N+1}$を、このl(x)へ総当りで代入し、結果が0になるかどうかを判定することで、エラー位置が見つけられます。これが`rsPositions(lx)`でやっていることです。

この展開した$l(x) = 1 + l_1x + l_2x^2 + ... + l_v x^v$の係数$l_k$を算出しl(x)を作る関数が、`rsLx(synds)`です。

たとえば、v=3のときは、l(x)の各係数は以下の値に相当します。

- $l_1 = -(a^{i_0}+a^{i_1}+a^{i_2})$
- $l_2 = a^{i_0+i_1}+a^{i_1+i_2}+a^{i_2+i_0}$
- $l_3 = -a^{i_0+i_1+i_2}$

この時点では、l(x)で求める具体的な$i_k$は不明なので、$a^{i_k}$を含んでいるシンドローム値を用いることで、各$l_k$を求めます。

この$l_k$の算出では、以下のシンドローム値を係数としたv個の方程式でなる連立方程式を解くことで行います(各項でかけるsとlのインデックスが逆順関係である点に注意)。

- $s_0 l_v + s_1 l_{v-1} + ... + s_{v-1} l_1 + s_v = 0$
- $s_1 l_v + s_2 l_{v-1} + ... + s_{v} l_1 + s_{v+1} = 0$
- ...
- $s_{v-1} l_v + s_v l_{v-1} + ... + s_{2v-2} l_1 + s_{2v-1} = 0$

導出時はまず、v=tと仮定して、この連立方程式を作ります。すると、この結果得られた連立方程式の係数行列のランクはvになります。
得られたランクをvとして、改めて連立方程式を作り、解くことでl(x)の係数が求まります。

v=3なら以下になります。

- $s_0 l_3 + s_1 l_2 + s_2 l_1 + s_3 = 0$
- $s_1 l_3 + s_2 l_2 + s_3 l_1 + s_4 = 0$
- $s_2 l_3 + s_3 l_2 + s_4 l_1 + s_5 = 0$



<details>
<summary>参考: この連立方程式の意味</summary>
<div>
> v=3の例で、一番上の等式について分析します。

> $e_0$の項だけにフォーカスすると、$s_0 = e_0a^{0i_0}$、$s_1 = e_0a^{1i_0}$、$s_2 = e_0a^{2i_0}$、$s_3 = e_0a^{3i_0}$です。

> そして、$e_0$のエラー位置$i_0$にフォーカスして、つぎに$e_0$の各$l_k$部について考えます。

> まず$l_2$に$s_1$中の$a^{i_0}$をかけると、項の一つが$l_3=a^{i_0+i_1+i_2}$になります。
> その余りの$a^{2i_0+i_1} + a^{2i_0+i_2}$は、$l_1$に$s_2$中の$a^{2i_0}$をかけたものの中に含まれます。
> 最後に余るのが、$a^{3i_0}$だけになり、つまり$s_3$を構成する$e_0a^{3i_0}$と同じものとなります。

> このため、$e_0$でくくると、その中は0になります。$e_1$、$e_2$でも同様にくるることで中は0になり、等式が成立します。

> 他の等式はs側で$a^{i_0}$をかける数が一つづつ増えただけの関係であり、同様にすべて0になります。

> そして、この関係は、一般のvについても同様に成立します。
</div>
</details>

<details>
<summary>参考: エラー個数とこの連立方程式のランクの関係</summary>
<div>
> 各方程式を$e_k$ごとにまとめると、b=0、t=3で実際のエラー数v=2なら

> - $e_0(l_3+a^{i_0}l_2+a^{2i_0}l_1) + e_1(l_3+a^{i_1}l_2+a^{2i_1}l_1) + s_3 = 0$
> - $e_0a^{i_0}(l_3+a^{i_0}l_2+a^{2i_0}l_1) + e_1a^{i_1}(l_3+a^{i_1}l_2+a^{2i_1}l_1) + s_3 = 0$
> - $e_0a^{2i_0}(l_3+a^{i_0}l_2+a^{2i_0}l_1) + e_1a^{2i_1}(l_3+a^{i_1}l_2+a^{2i_1}l_1) + s_3 = 0$

> となり、実質

> - $L_0 = l_3+a^{i_0}l_2+a^{2i_0}l_1$
> - $L_1 = l_3+a^{i_1}l_2+a^{2i_1}l_1$

> の2変数での連立方程式となるため、連立方程式のランクは2になります。

</div>
</details>

この連立方程式による計算方法は、"Peterson–Gorenstein–Zierler algorithm"(PGZアルゴリズム)といいます。

ちなみに、ハードウェア実装の場合は連立方程式を解くのに、[クラメルの公式](https://ja.wikipedia.org/wiki/%E3%82%AF%E3%83%A9%E3%83%A1%E3%83%AB%E3%81%AE%E5%85%AC%E5%BC%8F)で、[行列式](https://ja.wikipedia.org/wiki/%E8%A1%8C%E5%88%97%E5%BC%8F)を用いて解く実装が選ばれます。数あるPGZアルゴリズムの説明記事の中で、$l_k$をシンドローム値だけの式でいきなり表現しているものは、だいたいは行列式を展開したたものです。

<details>
<summary>参考: Berlekamp-Masseyアルゴリズムによるl(x)の算出</summary>
<div>

> PGZアルゴリズムでの$l_k$を求める連立方程式の方程式はすべて、

> - $s_k + s_{k-1}l_1 + s_{k-2}l_{2} + \dots + s_{k-v}l_v = 0$

> という１つの漸化式を満たすものです。
> つまり、シンドローム$s_k$全体は、この漸化式を満たす数列であるとみなせます。

> このため、中編にて`PolynomialUtils`で実装した、数列sからその線形漸化式を求める`findLinearRecurrence(s)`を用いることでも、l(x)を作ることが可能です。

> `findLinearRecurrence(s)`は、

> - $s_n + s_{n-1}c_1 - s_{n-2}c_{2} + \dots + s_{n-n}c_n = 0$

> が成立する多項式$c(x) = 1 + \sum_{k=1}^n c_k x^k$として返すため、`findLinearRecurrence(synds)`の結果は、l(x)そのものとなります。

</div>
</details>

注意する点は、もし実際の**エラー個数vがtより多い場合、デコード結果は不定**である点です。エラーを加えた結果が、**別の符号データからの訂正可能エラーとなりうる**場合も起こるからです。

ただし、どの符号データからも訂正可能でなければ、求めたl(x)は、どの$a^{-0}$, ..., $a^{-N+1}$を入れても0にならない式となり、この場合のみ、エラーが多すぎることが検知可能です。

つまり、

- 発生したエラーの個数が多いと、$エラー訂正できたとしても、もともとのデータに訂正されるとは限らない$

です。リードソロモン符号に限らない話ではありますが、うまく利用するには、用いる状況に応じた、検出可能エラー数を増やす、エラーの比率を下げる、などの設定が必要です。

## エラー値の算出


```js
  const rsOx = (sx, lx) =>
        poly.mod(poly.mul(sx, lx), poly.carry(poly.one(), d));
  const rsDLx = lx => poly.diff(lx);
  const rsErrors = (positions, ox, dlx) => positions.map(k => {
    const akinv = FF.alpha(-k);
    const oAkinv = poly.apply(ox, akinv);
    const dlAkinv = poly.apply(dlx, akinv);
    const ak1b = FF.alpha(k * (1 - b));
    return FF.neg(FF.mul(ak1b, FF.div(oAkinv, dlAkinv)));
  });
```

シンドローム多項式s(x)とエラー位置多項式l(x)の積のd未満の項でなす多項式o(x)、l(x)を微分した多項式dl(x)とします。

- $o(x) = s(x)l(x) \mod x^{d}$
- $dl(x) = \frac{d}{dx}l(x) = \sum_{k=0}^{v-1}{(-a^{i_k}\prod_{j\ne k}{(1-a^{i_j}x)})}$

これらの多項式を使った以下の式で、エラー位置$i_k$からエラー値$e_k$が求まります。

- $e_k = \frac{-a^{i_k(1-b)}o(a^{-i_k})}{dl(a^{-i_k})}$

パラメータb = 1のときは、oをdlで割って負にするだけになります。

- $e_k = \frac{-o(a^{-i_k})}{dl(a^{-i_k})}$

<details>
<summary>参考: エラー値の等式の意味</summary>
<div>

> 前提知識として、$x^4-1 = (x-1)(x^3+x^2+x+1)$より、$1+x+x^2+x^3 = \frac{x^4-1}{x-1}$となります。
> これは一般的に、$\sum_{i=0}^{n-1}x^i = 1+x+...+x^{n-1} = \frac{x^n - 1}{x - 1}$ の等式であり、有限体係数であっても成立する関係です。

> まず単純化のため、パラメータb=0で考えます。

> シンドローム$s_k = e_0a^{ki_0} + e_1a^{ki_1} + e_{v-1}a^{ki_v}$で、kは0からd-1までのd個あります。
そして、前述のシンドローム関数s(x)は、エラー値$e_k$ごとにまとめると、以下のようになります。

> - $s(x) = e_0(1+a^{i_0}x+(a^{i_0}x)^2 + \dots + (a^{i_0}x)^{d-1}) + e_1(1+a^{i_1}x+(a^{i_1}x)^2 + \dots + (a^{i_1}x)^{d-1}) + \dots +  e_{v-1}(1+a^{i_{v-1}}x + \dots + (a^{i_{v-1}}x)^{d-1})$

> これに、上述の等式を適用することで、

> - $s(x) = e_0\frac{(a^{i_0}x)^d-1}{a^{i_0}x-1} + e_1\frac{(a^{i_1}x)^d-1}{a^{i_1}x-1} + \dots + e_{v-1}\frac{(a^{i_{v-1}}x)^d-1}{a^{i_{v-1}}x-1}$

> となります。この結果の式に、$l(x) = \prod_{k=0}^{v-1}(1-a^{i_k}x)$をかけると、各項にある分母が消え-1になり、

> - $s(x)l(x) = -e_0 \left( ((a^{i_0}x)^d-1)\prod_{k \ne 0}(1-a^{i_k}x) \right) - e_1\left( ((a^{i_1}x)^d-1)\prod_{k \ne 1}(1-a^{i_k}x)\right) - \dots - e_{v-1}\left( ((a^{i_{v-1}}x)^d-1)\prod_{k \ne {v-1}}(1-a^{i_k}x)\right) $ 

> となります。

> ここで、$(a^{i_k}x)^d-1 \mod x^d = -1$となることを利用します。

> $o(x) = s(x)l(x) \mod x^d$より、

> - $o(x) = e_0\prod_{k \ne 0}(1-a^{i_k}x) + e_1\prod_{k \ne 1}(1-a^{i_k}x) + \dots + e_{v-1}\prod_{k \ne {v-1}}(1-a^{i_k}x) = \sum_{j=0}^{v-1}(e_j\prod{k\ne j}(1-a^{i_k}x))$

> となります(参考: 結果のo(x)は$v - 1$次多項式になるため、$\mod x^d$でなくても、$\mod x^v$でもよい)。

> このo(x)は、たとえばエラー位置$i_0$に対応する$a^{-i_0}$を代入すると、$(1-a^{i_0}x)$を含む$e_0$以外の項が0になります。

> - $o(a^{-i_0}) = e_0 \prod_{k \ne 0}(1-a^{i_k - i_0})$

> となります。一般のb != 0の場合は、

> - $o(a^{-i_0}) = e_0a^{bi_0} \prod_{k \ne 0}(1-a^{i_k - i_0})$

> です。

> そして、$a^{-i_0}$を代入すると、$(1-a^{i_0}x)$を含む項が0になるのは、l(x)の微分多項式dl(x)でも同様です。

> - $dl(a^{-i_0}) = -a^{i_0}\prod_{k \ne 0}(1-a^{i_k-i_0}) $

> よって、

> - $\frac{o(a^{-i_0})}{dl(a^{-i_0})} = -e_0a^{i_0(b-1)} $

> $i_0$、$e_0$を一般化すると、

> - $\frac{o(a^{-i_k})}{dl(a^{-i_k})} = -e_ka^{i_k(b-1)} $


> となります。

> この右辺の$e_k$以外の部分を打ち消すために、$-a^{i_k(1-b)}$をかけることで、$e_k$が得られることになります。

> - $-a^{i_k(1-b)}\frac{o(a^{-i_k})}{dl(a^{-i_k})} = e_k $

</div>
</details>

このアルゴリズムは "Forney algorithm"(もしくはForney's formula)といいます。

このForneyアルゴリズムを使わなくても、エラー位置からすべての$a^{i_k}$がわかっているので、エラー個数と同じv個のシンドロームの等式で、$e_k$を解くための連立方程式を作って、解くことでもエラー値を求めることができます。

## エラー値の多項式

```js
  const rsEx = (positions, errors) =>
        poly.sum(positions.map((k, i) => poly.monomial(errors[i], k)));
```

`rsEx(positions, erros)`は、エラー位置とそれに対応するエラー値のリストから、エラーの多項式e(x)を復元する関数です。

このe(x)を入力多項式r(x)から引くことで、補正された多項式c(x)が求まります。

## デコード関数本体

```js
  const decode = code => {
    const rx = rev2p(code.map(FF.fromNum));
    const synds = rsSyndromes(rx);
    if (!hasError(synds)) return code.slice(0, K);

    const lx = rsLx(synds);
    //const lx = PolynomialUtils(poly).findLinearRecurrence(synds);
    const positions = rsPositions(lx);
    if (positions.length === 0) {
      throw Error(`Cannot recover: error count > ${t}`);
    }

    const sx = a2p(synds);
    const ox = rsOx(sx, lx);
    const dlx = rsDLx(lx);
    const errors = rsErrors(positions, ox, dlx);
    const ex = rsEx(positions, errors);
    const cx = poly.sub(rx, ex);
    return p2rev(cx, N).slice(0, K).map(FF.toNum);
  };
```

`decode(code)`は、ここまで実装してきた

- シンドローム算出
- エラー位置算出
- エラー値算出
- エラー多項式で補正

を順に行い、その結果の補正済み多項式c(x)から、メッセージ部分だけを切り出して返す関数です。

## BCH符号

リードソロモン符号と似たエラー訂正符号として、[BCH符号](https://ja.wikipedia.org/wiki/BCH%E7%AC%A6%E5%8F%B7)があります。

### リードソロモン符号からみたBCH符号の仕組み

BCH符号も、生成多項式で割り切れるようにするエンコード、シンドロームを計算し、エラー位置とエラー値を求めるデコードといった、大枠はリードソロモン符号と一緒です。

しかし、BCH符号では、**デコードでのシンドロームやエラー位置の計算はGF($p^n$)係数の多項式**を用いるに対し、**入出力データは、その部分体GF(p)係数の多項式**として扱うものとなっています。
つまり、GF($2^n$)なBCH符号の場合では、入出力データはバイト列ではなく、ビット列(の多項式)を扱うものとなります。

- BCH符号の符号化データのビット数N: $p^n-1$までの整数
- BCH符号のエラー個数: t (BCH符号ではtを指定する)
- BCH符号の生成多項式: GF($p^n$)の原始元aの$a^1$,...,$a^{2t}$を解にもつ、**既約多項式(GF(p)係数の多項式)の最公倍多項式**
- BCH符号のエンコード: 生成多項式で割り切れるためのGF(p)係数多項式でパリティデータ列を付加
- BCH符号のシンドローム: **GF(p)係数値をGF($p^n$)係数値として解釈**した符号化データ多項式に、生成多項式を0にするそれぞれ値を代入した結果の**GF($p^n$)値**
- BCH符号のエラー位置多項式: PGZアルゴリズムの連立方程式で求めるG($p^n$)係数の多項式
- BCH符号のエラー多項式: GF($p^n$)のForneyアルゴリズムで部分体GF(p)値となる係数値を得るG(p)の多項式(ただし、GF(2)係数ならエラー値は1なので自明)

ポイントは、BCH符号の生成多項式g(x)が、GF($p^n$)の値を係数にする多項式ではなく、**GF($p^n$)の値を解に持つGF(p)を係数とする既約多項式**の積で構築する点です。
この結果、g(x)のモジュロ算を行うエンコードではGF(p)係数の多項式として扱う一方、デコードでは**GF(p)の値はGF($p^n$)の部分体GF(p)の値である**とし、入力データをG($p^n$)の多項式として扱うことで、リードソロモン符号と同様に行います。

ただし、求めるエラー位置はビット列のインデックスになり、エラー値はビットインデックスでのビット値になります。そして**GF(2)の多項式であれば、エラー値は1のみ**になり、Forneyアルゴリズム等を用いる必要がなくなります。

一方で、リードソロモン符号より係数値のバリエーションが少ないので、混入したエラー個数が指定のエラー個数を超えるとき、別の符号データの訂正可能データとなる可能性は高くなるでしょう。

### BCH符号の実装

BCH符号の共有パラメータは、素数べき有限体GF(2,n,f)、符号化データビット長N、エラー訂正数t、非負整数値cで構成します。パラメータcは、リードソロモン符号におけるbと同じものですがc=1が普通のようです。

与えられたエラー訂正数tから、生成多項式の次数、つまりパリティビット長Dが求まるため、素データビット数はNとtから算出されます。



<details>
<summary>`bchcode.js`コード全体</summary>
<div>

注: このコードは、有限体係数多項式での計算によるアルゴリズムの実装です。よりビット列に特化させた実装は可能です。

```js:bchcode.js
import {range, uniq, FieldUtils} from "./field-utils.js";
import {Polynomial, PolynomialUtils} from "./polynomial.js";
import {GFUtils} from "./gf.js";

export const BCHCode = (gf, N, t, c) => {
  const gfutils = GFUtils(gf), gfPoly = Polynomial(gf);
  const futils = FieldUtils(gf);
  const pfPoly = gf.poly;

  const pfPoly2gfPoly = (pfp, n) =>
        gfPoly.fromArray(pfPoly.toArray(pfp, n).map(gf.fromSF));
  const gfPoly2pfPoly = (gfp, n) =>
        pfPoly.fromArray(gfPoly.toArray(gfp, n).map(gf.toSF));

  // encode
  const bchGenerator = () => {
    const polys = uniq(range(2 * t).map(
      k => gfutils.findIrreducible(gf.alpha(c + k))), pfPoly.eql);
    return pfPoly.prod(polys);
  };
  const gen = bchGenerator();
  const d = pfPoly.order(gen), K = N - d;

  const encode = msg => {
    const m = pfPoly.carry(msg, d);
    const r = pfPoly.mod(m, gen);
    return pfPoly.sub(m, r);
  };

  // decode
  const bchSyndromes = rx =>
        range(2 * t).map(k => gfPoly.apply(rx, gf.alpha(c + k)));
  const hasError = synds => !synds.every(gf.isZero);

  const bchLx = synds => {
    const leqs0 = range(t).map(i => range(t + 1).map(j => synds[i + j]));
    const v = futils.rank(leqs0);
    const leqs = v === t ? leqs0 :
          range(v).map(i => range(v + 1).map(j => synds[i + j]));
    const rlx = futils.solve(leqs);
    return [gf.one(), ...rlx.reverse()];
  };
  const bchPositions = lx => range(N).filter(
    k => gf.isZero(gfPoly.apply(lx, gf.alpha(-k))));

  const bchOx = (sx, lx) =>
        gfPoly.mod(gfPoly.mul(sx, lx), gfPoly.carry(gfPoly.one(), t));
  const bchDLx = lx => gfPoly.diff(lx);

  const bchErrors = (positions, ox, dlx) => positions.map(k => {
    const akinv = gf.alpha(-k);
    const oAkinv = gfPoly.apply(ox, akinv);
    const dlAkinv = gfPoly.apply(dlx, akinv);
    const ak1b = gf.alpha(k * (1 - c));
    return gf.neg(gf.mul(ak1b, gf.div(oAkinv, dlAkinv)));
  });

  const bchErrorsPF2 = positions => gfPoly.sum(positions.map(
    k => gfPoly.carry(gfPoly.one(), k)));

  const bchEx = (positions, errors) =>
        gfPoly.sum(positions.map((k, i) => gfPoly.monomial(errors[i], k)));

  const decode = code => {
    const rx = pfPoly2gfPoly(code, N);
    const synds = bchSyndromes(rx);
    if (!hasError(synds)) return gfPoly2pfPoly(rx.slice(d), K);

    //const lx = bchLx(synds);
    const lx = PolynomialUtils(gfPoly).findLinearRecurrence(synds);
    const positions = bchPositions(lx);
    if (positions.length === 0) {
      throw Error(`Cannot recover: error count > ${t}`);
    }

    if (pfPoly.K.size() === 2) {
      const ex = bchErrorsPF2(positions);
      const cx = gfPoly.sub(rx, ex);
      return gfPoly2pfPoly(cx.slice(d), K);
    } else {
      const sx = synds;
      const ox = bchOx(sx, lx);
      const dlx = bchDLx(lx);
      const errors = bchErrors(positions, ox, dlx);

      const ex = bchEx(positions, errors);
      const cx = gfPoly.sub(rx, ex);
      return gfPoly2pfPoly(cx.slice(d), K);
    }
  };

  return {K, d, gen, encode, decode};
};
```

デコードでは、Berlekamp-Masseyアルゴリズムでlxを求める実装を用いるようにしていますが、PGZアルゴリズムのコード`bchLx()`も載せてあります。


</div>
</details>

<details>
<summary>`bchcode.js`利用コード例</summary>
<div>

```js:bchcode-example.js
{
  console.log("[BCHCode with GF2n]");
  const gf = GF2n(4, 0b10011);
  const N = 15, t = 2, c = 1;
  const bch = BCHCode(gf, N, t, c);

  const msg = 0b0000101;
  console.log("msg", msg.toString(2).padStart(bch.K, "0"));
  let coded = bch.encode(msg);
  console.log("coded", coded.toString(2).padStart(N, "0"));
  coded = coded ^ (1 << 14);
  coded = coded ^ (1 << 12);
  //coded = coded ^ (1 << 1);
  console.log("error coded", coded.toString(2).padStart(N, "0"));

  const fixed = bch.decode(coded);
  console.log("fixed", fixed.toString(2).padStart(bch.K, "0"));
  console.log();
}
```

実行結果

```
[BCHCode with GF2n]
msg 0000101
coded 000010100110111
error coded 101010100110111
fixed 0000101
```

`BCHCode'では、素数べき要素数の有限体として`GF2n'だけでなく、`GF'も利用可能ですが、入出力データはGF内の要素の表現形式と同じものになります。
つまり、`GF2n`であればビット列の`PF2Polynomial`であり、`GF`なら先頭を下位とする整数配列の`Polynomial(PF(p))`になります。

</div>
</details>

## 一般化されたBCH符号とリードソロモン符号

BCH符号では、GF(p)係数の入出力データと、その拡大体であるGF($p^n$)係数でのエラー補正計算を行うものでした。

前編では、GF($(p^n)^m$)を、GF($p^n$)係数の多項式を用いて構築しました。
BCH符号の仕組みでは、GF($p^n$)係数の入出力データと、その拡大体であるGF($(p^n)^m$)係数でのエラー補正計算を行うことも可能です。

この、素数要素数に限らない有限体GF(q)とその拡大体GF($q^m$)とで、エンコードとデコードを構成するよう一般化したものが、「一般化されたBCH符号」となります。
一般化されたBCH符号でのn=1、m>1のケースが、ビット列を扱う狭義のBCH符号になります。

そして、この「一般化されたBHC符号」の上で、m=1のケースが、リードソロモン符号相当になります。
GF(q)の値sを解とする一次多項式x-sは、GF(q)係数での既約多項式であり、その総乗は最小公倍多項式です。



# 参考: 巡回符号

リードソロモン符号やBCH符号は[巡回符号](https://ja.wikipedia.org/wiki/%E5%B7%A1%E5%9B%9E%E7%AC%A6%E5%8F%B7)の一種です。

巡回符号とは、符号データを、循環シフトさせたデータも必ず符号データになるタイプの符号システムのことです。

前編の最後で、GF($p^n$)では、$x^{p^n-1}-1 = \prod_{k=0}^{p^n-2}(x-a^k)$が成立する、ことを示しました。

これはいくつかの$(x-a^k)$の積であるリードソロモン符号やBCH符号の生成多項式g(x)はともに、**$x^{p^n-1}-1$を割り切る**ことを意味します。

- $x^{p^n-1}-1 = b(x)g(x)$

ある符号データ$c(x)=d(x)g(x)$とし、その最上位係数を$v$とすると、1つ巡回したデータは

- $c'(x) = xc(x)+v - vx^{p^n-1} = xc(x) - v(x^{p^n-1}-1)$

であり、c(x)、$x^{p^n-1}-1$それぞれのg(x)での表現の式を代入すれば、

- $c'(x) = xd(x)g(x)-vb(x)g(x) = (xd(x)-vb(x))g(x)$

となり、必ずg(x)で割り切れることが示されました。

一般化すると、**$x^{N}-1$を割り切るg(x)**で割り切れる、任意のN-1次データc(x)は、**巡回させたc'(x)も、必ずg(x)で割り切れるデータになる**、という性質を持ちます。

そして、符号化データがこのタイプのデータとなる符号システムのことを、[巡回符号](https://ja.wikipedia.org/wiki/%E5%B7%A1%E5%9B%9E%E7%AC%A6%E5%8F%B7)と呼びます。

チェックサムの[CRC](https://ja.wikipedia.org/wiki/%E5%B7%A1%E5%9B%9E%E5%86%97%E9%95%B7%E6%A4%9C%E6%9F%BB)も、巡回符号の一種となります。

巡回符号な符号化システムでは、桁上げ(シフト)と(多項式での)モジュロ算で、符号化(エンコード)が実装できる、とのことです。