2元BCH符号を実装してみた（JavaScript）

勉強したので習作。「[巡回ハミング符号を実装してみた](https://qiita.com/41semicolon/items/083924c62276bcdc2e07)」の知識を一部使う（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を法とする算術を定義する:

```js
// 生成多項式のベクトル表現を二進数として整数変換した値を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$ をかけたものを生成多項式で割った余り（に対応する整数）が検査ビットになる。よって、以下のようにエンコード関数が実装できる:

```js
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 を代入した数値であるので以下のように実装可能

```js
// 拡大元α の α^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になるのが二つだけ見つかるので、それが誤り位置だと気づく。コードに落とす:

```js
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通りの無意味な符号語をデコーダに渡して変なことが起きないこと（例外がスローされる）も確かめよう:

```js
// 動作確認

// エンコードされる符号語. 全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 個のグループに分類でき、各グループの中では巡回操作によって重なり合うことがわかる:

```js
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
*/
```

# コード全体

```js
// 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')));
```

# 参考
以下を特に参照した

- [Wikipedia](https://ja.wikipedia.org/wiki/BCH%E7%AC%A6%E5%8F%B7)
- [富山大幸山研](https://kouyama.sci.u-toyama.ac.jp/main/computer/index.htm)が公開している[PDF](https://kouyama.sci.u-toyama.ac.jp/main/education/2012/infomath/pdf/text/codeBCH.pdf)
- [京大の講義資料](http://www.geocities.jp/kashi_pong/course_information_theory_2014.html)