勉強したので習作。「巡回ハミング符号を実装してみた」の知識を一部使う（2を法とした算術・シンドロームの概念等）

1. 概念理解

解きたい課題: 長さ7のビット列を転送したいのだが、転送中のノイズによりビット反転する可能性がある。そこで検査用に長さ8 のビット列を付け加え、全体で長さ 15のビット列を送信することにしよう。15 ビット中 2ビットのビット反転があったとしても、それを検知・訂正したい。

トレードオフを知る: 符号語（＝データ用ビット列＋検査用ビット列）の長さを固定した場合、検査用ビット列の長さを長くすればより多くの誤り訂正ができるが、その一方データビット列の長さが短くなるので転送効率が悪くなる。このように、誤り訂正と転送効率のトレードオフがある。

カスタマイズ性: BCH 符号はこのトレードオフに対して、柔軟なカスタマイズを提供。今回は検査用に長さ 8のビット列を採用し、2つの誤り訂正を実現する。検査用のビット列の長さを増減させ、誤り訂正機能の強弱を変えた符号も実装できる。

イメージ図:

送信したいビット列:           1101000
検査ビット列の付与:           110100010000001
-----送信--------------
ノイズによるビット反転:       110101010001001  (下から4, 10番目が反転)
-----受信--------------
検査:                       (誤りあり, 場所: 4, 10)
訂正後のビット列             110100010000001
検査ビットの剥離             1101000

2. 事前準備

生成多項式: $x^8 + x^7 + x^6 + x^4 + 1$ のベクトル表現を二進数として整数変換した 465 を使って、2を法とする算術を定義する:

// 生成多項式のベクトル表現を二進数として整数変換した値をgとして各種演算を定義
function mod2 (g) {
  const plus = (...args) => args.reduce((a, x) => a ^ x, 0);
  function remG (x) { // 入力: A(α)の整数表現, 出力: R(α)の整数表現 (R(x)は A(x)を生成多項式で割った余り)
    let result = x;
    while (true) {
      const shift = Math.clz32(g) - Math.clz32(result);
      if (shift < 0) break;
      result ^= g << shift;
    }
    return result;
  }
  function mul2 (x, y) {
    let result = 0;
    for (let pos = 0; pos < 15; pos++) {
      if (y & (1 << pos)) {
        result = plus(result, remG(x << (pos % 15)));
      }
    }
    return result;
  }
  return {
    remG,
    plus,
    mul: (...args) => args.reduce((a, x) => mul2(a, x), 1)
  };
}

const { remG, plus, mul } = mod2(465);

3. エンコード

巡回ハミング符号と全く同様。

検査ビット列の求め方：入力多項式に $x^8$ をかけたものを生成多項式で割った余り（に対応する整数）が検査ビットになる。よって、以下のようにエンコード関数が実装できる:

function encode (num) {
  const shifted = num << 8;
  return shifted + remG(shifted);
}

// 具体例
[0b101, 0b001, 0b10010].forEach(n => console.log(encode(n).toString(2).padStart(15, '0')));
/*
000010100110111
000000111010001
001001001001001
*/

4. デコード

大まかな流れは、「シンドローム計算」をしてから「誤り位置の計算」。誤り位置がわかれば復号は容易。

シンドローム計算: 4つのシンドローム値がある（数学的には 0,2番目から 1,3番目が導けるので実質的には2つ）。それぞれの値は、受信ビット列の多項式にα or α^2, α^3, α^4 を代入した数値であるので以下のように実装可能

// 拡大元α の α^n の値. α^15 = 1
const powAlpha = n => remG(1 << (n % 15));

// シンドローム計算: シンドロームは4つ。受信多項式にα,α^2, α^3, α^4を代入した値
function syndromeOf (num) {
  let syns = [0, 0, 0, 0];
  for (let pos = 0; pos < 15; pos++) {
    if (num & (1 << pos)) {
      syns = syns.map((s, i) => s ^ powAlpha((i + 1) * pos) );
    }
  }
  return syns;
}

誤り位置計算: 誤りがない場合には、4つのシンドローム値は全て0 になる。誤りが一つの場合には（巡回ハミング符号と同様）0番目のシンドローム値が、α^i (ただし i = 0..15) のどれかに一致し、その i こそが誤り位置になる。以下誤りが二つの場合について記述。

