本記事で説明する内容:

「2 元原始狭義 BCH 符号」の例
符号化
総当たりによる誤り訂正

本記事で説明しない内容:

一般の BCH 符号
その他の方法による誤り訂正

2 元 BCH 符号自体の説明は別記事にしました。

参考「2 元 BCH 符号入門: 符号化 - Qiita」

1. 有限体上の計算の準備

1.1. 指数表および対数表を作成する

有限体上の掛け算や割り算を効率良くするために、$\alpha^i$ の値を対応させる指数表および対数表を作成します。

ここでは、原始多項式を $x^4 + x + 1$ とする有限体 $GF(2^4)$ を使用することにします。

const GF_LOG_NEGATIVE_INFINITY = 15;

// 指数表を作成
const generateExponents = () => {

	const exponents = new Uint8Array(16);

	exponents[0] = 1;

	for (let log = 1; log < 15; log++) {
		const antilog = exponents[log - 1] << 1;
		exponents[log] = (antilog >>> 4 !== 0 ? antilog ^ 0b10011 : antilog);
	}

	exponents[GF_LOG_NEGATIVE_INFINITY] = 0;

	return exponents;

};

const GF_EXP = generateExponents();

// 対数表を作成
const generateLogarithms = exponents => {

	const logarithms = new Uint8Array(16);

	for (let i = 0; i <= 15; i++) {
		logarithms[exponents[i]] = i;
	}

	return logarithms;

};

const GF_LOG = generateLogarithms(GF_EXP);

// TODO: 
console.log(Array.from(
	GF_EXP,
	(antilog, i) => {
		const log = (i !== GF_LOG_NEGATIVE_INFINITY ? String(i).padEnd(3) : '-∞');
		const antilogBinary = antilog.toString(2).padStart(4, '0');
		return `α^${log} = 0b${antilogBinary} = ${antilog}`;
	},
).join('\n'));

実行結果

α^0   = 0b0001 = 1
α^1   = 0b0010 = 2
α^2   = 0b0100 = 4
α^3   = 0b1000 = 8
α^4   = 0b0011 = 3
α^5   = 0b0110 = 6
α^6   = 0b1100 = 12
α^7   = 0b1011 = 11
α^8   = 0b0101 = 5
α^9   = 0b1010 = 10
α^10  = 0b0111 = 7
α^11  = 0b1110 = 14
α^12  = 0b1111 = 15
α^13  = 0b1101 = 13
α^14  = 0b1001 = 9
α^-∞ = 0b0000 = 0

JavaScript では負の無限を -Infinity で表すことができますが Uint8Array 等では扱えないため、ここでは他の数を負の無限として扱うように定数 (GF_LOG_NEGATIVE_INFINITY = 15) を定義します。

ここで、GF_EXP[i] は $x = \alpha^i$ を表し、GF_LOG[x] は $i = \log_\alpha x$ を表します。ただし GF_EXP[GF_LOG_NEGATIVE_INFINITY] === 0 および GF_LOG[0] === GF_LOG_NEGATIVE_INFINITY です。

ここでは

\begin{align}
\alpha^0 &= 1 \\
\alpha^i &= (\alpha^{i-1} \cdot \alpha) \bmod (\alpha^4 + \alpha + 1)
\end{align}

のようにして計算を行います。

原始多項式 $x^4 + x + 1$ の係数を 2 進数で表すと 10011 になります。

1.2. 有限体上の四則演算

ここで使用する有限体の性質より、以下のようにして有限体上の四則演算を行えます:

足し算 $x + y$
- x ^ y
引き算 $x - y$
- x ^ y
掛け算 $x \times y$
- $x \neq 0 \land y \neq 0$ のとき GF_EXP[(GF_LOG[x] + GF_LOG[y]) % 15]
- $x = 0 \lor y = 0$ のとき 0
割り算 $x \div y$
- $x \neq 0 \land y \neq 0$ のとき GF_EXP[(15 + GF_LOG[x] - GF_LOG[y]) % 15]
- $x = 0$ のとき 0

2. 最小多項式を計算する

生成多項式を計算する準備として最小多項式の一覧を作成します。

// 最小多項式の一覧を作成
/**
 * 共役数を計算する
 */
