前編内容
- 素数べき要素数の有限体のための基礎知識
- 素数べき要素数の有限体とその作り方
- 素数べき要素数の有限体の例
- 素数べき要素数の有限体のガロア群
中編内容
- JavaScriptで実装する素数べき要素数有限体
- 2べき要素数の有限体に特化した実装
後編内容
- 有限体の応用: リードソロモン符号とBCH符号
注: 記事内のJavaScriptコードについて
この記事内のコードで用いるJavaScriptの仕様は、Array.flatMap()のあるECMAScript2019仕様までを用います。JSファイルは、すべて拡張子は.jsのES moduleです。
非ECMAScriptな要素として、console.logとconsole.assertのみを用います。その他のWeb APIやnode.jsビルトイン等は一切使用していません。
この記事内のコードは、以下のgistに置いています。
多項式や有限体等の要素は、数値と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`コード全体
// 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コードでのリードソロモン符号の例として、
の73ページの例を利用しました。
コード
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$
参考: この連立方程式の意味
$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_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`コード全体
注: このコードは、有限体係数多項式での計算によるアルゴリズムの実装です。よりビット列に特化させた実装は可能です。
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`利用コード例
{
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も、巡回符号の一種となります。
巡回符号な符号化システムでは、桁上げ(シフト)と(多項式での)モジュロ算で、符号化(エンコード)が実装できる、とのことです。