誤り多項式を $\sigma(x) = (x+\alpha^i)(x+\alpha^j)$ と定義しよう(i, j は誤り箇所）。すると $S_1\sigma(x) = S_1x^2 + S_1^2x + S_1^3 + S_3$ が容易に導ける（$S_1$は0番目のシンドローム値、$S_3$は2番目のシンドローム値）。この多項式に、$x = \alpha^0, \alpha^1, \alpha^2, ...\alpha^{14}$ を代入していくと、値が0になるのが二つだけ見つかるので、それが誤り位置だと気づく。コードに落とす:

function decode (num) {
  const syn = syndromeOf(num);
  // no error
  if (syn.every(x => x === 0)) {
    return num >> 8;
  }
  // one error
  const aPowArr = range(15).map(powAlpha);
  const oneError = aPowArr.findIndex(x => x === syn[0]);
  if (oneError !== -1) {
    const corrected = num ^ (1 << oneError);
    return corrected >> 8;
  }
  // two error
  const [S, , T] = syn;
  const sigma = x => plus(mul(S, x, x), mul(S, S, x), mul(S, S, S), T);
  const locations = aPowArr.map((x, i) => sigma(x) === 0 ? i : null).filter(x => x !== null);
  // console.log(locations);
  if (locations.length !== 2) {
    throw Error(`decode error: ${locations}`);
  }
  const corrected = locations.reduce((a, loc) => a ^ (1 << loc), num);
  return corrected >> 8;
}

5. 動作確認・テスト

入力は 7ビットなので、エンコーダは128種類の符号語を生成する。このそれぞれの符号語に対して、エラーなし(1パターン）、1ビット誤り（15パターン）、2ビット誤り（105パターン）が想定される。そのためテストするケースは、 128*(1+15+105) = 15488通りある。それぞれについて、符号化・復号が可能であることを確認する。また、符号語長は 15bit なので32768通りの符号語があることになる。引き算をして 17280通りの無意味な符号語をデコーダに渡して変なことが起きないこと（例外がスローされる）も確かめよう:

// 動作確認

// エンコードされる符号語. 全128種
const WORDS = range(1 << 7).map(x => [x, encode(x)]);

// 転送によって 1ビットの誤りを含む符号語. 全128*15種
const EWORDS = WORDS.reduce((a, w) => a.concat(range(15).map(loc => [w[0], w[1] ^ (1 << loc)])), [])

// 転送によって 1ビットの誤りを含む符号語. 全128*105種
const E2WORDS = WORDS.reduce((acc, word) => {
  let result = [];
  range(15).forEach(loc1 => {
    range(15).forEach(loc2 => {
      if (loc1 >= loc2) return;
      result.push([word[0], word[1] ^ (1 << loc1) ^ (1 << loc2)]);
    });
  });
  return acc.concat(result);
}, []);

const ALLWORDS = [...WORDS, ...EWORDS, ...E2WORDS];
console.log(ALLWORDS.every(([orig, word]) => orig === decode(word))); // -> true

const INVALIDWORDS = range(1 << 15).filter(x => !ALLWORDS.map(x => x[1]).includes(x));
console.log(1 << 15, ALLWORDS.length, INVALIDWORDS.length) // -> 32768 15488 17280

let [cnt1, cnt2] = [0, 0];
INVALIDWORDS.forEach(w => {
  try {
    const x = decode(w);
    cnt1 += 1;
  } catch (e) {
    cnt2 += 1;
  }
});
console.log(cnt1, cnt2); // -> 0 17280

6. 遊んでみた

6-1. 巡回符号であることの確認

巡回符号であることを確認してみよう。入力は7ビットなので 128種類の符号語にエンコードされるのだが、12 個のグループに分類でき、各グループの中では巡回操作によって重なり合うことがわかる:

const cycle15 = n => ((n << 1) + (1 & (n >> 14))) & ((1 << 15) - 1);

// 巡回操作で入力値と重なる数値の一覧
function groupOf (n) {
  const result = new Set();
  let val = n;
  range(15).forEach(i => {
    result.add(val);
    val = cycle15(val);
  });
  return [...result].sort((a, b) => a - b);
};

// 巡回符号であることを確認
const words = range(1 << 7).map(x => encode(x));
console.log(words.every(word => {
  const group = groupOf(word); // 自分を含むグループを一覧化し
  return group.every(elem => words.includes(elem)); // それらが実際に符号化されているか
}));

// いくつのグループがあるか -> 12
const names = new Set(words.map(x => groupOf(x)[0])); // ソートしているので先頭を代表とする
console.log(names.size);

names.forEach(n => console.log(n.toString(2).padStart(15, '0')));
/*
000000000000000
000000111010001
000001001110011
000010100110111
000011010010101
000101110111111
000110011111011
000111101011001
001001001001001
001011010101111
011011011011011
111111111111111
*/

コード全体

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

// 生成多項式のベクトル表現を二進数として整数変換した値をgとして各種演算を定義
function mod2 (g) {
  const plus = (...args) => args.reduce((a, x) => a ^ x, 0);
  function remG (x) { // 入力: A(α)の整数表現, 出力: R(α)の整数表現 (R(x)は A(x)を生成多項式で割った余り)
    let result = x;
    while (true) {
      const shift = Math.clz32(g) - Math.clz32(result);
      if (shift < 0) break;
      result ^= g << shift;
    }
    return result;
  }
  function mul2 (x, y) {
    let result = 0;
    for (let pos = 0; pos < 15; pos++) {
      if (y & (1 << pos)) {
        result = plus(result, remG(x << (pos % 15)));
      }
    }
    return result;
  }
  return {
    remG,
    plus,
    mul: (...args) => args.reduce((a, x) => mul2(a, x), 1)
  };
}

const { remG, plus, mul } = mod2(465);

function encode (num) {
  const shifted = num << 8;
  return shifted + remG(shifted);
}


// 拡大元α の α^n の値. α^15 = 1
const powAlpha = n => remG(1 << (n % 15));

// シンドローム計算: シンドロームは4つ。受信多項式にα,α^2, α^3, α^4を代入した値
function syndromeOf (num) {
  let syns = [0, 0, 0, 0];
  for (let pos = 0; pos < 15; pos++) {
    if (num & (1 << pos)) {
      syns = syns.map((s, i) => s ^ powAlpha((i + 1) * pos) );
    }
  }
  return syns;
}

function decode (num) {
  const syn = syndromeOf(num);
  // no error
  if (syn.every(x => x === 0)) {
    return num >> 8;
  }
  // one error
  const aPowArr = range(15).map(powAlpha);
  const oneError = aPowArr.findIndex(x => x === syn[0]);
  if (oneError !== -1) {
    const corrected = num ^ (1 << oneError);
    return corrected >> 8;
  }
  // two error
  const [S, , T] = syn;
  const sigma = x => plus(mul(S, x, x), mul(S, S, x), mul(S, S, S), T);
  const locations = aPowArr.map((x, i) => sigma(x) === 0 ? i : null).filter(x => x !== null);
  // console.log(locations);
  if (locations.length !== 2) {
    throw Error(`decode error: ${locations}`);
  }
  const corrected = locations.reduce((a, loc) => a ^ (1 << loc), num);
  return corrected >> 8;
}

// 動作確認

// エンコードされる符号語. 全128種
const WORDS = range(1 << 7).map(x => [x, encode(x)]);

// 転送によって 1ビットの誤りを含む符号語. 全128*15種
const EWORDS = WORDS.reduce((a, w) => a.concat(range(15).map(loc => [w[0], w[1] ^ (1 << loc)])), [])

// 転送によって 1ビットの誤りを含む符号語. 全128*105種
const E2WORDS = WORDS.reduce((acc, word) => {
  let result = [];
  range(15).forEach(loc1 => {
    range(15).forEach(loc2 => {
      if (loc1 >= loc2) return;
      result.push([word[0], word[1] ^ (1 << loc1) ^ (1 << loc2)]);
    });
  });
  return acc.concat(result);
}, []);

const ALLWORDS = [...WORDS, ...EWORDS, ...E2WORDS];
console.log(ALLWORDS.every(([orig, word]) => orig === decode(word)));

const INVALIDWORDS = range(1 << 15).filter(x => !ALLWORDS.map(x => x[1]).includes(x));
console.log(1 << 15, ALLWORDS.length, INVALIDWORDS.length)

let [cnt1, cnt2] = [0, 0];
INVALIDWORDS.forEach(w => {
  try {
    const x = decode(w);
    cnt1 += 1;
  } catch (e) {
    cnt2 += 1;
  }
});
console.log(cnt1, cnt2);

const cycle15 = n => ((n << 1) + (1 & (n >> 14))) & ((1 << 15) - 1);

// 巡回操作で入力値と重なる数値の一覧
function groupOf (n) {
  const result = new Set();
  let val = n;
  range(15).forEach(i => {
    result.add(val);
    val = cycle15(val);
  });
  return [...result].sort((a, b) => a - b);
};

// 巡回符号であることを確認
const words = range(1 << 7).map(x => encode(x));
console.log(words.every(word => {
  const group = groupOf(word); // 自分を含むグループを一覧化し
  return group.every(elem => words.includes(elem)); // それらが実際に符号化されているか
}));

// いくつのグループがあるか -> 12
const names = new Set(words.map(x => groupOf(x)[0])); // ソートしているので先頭を代表とする
console.log(names.size);

names.forEach(n => console.log(n.toString(2).padStart(15, '0')));

参考

以下を特に参照した