// 意味的には
// {α^i, α^{2i}, α^{4i}, α^{8i}, ...}
const generateConjugations = i => {

	const conjugationsSet = new Set();

	for (let j = 0; j < 4; j++) {
		conjugationsSet.add((i << j) % 15);
	}

	const conjugations = new Uint8Array(conjugationsSet);

	return conjugations;

};

/**
 * 最小多項式を計算する
 */
// 意味的には
// m_i(x) = (x - α^i) (x - α^{2i}) (x - α^{4i}) (x - α^{8i}) ...
const generateMinimumPolynomial = conjugations => {

	const minimumPolynomialTypedArray = new Uint8Array(conjugations.length + 1);

	// 意味的には m_{i,-1}(x) = 1
	minimumPolynomialTypedArray[0] = 1;

	for (let j = 0; j < conjugations.length; j++) {

		// 意味的には
		// m_{i,j}(x) = m_{i,j-1}(x) * (x + α^{conjugations[j]})
		//            = m_{i,j-1}(x) * x + m_{i,j-1}(x) * α^{conjugations[j]}
		for (let k = j + 1; k > 0; k--) {

			if ( minimumPolynomialTypedArray[k - 1] === 0 ) continue;

			minimumPolynomialTypedArray[k] ^= GF_EXP[(
				GF_LOG[minimumPolynomialTypedArray[k - 1]] + conjugations[j]
			) % 15];

		}

	}

	// 意味的には
	// m_i(x) = m_{i,conjugations.length-1}(x)

	// 係数が 0 または 1 か確認
	if ( minimumPolynomialTypedArray.some(coefficient => 0 !== coefficient && 1 !== coefficient) ) {
		throw new Error('Invalid minimum polynomial');
	}

	// 配列を 2 進数として単一の数に変換
	const minimumPolynomial = minimumPolynomialTypedArray
		.reduce((binary, coefficient) => binary << 1 ^ coefficient);

	return minimumPolynomial;

};

const generateMinimumPolynomials = () => {

	const minimumPolynomials = new Uint8Array(15);

	for (let i = 0; i < 15; i++) {

		if ( 0 !== minimumPolynomials[i] ) continue;

		const conjugations = generateConjugations(i);
		const minimumPolynomial = generateMinimumPolynomial(conjugations);

		// 共役な α^i は最小多項式が等しい
		for (const conjugation of conjugations) {
			minimumPolynomials[conjugation] = minimumPolynomial;
		}

	}

	return minimumPolynomials;

};

const BCH_MINIMUM_POLYNOMIALS = generateMinimumPolynomials();

// TODO: 
console.log(Array.from(
	BCH_MINIMUM_POLYNOMIALS,
	(minimumPolynomial, i) => {
		const minimumPolynomialBinary = minimumPolynomial.toString(2).padStart(5, '0');
		return `m_${String(i).padEnd(2)}(x): 0b${minimumPolynomialBinary}`;
	},
).join('\n'));

実行結果

m_0 (x): 0b00011
m_1 (x): 0b10011
m_2 (x): 0b10011
m_3 (x): 0b11111
m_4 (x): 0b10011
m_5 (x): 0b00111
m_6 (x): 0b11111
m_7 (x): 0b11001
m_8 (x): 0b10011
m_9 (x): 0b11111
m_10(x): 0b00111
m_11(x): 0b11001
m_12(x): 0b11111
m_13(x): 0b11001
m_14(x): 0b11001

最小多項式を計算するためにまず数 $\alpha^i$ に対する共役数 $\alpha^{2i}, \alpha^{4i}, \alpha^{8i}, \dots$ を求めます。

ここでは $\alpha^{15} = 1$ より $\alpha^{2^4} = \alpha^{16i} = \alpha^{15i} \alpha^i = \alpha^i$ となるため、各 $i$ に対して最大 4 つの数 $\alpha^i, \alpha^{2i}, \alpha^{4i}, \alpha^{8i}$ を計算すればいいことになります。

実際に $i$ に値を代入すると以下のように共役の関係である数を求められます:

$\{\alpha^1, \alpha^{2\phantom{0}}, \alpha^{4\phantom{0}}, \alpha^{8\phantom{0}}\}$
$\{\alpha^3, \alpha^{6\phantom{0}}, \alpha^{12}, \alpha^{9\phantom{0}}\}$
$\{\alpha^5, \alpha^{10}\}$
$\{\alpha^7, \alpha^{14}, \alpha^{13}, \alpha^{11}\}$

