LoginSignup
9
7

More than 5 years have passed since last update.

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

Posted at

勉強したので習作。誤り訂正符号に関する勉強の集大成。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. αを巡る計算ドリル

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

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

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

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

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

  • 性質3:$α^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符号と同じように考える。誤り訂正が3の場合、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符号と全く同じ。入力ビット列に対応する多項式表現を求めて、検査分(今回は8)だけ多項式の次数をあげて、生成多項式で剰余算をして、出てきた多項式に対応するビット列が検査ビットとなる。

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);
*/

参考:

9
7
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
9
7