リードソロモン符号(= RS符号) を実装してみた（JavaScript）

勉強したので習作。誤り訂正符号に関する勉強の集大成。QRコードで利用される。エンコードまではできたけど復号はできていない。

1. 全体像

解きたい課題: ノイズにより誤る可能性があるので 検査用に長さ7 のバイト列を付け加えることで、最大で247バイトのデータ（つまり、符号長の最大長は 254バイト）を送信する。符号の内 3つのバイトが壊れてもそれを検知・訂正する（バイトが壊れるとは、バイトを構成する8つのビットのうち、どれか 1つ以上のビットが反転する状態のこと）。

本稿での実践: 参考資料p73 の再現。19バイトの情報に 7バイトをくっつけて 26バイトの符号を作る

送りたい情報19バイト:
0x80, 0x44, 0x85, 0xA7, 0x49, 0xA7, 0x8B, 0x6C, 0x00, 0xEC,
0x11, 0xEC, 0x11, 0xEC, 0x11, 0xEC, 0x11, 0xEC, 0x11,
検査用の符号をくっつけて26バイトにする:
0x80, 0x44, 0x85, 0xA7, 0x49, 0xA7, 0x8B, 0x6C, 0x00, 0xEC,
0x11, 0xEC, 0x11, 0xEC, 0x11, 0xEC, 0x11, 0xEC, 0x11, 0xF9,
0xBB, 0x0B, 0xA1, 0x4B, 0x45, 0xF4

1. αを巡る計算ドリル

αの正体を知ることよりも重要なのは、 αの性質を知ること:

• 性質１：$α^0, α^1, α^2, ..., α^{254}$ の255種類の数は全部違う（同じものがない）
• 性質２：$α^{255} = 1$ が成立するので、後は何乗しても種類が増えない

僕らは、1バイト＝256種類のデータを何らかの形で数学の世界に持ち込みたい。そのまま自然数、つまり 1~256 に対応付けるというのも一つの方法なんだけど、 $\lbrace {α^i | i = 0, 1, ...254 } \rbrace$ に 0 を付け加えた集合に対応付けたほうが色々やりやすい（と数学者がいってるのでそれを信じる）。

α が住んでいる世界には、このようにαを何回もかけてできる 255種類の住人 と 0 という住人がいて、それで全員である。αのべき乗じゃないやつ、例えば、$1+α$ は住人じゃないの？というかもしれないけれど、実はそれも住人で、計算すればわかるんだけど $1+α = α^{25}$ が成立する。こんな感じで、任意のαの多項式が住人であるというすごく閉じた世界を作っている。

というか、そうなるように αが住む世界での足し算や掛け算を定義した、という方が正しくて、 それを2を法とした計算と呼んでいる。例えば、１＋１＝０といった感じ。このような計算ルールに加えて、αには、もう一つ重要な性質がある:

• 性質３：$α^8 + α^{4} + α^3 + α^2 + 1 = 0$

この式を再帰的に使うと αの次数をどんどん小さくできて、任意の次数のαの多項式を$α^7$ 以下の次数にまで変形できる。これを使うと、さっき出てきた $1+α = α^{25}$ が導ける、というわけ。

上記のような小さな世界 $GF(256)$ でよく使う関数を定義した（コードは末尾）:

• ELEMENTS: 長さ256の配列で、i 番目には $α^i$ の整数表現（$α^i$ を αの7次以下の多項式に変換したとき、各係数をビットとみなした8ビットを整数に変換した値）を保持。最後の255番目には 0 を格納。
• getIndexOf: $α$の多項式（の整数表現）が与えられたときに、それが、$α$ の何乗であるか？を調べてくれる関数。
• aPlus: 任意の整数 i, j に対して、$α^i + α^j = α^k$ が成立する。k の値を計算する関数。

2. 数学

バイト列と多項式表現: 入力を 19バイトとしよう。これは「19個の0..255の整数」を係数に持つ18次の多項式に対応付けることも可能ではある。例えば冒頭の 19バイトの入力は $128x^{18}+68x^{17}+...+236x+17$に対応付けることもできる。がしかし、それをやめよう。そうではなく、$\alpha^{7}x^{18}+\alpha^{102}x^{17}+...+\alpha^{122}x+\alpha^{100}$に対応付ける。このように入力ビット列を GF(256)の係数を持つ x の多項式に「翻訳」する。翻訳さえしてしまえば、数学の力を存分に発揮できる。