例えば $i = 1$ のとき最小多項式は

m_1(x) = \prod_{i \in \{1, 2, 4, 8\}} (x - \alpha^i)

となり、

m_1(x) = m_2(x) = m_4(x) = m_8(x)

が成り立ちます。

ここでは JavaScript で計算する都合上

\begin{alignat}{2}
&m_{1,-1} & (x) &= 1 \\
&m_{1,0} & (x) &= (x + \alpha^1) \\
&m_{1,1} & (x) &= (x + \alpha^1) (x + \alpha^2) \\
&m_{1,2} & (x) &= (x + \alpha^1) (x + \alpha^2) (x + \alpha^4) \\
&m_{1,3} & (x) &= (x + \alpha^1) (x + \alpha^2) (x + \alpha^4) (x + \alpha^8) \\
& & &= m_1(x)
\end{alignat}

のように計算します。

ここでは最小多項式の係数を 2 進数で表し、例えば最小多項式

\begin{align}
m_1(x) &= x^4 + x + 1 \\
&= 1 x^4 + 0 x^3 + 0 x^2 + 1 x + 1
\end{align}

は 0b10011 になります。

0 または 1 のみを含む配列を 2 進数として単一の数に変換するのは、例えば配列 [1, 0, 0, 1, 1] の場合は以下のように行えます。

(0b00001, 0) => 0b00001 << 1 ^ 0; // 0b00010
(0b00010, 0) => 0b00010 << 1 ^ 0; // 0b00100
(0b00100, 1) => 0b00100 << 1 ^ 1; // 0b01001
(0b01001, 1) => 0b01001 << 1 ^ 1; // 0b10011

3. 生成多項式を計算する

最小多項式から生成多項式を計算します。

/**
 * 生成多項式を計算する
 */
// 意味的には
// G_d(x) = LCM(m_1(x), m_2(x), ... , m_{d-2}(x), m_{d-1}(x))
const bchGenerator = d => {

	// 最小多項式の重複を除去する
	// 意味的には最小公倍数になる
	const factors = new Uint8Array(new Set(
		BCH_MINIMUM_POLYNOMIALS.subarray(1, d - 1),
	));

	// TODO: 
	console.log(Array.from(factors, factor => `0b${factor.toString(2).padStart(5, '0')}`).join(', '));

	// 因数を有限体上で総乗する
	const generator = factors.reduce((generatorPrev, factor) => {

		let generator = 0;

		// 意味的には
		// generator = factor generatorPrev
		for (
			let multiplier = factor, addend = generatorPrev;
			multiplier !== 0;
			multiplier >>>= 1, addend <<= 1
		) {
			if ( multiplier & 1 !== 0 ) {
				generator ^= addend;
			}
		}

		return generator;

	});

	return generator;

};

例

// 
const generator = bchGenerator(7);

console.log(`0b${generator.toString(2)}`);

// 
const eccLength = 31 - Math.clz32(generator);

const n = 15;
const k = n - eccLength;

console.log(`(${n}, ${k})`);

例の実行結果

0b10011, 0b11111, 0b00111
0b10100110111
(15, 5)

ここでは生成多項式は

G_d(x) = \operatorname{LCM} (m_1(x), m_2(x), \dots, m_{d-2}(x), m_{d-1}(x))

より、例えば $d = 7$ のとき

\begin{align}
G_7(x) &= \operatorname{LCM} (m_1(x), m_2(x), m_3(x), m_4(x), m_5(x), m_6(x)) \\
&= m_1(x) m_3(x) m_5(x) \\
&= (x^4 + x + 1) (x^4 + x^3 + x^2 + x + 1) (x^2 + x + 1) \\
&= x^{10} + x^8 + x^5 + x^4 + x^2 + x + 1
\end{align}

のように計算します。

生成多項式も最小多項式と同様に係数を 2 進数で表し、例えば生成多項式

\begin{align}
G_7(x) &= x^{10} + x^8 + x^5 + x^4 + x^2 + x + 1 \\
&= 1 x^{10} + 0 x^9 + 1 x^8 + 0 x^7 + 0 x^6 + 1 x^5 + 1 x^4 + 0 x^3 + 1 x^2 + 1 x + 1
\end{align}

