前編内容

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

中編内容

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

後編内容

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

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

この記事内のコードで用いるJavaScriptの仕様は、Array.flatMap()のあるECMAScript2019仕様までを用います。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]);

リードソロモン符号とは

有限体の応用として、エラー検出だけでなくエラー訂正も可能なリードソロモン符号があります。

このリードソロモン符号はCDやMPEG2、ORコードなどでも使われている機能です。

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

リードソロモン符号では、有限体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を使用します。

`rscode.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};
};

`rscode.js`利用例と結果

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

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

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

コード

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は配列表現です)。

以下、このrscode.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で算出できる値です。

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

  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も用意しています。

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

  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$となる多項式である、ということです。

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

  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)$との区別ができず、エラー位置の特定が不可能になります。

デコード

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

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

シンドローム値の算出

  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}$

エラー位置の算出

  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のときは、

\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$

参考: この連立方程式の意味

> 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についても同様に成立します。

参考: エラー個数とこの連立方程式のランクの関係

> 各方程式を$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になります。

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

ちなみに、ハードウェア実装の場合は連立方程式を解くのに、クラメルの公式で、行列式を用いて解く実装が選ばれます。数あるPGZアルゴリズムの説明記事の中で、$l_k$をシンドローム値だけの式でいきなり表現しているものは、だいたいは行列式を展開したたものです。

参考: Berlekamp-Masseyアルゴリズムによるl(x)の算出

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)そのものとなります。

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

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

つまり、

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

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

エラー値の算出

  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})}$

参考: エラー値の等式の意味

前提知識として、$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 $

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

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

エラー値の多項式

  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)が求まります。

デコード関数本体

  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符号があります。

リードソロモン符号からみた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から算出されます。

`bchcode.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()も載せてあります。

`bchcode.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))`になります。

一般化された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符号は巡回符号の一種です。

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

前編の最後で、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)で割り切れるデータになる、という性質を持ちます。

そして、符号化データがこのタイプのデータとなる符号システムのことを、巡回符号と呼びます。

チェックサムのCRCも、巡回符号の一種となります。

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