翻訳例.txt

<例: x^17の係数（多項式の前から二つ目の項の係数）を求める>
二つ目のバイト: 0x44 = 0b01000100
αの多項式表現と思い込むと: a^6 + a^2
これに対応する a^k は？: k = 102 (a^102 = a^6 + a^2)
x^17 の係数は: a^102

体を作りたい: 上記の方法で、任意の入力ビット列を変換して、「αの多項式を係数にもつ xの多項式」を得られた。やるべきは これを要素とするような体を作ること。

和の定義: 和とは、$A(x) + B(x) = C(x)$ となる C を求めることだが、x の次数が同じ者同士をまとめるという意味で「普通に」定義すると、結局は同じ次数同士の係数の足し算: $f(α)x^4 + g(α)x^4 = h(α)x^4$ となるhをどう定義しようか？ということに帰着する。この係数の足し算、つまりαの多項式同士の足し算として、GF(256)の演算を採用する。

積の定義: 積とは、$A(x)B(x) = C(x)$ となるCを求めることだけど、分配則を使って展開して、同じ係数でまとめてというように「普通に」定義すると、結局 $f(α)g(α) = h(α)$となる h をどう定義しようか？ということに帰着する。αの多項式同士の掛け算をGF(256)での掛け算として定義する。

剰余算: ハミング符号やBCH符号と同じように考える。誤り訂正が３の場合、7次の生成多項式 $G(x) = x^7 + α^{87} x^6 + α^{229} x^5 + α^{146} x^4 + α^{149} x^3 + α^{238} x^2 + α^{102} x + α^{21}$ となる（係数の計算は、前述のaPlus関数で求められる）。多項式に対して、G で割った剰余多項式を法とする（同じ剰余多項式を持つ多項式同士を同一視する）。

getIndexOf(0b1111111); // -> 87: x^6 の係数
getIndexOf(0b110011100110); // -> 229: x^5 の係数

3. エンコード

ハミング符号・BCH符号と全く同じ。入力ビット列に対応する多項式表現を求めて、検査分（今回は８）だけ多項式の次数をあげて、生成多項式で剰余算をして、出てきた多項式に対応するビット列が検査ビットとなる。

const G = [0, 87, 229, 146, 149, 238, 102, 21];

function encode (arr) {
  const aArr = arr.map(getIndexOf);
  const shifted = [...aArr, ...Array(G.length - 1).fill(255)];
  return [...aArr, ...rem256(shifted, G)].map(x => ELEMENTS[x]);
}

const data = [
  0x80, 0x44, 0x85, 0xA7, 0x49, 0xA7, 0x8B, 0x6C, 0x00, 0xEC,
  0x11, 0xEC, 0x11, 0xEC, 0x11, 0xEC, 0x11, 0xEC, 0x11,
];

console.log(encode(data).map(x => '0x' +  x.toString(16).toUpperCase()));
/*
[ '0x80', '0x44', '0x85', '0xA7', '0x49', '0xA7', '0x8B', '0x6C', '0x0', '0xEC',
  '0x11', '0xEC', '0x11', '0xEC', '0x11', '0xEC', '0x11', '0xEC', '0x11', '0xF9',
  '0xBB', '0xB', '0xA1', '0x4B', '0x45', '0xF4' ]
*/

とりあえず参考資料を再現できた。

4. デコード

ざっくり言って以下の手順

• シンドロームの計算
• 誤り数の判定
• 誤り箇所の判定
• 誤りの訂正

力尽きた。気が向いたらそのうちやろう。

• エラーなしの場合に Y(β) = 0 となることは確認
• エラーありの場合に Y(β) ≠ 0 となることは確認

なぜ力尽きたか

• 参考資料はエンコードまでしかやってなく、私に読める他の資料をさがしたが、見当たらなかった
• 本を買った方がよさそう
• 難しそう

5. コード

const assert = require('assert');

// utility
const range = end => [...Array(end).keys()];