は 0b10100110111 になります。

ここでは符号長は $n = 15$ であり、$d = 7$ のとき生成多項式 $G_7(x)$ は 10 次多項式よりデータ長は

\begin{align}
k &= 15 - 10 \\
&= 5
\end{align}

になります。

$(n, k)$	生成多項式	最小ハミング距離	誤り訂正能力	誤り訂正符号の長さ
$(15,11)$	`0b10011`	3	1	4
$(15, 7)$	`0b111010001`	5	2	8
$(15, 5)$	`0b10100110111`	7	3	10
$(15, 1)$	`0b111111111111111`	15	7	14

4. 符号化

有限体上で情報多項式を生成多項式で割った余りを計算すると、その余りの多項式の係数が誤り訂正符号になります。

// 
const encodeBinaryCyclic = (n, k, generator) => {

	const codeLastOffset = n - 1;
	const dataLastOffset = k - 1;
	const eccLength = n - k;

	return data => {

		const shifted = data << eccLength;

		// 有限体 GF(2^m) 上で多項式の割り算の余りを計算する
		let remainder = shifted;

		for (
			let i = codeLastOffset, subtrahend = generator << dataLastOffset;
			i >= eccLength;
			i--, subtrahend >>>= 1
		) {
			if ( remainder >>> i !== 0 ) {
				remainder ^= subtrahend;
			}
		}

		// data と remainder を合わせて符号とする
		// shifted ^ remainder === data << eccLength ^ remainder
		const code = shifted ^ remainder;

		return code;

	};

};

生成多項式を 0b10100110111 とする (15, 5) 2 元 BCH 符号の場合は、以下のようにして符号化できます。

例

// 
const encodeBinaryBCH = encodeBinaryCyclic(15, 5, 0b10100110111);

// 
const data = 0b00101;
const code = encodeBinaryBCH(data);

console.log(`0b${code.toString(2).padStart(15, '0')}`);

例の実行結果

0b001010011011100

上記の例では、以下のようにして割り算の余りを計算しています。

0b00101_0000000000 ^ 0b1_0100110111 << 2 === 0b00000_0011011100;

5. 総当たりによる誤り訂正

BCH 符号の誤り訂正はいくつかの方法がありますが、データ長が短い場合は総当たりで誤り訂正する方が計算量を減らせます。

// 
const popCount = x => {
	const a = x - (x >>> 1 & 0x55555555);
	const b = (a & 0x33333333) + (a >>> 2 & 0x33333333);
	const c = b + (b >>> 4) & 0x0f0f0f0f;
	const y = c * 0x01010101 >>> 24;
	return y;
};

const correctBinaryCyclic = (t, codes) => code => (
	codes.find(codeCorrected => popCount(code ^ codeCorrected) <= t)
);

例

// 
const generateCodes = () => {

	const codes = new Uint16Array(32);

	for (let data = 0; data < 32; data++) {
		codes[data] = encodeBinaryBCH(data);
	}

	return codes;

};

const correctBinaryBCH = correctBinaryCyclic(3, generateCodes());

// 
const data = 0b00101;
const code = encodeBinaryBCH(data);

console.log(`0b${code.toString(2).padStart(15, '0')}`);

// 
const codeCorrupted = code ^ 0b000100010001000;
const codeCorrected = correctBinaryBCH(codeCorrupted);

console.log(`0b${codeCorrupted.toString(2).padStart(15, '0')}`);
console.log(`0b${codeCorrected.toString(2).padStart(15, '0')}`);

例の実行結果

0b001010011011100
0b001110001010100
0b001010011011100

正しい符号の一覧を作成し、受け取った符号とのハミング距離が誤り訂正能力 (ここでは 3) 以下の正しい符号を探します。

2 つの値のハミング距離は、2 つの値の排他的論理和のハミング重み (ビットカウント) と等しいです。

参考「JavaScript でなるべく速くビットカウントする - Qiita」

6. 関連記事

参考「JavaScript でリード・ソロモン符号: 符号化と誤り検出 - Qiita」
参考「JavaScript で 2 元ゴレイ符号: 符号化と総当たりによる誤り訂正 - Qiita」