// 2を法とした算術において 多項式f を多項式g で割った余り.入出力は対応する整数値
function rem2 (f, g) {
  let result = f;
  while (true) {
    const shift = Math.clz32(g) - Math.clz32(result);
    if (shift < 0) break;
    result ^= g << shift;
  }
  return result;
}

// α^n (n= 0, 1, ... n-1) を求める
function alphaPowersOf (n, g) {
  const result = [];
  let val = 1;
  range(n).forEach(() => {
    val = rem2(val, g);
    result.push(val);
    val = val << 1;
  });
  return result;
}

const F = 0b100011101; // αの原始多項式 z8+z4+z3+z2+1 に対応

// GF(255)の要素。順番にこそ意味があって i 番目は α^i に対応
// 最後、つまり 255番目に 0を加えて 256種類 = 1Byteで表現できる多項式全体をカバー
const ELEMENTS = [...alphaPowersOf(255, F), 0];

// 与えられたαの多項式が元の何番目か？（＝αの何乗か？）
// 例えば a12 + a7の場合 (1<<12) & (1<<7) を入力する
function getIndexOf(n) {
  return ELEMENTS.findIndex(x => x === rem2(n, F));
}

// α^i + α^j = α^k
function aPlus2 (i, j) {
  assert(i < 256 || j < 256);
  if(i === j) return 255; // 255番目は 0
  if(i === 255) return j;
  if(j === 255) return i;
  return getIndexOf(ELEMENTS[i] ^ ELEMENTS[j]);
}

const aPlus = (...args) => args.reduce(aPlus2, 255);

// 256を法とした算術において 多項式f を多項式g で割った余り.入出力は対応する配列
// 注！: 簡単のため g の最高次数の係数が 1 (gArr[0] === 0)であることを前提にしている(G は実際そう)
// 注！: 簡単のため f の各係数は 0..255 であることを前提（実際そう）
// rem2 関数と比較すると実装を理解しやすい
function rem256(fArr, gArr) {
  assert(gArr[0] === 0);
  assert(fArr.every(n => n < 256));
  let result = [...fArr]; // copy
  while (true) {
    const shift = result.length - gArr.length;
    if (shift < 0) break;
    // 商を α^q とする. G の最高次数の係数が1を仮定しているので簡単
    const q = result[0];
    // console.log(q, JSON.stringify(result.map(n => ELEMENTS[n])));
    range(gArr.length).forEach(i => {
      // 割る側の係数が 0 なら何もしない
      // 割る側がx3 + x + 1 のように穴あきの場合に考える必要
      if (gArr[i] === 255) {
        assert(false); // 我々の利用ケースでは来ないはず
        return;
      }
      result[i] = aPlus(result[i], (q + gArr[i]) % 255);
    });
    // 先頭の255(係数0に等しい)を取る [255, 255, 1,2,3] -> [1,2,3]
    assert(result[0] === 255);
    const where = result.findIndex(x => x !== 255);
    if (where === -1) {
      result = []; // 余りなし
    } else {
      result = result.slice(where);
    }
  }
  return result;
}

//const G = [0, 43, 139, 206, 78, 43, 239, 123, 206, 214, 147, 99, 150, 39, 243, 163, 136];
//const G = [0, 75, 249, 78, 6]; // 生成多項式 G: x4+ a75 x3 + a249 x2 + a75 x + a6
const G = [0, 87, 229, 146, 149, 238, 102, 21];

function encode (arr) {
  const aArr = arr.map(getIndexOf);
  const shifted = [...aArr, ...Array(G.length - 1).fill(255)];
  return [...aArr, ...rem256(shifted, G)].map(x => ELEMENTS[x]);
}

const data = [
  0x80, 0x44, 0x85, 0xA7, 0x49, 0xA7, 0x8B, 0x6C, 0x00, 0xEC,
  0x11, 0xEC, 0x11, 0xEC, 0x11, 0xEC, 0x11, 0xEC, 0x11,
];

/*
const word = encode(data);
const eword = [...word];
eword[3] = 0x10;

const A = ewords.map(getIndexOf);
const s = rem256(A, G);


let S = 0;

console.log(A, r.length, r);
console.log(E, s.length, s);
*/

参考: