リードソロモン符号の理論 — バースト誤りに強い符号を数学から理解する

CDやDVDの表面に小さな傷がついても音楽は正しく再生されます。QRコードの一部がインクで汚れていても、スマートフォンはデータを正しく読み取ります。深宇宙探査機ボイジャーは、太陽系の果てから送られてきた極めて微弱な信号からでも、木星や土星の鮮明な画像を復元できました。これらの「奇跡」の裏側で働いている技術こそが、リードソロモン符号(Reed-Solomon code, RS符号)です。

リードソロモン符号が他の誤り訂正符号と一線を画す最大の特徴は、バースト誤りに強いことです。通信路や記録媒体で発生する誤りは、ランダムに散らばるとは限りません。CDの傷やディスクの汚れは連続した領域に影響し、無線通信のフェージングは一定時間にわたってビットを壊します。このような連続する誤り(バースト誤り)に対してRS符号は驚くほど効果的に機能します。

RS符号を理解すると、以下のような分野への応用が見えてきます。

  • 光ディスク・ストレージ: CD、DVD、Blu-rayのすべてでRS符号が使われており、傷や指紋による連続的な読み取りエラーを訂正しています
  • QRコード: RS符号の誤り訂正レベル(L/M/Q/H)の設定により、最大30%の損傷まで復元可能です
  • 衛星通信・深宇宙通信: NASAの深宇宙ネットワークではRS符号と畳み込み符号の連接符号が長年の標準でした
  • 地上デジタル放送: ISDB-T(日本の地デジ規格)はRS符号をインターリーバと組み合わせて使用しています
  • データ通信: RAID-6のストレージ冗長化やイーサネットの誤り訂正にもRS符号が活用されています

本記事の内容

  • ガロア体 $\mathrm{GF}(2^m)$ の基礎と有限体上の演算
  • 原始多項式と原始元の定義と役割
  • RS符号の構成と生成多項式による符号化
  • シンドロームの計算と誤り訂正の仕組み
  • 誤り訂正能力 $t = (n – k) / 2$ の意味
  • バースト誤りに強い理由の数学的説明
  • BCH符号との関係
  • Python による $\mathrm{GF}(2^8)$ 上のRS符号化・復号シミュレーション

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

ガロア体とは何か — RS符号の舞台を用意する

なぜ有限体が必要なのか

前提記事で扱ったハミング符号や線形符号は、ビットの世界(0と1の世界)で動いていました。加算は排他的論理和(XOR)、乗算は論理積(AND)であり、すべての計算は $\mathrm{GF}(2)$ という最も小さな有限体の上で行われていました。

しかし、RS符号はビット単位ではなくシンボル単位で動きます。1シンボルは $m$ ビットからなり、$2^m$ 通りの値を取ります。たとえば $m = 8$ なら1シンボルは1バイト(0から255の整数)です。RS符号がバースト誤りに強い根本的な理由は、この「シンボル単位の処理」にあります。連続する8ビットの誤りも、RS符号から見ればたった1シンボルの誤りに過ぎないのです。

シンボル同士で加減乗除を行うには、$2^m$ 個の元からなる体(たいすう的な構造)が必要です。これがガロア体 $\mathrm{GF}(2^m)$(Galois Field)です。ガロア体は19世紀のフランスの数学者エヴァリスト・ガロアにちなんで名付けられた有限体で、加減乗除のすべてが閉じた集合として定義されます。

有限体の定義

有限体 $\mathrm{GF}(q)$ とは、有限個の元からなる集合で、以下の演算規則を満たすものです。

$$ (\mathrm{GF}(q),\ +,\ \cdot) $$

  1. 加法に関する群: 加法の単位元 $0$ が存在し、すべての元に加法逆元が存在する
  2. 乗法に関する群: 乗法の単位元 $1$ が存在し、$0$ 以外のすべての元に乗法逆元が存在する(つまり割り算ができる)
  3. 分配法則: $a \cdot (b + c) = a \cdot b + a \cdot c$ が成り立つ

有限体が存在するのは、元の数 $q$ が素数の冪 $p^m$($p$ は素数、$m$ は正の整数)のときに限られるという重要な定理があります。RS符号で使うのは $q = 2^m$、つまり $p = 2$ の場合です。

GF(2^m) の構成方法

$\mathrm{GF}(2)$ は $\{0, 1\}$ の2元で構成される最も単純な体です。しかし $\mathrm{GF}(2^m)$ を作るには、単に $\{0, 1, 2, \dots, 2^m – 1\}$ に通常の整数演算を定義するだけではうまくいきません。たとえば $2 \times 2 = 4$ は $m = 2$ のとき $\mathrm{GF}(4) = \{0, 1, 2, 3\}$ の範囲に収まっても、$2 + 3 = 5$ は $\mathrm{GF}(4)$ からはみ出してしまいます。mod 4 を取ると $5 \bmod 4 = 1$ ですが、この場合 $2$ の逆元が存在しなくなり体の条件を満たせません。

正しい方法は、$\mathrm{GF}(2)$ 上の多項式の剰余環を使って構成することです。具体的には次の手順を踏みます。

  1. $\mathrm{GF}(2)$ 上の $m$ 次既約多項式 $p(x)$ を1つ選ぶ
  2. $\mathrm{GF}(2^m)$ の元を $m – 1$ 次以下の多項式 $a_{m-1}x^{m-1} + \cdots + a_1 x + a_0$(各 $a_i \in \{0, 1\}$)として表す
  3. 加算は多項式の各項を $\mathrm{GF}(2)$(mod 2)で行う
  4. 乗算は多項式の積を $p(x)$ で割った余りとする

この構成により、$\mathrm{GF}(2^m)$ はちょうど $2^m$ 個の元を持つ有限体となります。

GF(2^m) の加算

$\mathrm{GF}(2^m)$ の加算は非常にシンプルです。多項式の各係数を mod 2 で加えるため、これはビット列のXORに一致します。

たとえば $\mathrm{GF}(2^4)$ で、元 $\alpha^3 + \alpha + 1$(ビット表現: 1011)と元 $\alpha^3 + \alpha^2$(ビット表現: 1100)の和は次のようになります。

$$ (\alpha^3 + \alpha + 1) + (\alpha^3 + \alpha^2) = \alpha^2 + \alpha + 1 $$

ビット表現で見ると $1011 \oplus 1100 = 0111$ です。$\mathrm{GF}(2)$ の係数は mod 2 なので $1 + 1 = 0$ となることに注意してください。

また、$\mathrm{GF}(2^m)$ では各元の加法逆元は自分自身です。つまり $a + a = 0$ であり、減算と加算は全く同じ演算です。

$$ a – b = a + b = a \oplus b $$

GF(2^m) の乗算

乗算は加算よりも少し複雑です。まず通常の多項式乗算を行い、次に結果を既約多項式 $p(x)$ で割った余りを取ります。

たとえば $\mathrm{GF}(2^4)$ で既約多項式を $p(x) = x^4 + x + 1$ とすると、$\alpha^2$ と $\alpha^2 + \alpha$ の積は次のように計算します。

まず通常の多項式として掛け合わせます。

$$ \alpha^2 \cdot (\alpha^2 + \alpha) = \alpha^4 + \alpha^3 $$

結果に $\alpha^4$ が現れたため、$m = 4$ 次以上の項が存在します。ここで $p(\alpha) = 0$ の関係、すなわち $\alpha^4 = \alpha + 1$ を用いて $\alpha^4$ を置換します。

$$ \alpha^4 + \alpha^3 = (\alpha + 1) + \alpha^3 = \alpha^3 + \alpha + 1 $$

この結果は $m – 1 = 3$ 次以下の多項式なので、$\mathrm{GF}(2^4)$ の元として有効です。

乗算を効率的に行うために、実用的な実装では対数表(log table)逆対数表(antilog table / exp table)を事前に計算しておきます。$\mathrm{GF}(2^m)$ の $0$ 以外のすべての元は原始元 $\alpha$ の冪で表せるため、$a = \alpha^i$, $b = \alpha^j$ ならば次のように乗算を対数の加算に変換できます。

$$ a \cdot b = \alpha^i \cdot \alpha^j = \alpha^{(i + j) \bmod (2^m – 1)} $$

除算も同様に指数の減算で実現できます。

$$ a / b = \alpha^{(i – j) \bmod (2^m – 1)} $$

ここまでで、RS符号の舞台となるガロア体の演算規則がわかりました。次に、この有限体の構造を決定づける「原始多項式」と「原始元」の概念を見ていきましょう。

原始多項式と原始元

原始元とは何か

ガロア体 $\mathrm{GF}(2^m)$ には $2^m$ 個の元がありますが、$0$ を除く $2^m – 1$ 個の非零元はある特別な構造を持っています。この非零元の集合 $\mathrm{GF}(2^m)^*$ は巡回群(cyclic group)を成し、その生成元を原始元(primitive element)$\alpha$ と呼びます。

つまり、$\alpha$ の冪乗を取ると

$$ \alpha^0 = 1,\ \alpha^1,\ \alpha^2,\ \dots,\ \alpha^{2^m – 2} $$

という $2^m – 1$ 個の元が生成され、これらが $\mathrm{GF}(2^m)^*$ のすべての元を尽くします。さらに $\alpha^{2^m – 1} = 1$(単位元に戻る)となり、巡回構造が閉じます。

直感的に言えば、原始元はガロア体の「回転の種」のようなものです。原始元を繰り返し掛けていくことで、0以外のすべての元を生成できます。この性質があるからこそ、すべての非零元を $\alpha$ の冪として表現でき、前節で述べた「乗算を指数の加算に変換する」テクニックが使えるのです。

原始多項式の定義と役割

原始多項式(primitive polynomial)とは、$\mathrm{GF}(2)$ 上の $m$ 次既約多項式のうち、その根が原始元となるものです。もう少し正確に述べると、$\mathrm{GF}(2)$ 係数の $m$ 次多項式 $p(x)$ が次の条件を満たすとき、$p(x)$ は原始多項式です。

  1. $p(x)$ は $\mathrm{GF}(2)$ 上で既約(因数分解できない)
  2. $p(x)$ の根 $\alpha$($\mathrm{GF}(2^m)$ 内)が原始元である($\alpha$ の位数が $2^m – 1$)

原始多項式は $\mathrm{GF}(2^m)$ の構成そのものを決定します。異なる原始多項式を選ぶと、各元の具体的な表現は変わりますが、体としての代数的構造は同型です。

GF(2^8) の原始多項式の例

RS符号で最もよく使われるのは $\mathrm{GF}(2^8)$ です。1シンボルが1バイト(8ビット)に対応し、$\mathrm{GF}(2^8)$ は $256$ 個の元を持ちます。

$\mathrm{GF}(2^8)$ の代表的な原始多項式は次のものです。

$$ p(x) = x^8 + x^4 + x^3 + x^2 + 1 $$

これは10進表現では $\texttt{0x11D}$($= 285$)に対応します。この多項式を使うと、原始元 $\alpha$($= 2$ に対応)の冪乗は $\alpha^{254}$ まで255個の互いに異なる非零元を生成し、$\alpha^{255} = 1$ に戻ります。

この原始多項式は、QRコードの規格やDVDの誤り訂正などで広く採用されている標準的な選択です。

原始元の冪のテーブル

$\mathrm{GF}(2^8)$ の原始元 $\alpha$ の最初のいくつかの冪を具体的に計算してみましょう。$p(\alpha) = 0$ より $\alpha^8 = \alpha^4 + \alpha^3 + \alpha^2 + 1$ という関係が成り立ちます。

$$ \begin{align} \alpha^0 &= 1 & \text{(= 0x01)} \\ \alpha^1 &= \alpha & \text{(= 0x02)} \\ \alpha^2 &= \alpha^2 & \text{(= 0x04)} \\ \alpha^3 &= \alpha^3 & \text{(= 0x08)} \\ \alpha^7 &= \alpha^7 & \text{(= 0x80)} \\ \alpha^8 &= \alpha^4 + \alpha^3 + \alpha^2 + 1 & \text{(= 0x1D)} \end{align} $$

$\alpha^8$ の計算では、$\alpha^8$ がそのまま $\mathrm{GF}(2^8)$ の元としては存在しないため(8次の項があるため)、$\alpha^8 = \alpha^4 + \alpha^3 + \alpha^2 + 1$ という原始多項式から導かれる関係を使って8ビット以内に「折り返し」ています。これ以降の $\alpha^9, \alpha^{10}, \dots$ も、積を取るたびに8ビットを超えた分を $p(x)$ で割って余りを求めることで計算します。

ガロア体の構造と原始元の性質が理解できたところで、いよいよRS符号の本体に入ります。まずは「RS符号とは何か」をその定義から見ていきましょう。

リードソロモン符号の定義

RS符号の直感的な理解

リードソロモン符号を理解するための最もよいアナロジーは、多項式の補間です。

高校数学で学んだように、$k$ 個の点が与えられれば $k – 1$ 次の多項式が一意に決まります。2点があれば直線(1次多項式)が、3点があれば放物線(2次多項式)が定まるのと同じです。

RS符号のアイデアはこうです。$k$ 個の情報シンボルを $k – 1$ 次の多項式の係数と見なします。この多項式を $n$ 個の異なる点で評価して、$n$ 個のシンボルを得ます。$n > k$ なので $n – k$ 個の冗長シンボルが追加されたことになります。受信側では、$n$ 個のシンボルのうちいくつかが壊れていたとしても、壊れていない $k$ 個以上の点が残っていれば元の多項式を復元できます。

これがRS符号の根本的なアイデアです。「多項式は十分な数の点から一意に復元できる」という数学的事実を、誤り訂正に転用しているのです。

RS符号のパラメータ

RS符号 $\mathrm{RS}(n, k)$ は、$\mathrm{GF}(2^m)$ 上の線形符号で、以下のパラメータを持ちます。

$$ \begin{align} n &= 2^m – 1 & \text{(符号長: 1ブロックのシンボル数)} \\ k & & \text{(情報シンボル数: } 1 \leq k \leq n – 1 \text{)} \\ n – k &= 2t & \text{(冗長シンボル数: 検査シンボル数)} \\ t &= \frac{n – k}{2} & \text{(訂正可能なシンボル誤り数)} \end{align} $$

ここで重要なのは、$n$, $k$, $t$ の単位がすべてシンボルであることです。1シンボルは $m$ ビットに対応するため、ビットレベルでの符号長は $n \times m$ ビットとなります。

たとえば $\mathrm{GF}(2^8)$ 上の $\mathrm{RS}(255, 223)$ は次のようなパラメータを持ちます。

  • 符号長: $n = 255$ シンボル(255バイト)
  • 情報長: $k = 223$ シンボル(223バイト)
  • 冗長シンボル数: $n – k = 32$ シンボル
  • 訂正能力: $t = 16$ シンボル

つまり、255バイトのブロック中で最大16バイトが壊れても、元のデータを完全に復元できます。符号化率は $R = 223/255 \approx 0.875$ であり、約12.5%の冗長性で16シンボルの訂正能力を実現しています。

シングルトン限界とMDS符号

RS符号の最小距離 $d_{\min}$ は次の式で与えられます。

$$ d_{\min} = n – k + 1 = 2t + 1 $$

前提記事で述べたシングルトン限界(Singleton bound)を思い出しましょう。任意の $(n, k)$ 線形符号について、最小距離は次の上限に制限されます。

$$ d_{\min} \leq n – k + 1 $$

RS符号はこの上限を等号で達成します。つまり、与えられた符号長 $n$ と情報長 $k$ のもとで、理論的に達成可能な最大の最小距離を持っています。このような符号をMDS符号(Maximum Distance Separable code)と呼びます。

MDS符号であることの意味は、RS符号が冗長シンボルを「一切無駄なく」誤り訂正能力に変換しているということです。$n – k$ 個の冗長シンボルから得られる最大の訂正能力をフルに引き出しています。

ここまでRS符号の定義とパラメータを見てきました。次は、具体的にどのように符号化を行うのか、生成多項式の構成から見ていきましょう。

RS符号の符号化 — 生成多項式による構成

生成多項式の定義

RS符号の符号化は生成多項式 $g(x)$ に基づいて行われます。$t$ シンボルの誤りを訂正する $\mathrm{RS}(n, k)$ 符号の生成多項式は、原始元 $\alpha$ の連続する $2t$ 個の冪を根として持つ多項式です。

$$ g(x) = \prod_{i=1}^{2t} (x – \alpha^i) = (x – \alpha)(x – \alpha^2) \cdots (x – \alpha^{2t}) $$

$\mathrm{GF}(2^m)$ 上の多項式なので、$-\alpha^i = \alpha^i$($\mathrm{GF}(2^m)$ では $-1 = 1$)であることに注意してください。したがって

$$ g(x) = \prod_{i=1}^{2t} (x + \alpha^i) $$

とも書けます。この生成多項式の次数は $2t = n – k$ であり、$g(x)$ の係数はすべて $\mathrm{GF}(2^m)$ の元です。

生成多項式がなぜこの形なのかを直感的に理解しましょう。$g(x)$ の根が $\alpha, \alpha^2, \dots, \alpha^{2t}$ であるということは、$g(x)$ で割り切れるすべての多項式 $c(x)$ もこれらの根を持つということです。つまり、すべての符号語多項式は $c(\alpha^i) = 0$($i = 1, 2, \dots, 2t$)を満たします。この性質が、後に述べるシンドロームによる誤り検出の数学的基盤となります。

組織符号としての符号化

実用的な実装では、組織符号(systematic code)として符号化するのが一般的です。組織符号では、符号語の先頭 $k$ シンボルが元の情報シンボルそのままであり、末尾の $n – k$ シンボルが検査(パリティ)シンボルとなります。

符号化の手順は次の通りです。

ステップ1: $k$ 個の情報シンボル $m_0, m_1, \dots, m_{k-1}$ から情報多項式を作ります。

$$ m(x) = m_0 x^{k-1} + m_1 x^{k-2} + \cdots + m_{k-1} $$

ステップ2: 情報多項式を $x^{n-k}$ 倍して、検査シンボルの領域を空けます。

$$ x^{n-k} \cdot m(x) $$

ステップ3: $x^{n-k} \cdot m(x)$ を生成多項式 $g(x)$ で割り、余り $r(x)$ を求めます。

$$ x^{n-k} \cdot m(x) = q(x) \cdot g(x) + r(x) $$

ここで $r(x)$ の次数は $n – k – 1$ 以下です。

ステップ4: 符号語多項式 $c(x)$ を次のように構成します。

$$ c(x) = x^{n-k} \cdot m(x) + r(x) $$

$\mathrm{GF}(2^m)$ では $a + a = 0$(減算と加算が同じ)なので、$c(x) = x^{n-k} \cdot m(x) – r(x)$ とも書けます。この構成により $c(x) = q(x) \cdot g(x)$ となるため、$c(x)$ は $g(x)$ で割り切れます。

符号語の構造を見ると、$c(x)$ の上位 $k$ シンボルは $m(x)$ の係数そのもの(情報シンボル)であり、下位 $n – k$ シンボルが $r(x)$ の係数(検査シンボル)です。これが組織符号の意味です。

符号化の具体例

簡単のため、$\mathrm{GF}(2^3)$ 上の $\mathrm{RS}(7, 3)$ を考えましょう。$n = 7$, $k = 3$, $t = 2$ で、最大2シンボルの誤りを訂正できます。

$\mathrm{GF}(2^3)$ の原始多項式を $p(x) = x^3 + x + 1$ とし、原始元を $\alpha$($\alpha^3 = \alpha + 1$)とします。$\mathrm{GF}(2^3)$ の全元は

$$ \{0,\ 1,\ \alpha,\ \alpha^2,\ \alpha^3 = \alpha + 1,\ \alpha^4 = \alpha^2 + \alpha,\ \alpha^5 = \alpha^2 + \alpha + 1,\ \alpha^6 = \alpha^2 + 1\} $$

の8個です。$2t = 4$ なので生成多項式は

$$ g(x) = (x + \alpha)(x + \alpha^2)(x + \alpha^3)(x + \alpha^4) $$

この展開を実際に行ってみましょう。

まず $(x + \alpha)(x + \alpha^2)$ を計算します。

$$ (x + \alpha)(x + \alpha^2) = x^2 + (\alpha^2 + \alpha)x + \alpha^3 = x^2 + \alpha^4 x + \alpha^3 $$

$\alpha + \alpha^2 = \alpha^4$、$\alpha \cdot \alpha^2 = \alpha^3$ という $\mathrm{GF}(2^3)$ の演算を使っています。

次に $(x + \alpha^3)(x + \alpha^4)$ を計算します。

$$ (x + \alpha^3)(x + \alpha^4) = x^2 + (\alpha^3 + \alpha^4)x + \alpha^7 = x^2 + \alpha^5 \cdot x + 1 \cdot x^0 $$

$\alpha^3 + \alpha^4 = (\alpha + 1) + (\alpha^2 + \alpha) = \alpha^2 + 1 = \alpha^6$… いえ、ここを丁寧にやり直しましょう。$\alpha^3 = \alpha + 1$, $\alpha^4 = \alpha^2 + \alpha$ なので

$$ \alpha^3 + \alpha^4 = (\alpha + 1) + (\alpha^2 + \alpha) = \alpha^2 + 1 = \alpha^6 $$

また $\alpha^7 = \alpha^0 = 1$($2^3 – 1 = 7$ なので $\alpha^7 = 1$)です。したがって

$$ (x + \alpha^3)(x + \alpha^4) = x^2 + \alpha^6 x + 1 $$

最後に2つの結果を掛け合わせると4次の生成多項式が得られます。具体的な係数の計算は煩雑ですが、結果は次のようになります。

$$ g(x) = x^4 + \alpha^3 x^3 + x^2 + \alpha x + \alpha^3 $$

この生成多項式を使って符号化を行えば、7シンボルの符号語が得られます。

ここまでで符号化の手順が分かりました。符号語が正しく構成されていれば $g(x)$ で割り切れるという性質がありました。この性質を利用して、受信語に誤りがあるかどうかを判定するのが「シンドローム」です。

シンドロームの計算

受信語と誤りの表現

通信路を通して符号語 $c(x)$ を送信し、受信側で受け取った多項式を $r(x)$ とします。通信路で誤りが発生した場合、受信語は次のように表せます。

$$ r(x) = c(x) + e(x) $$

ここで $e(x)$ は誤り多項式(error polynomial)であり、誤りが発生した位置に対応する係数が非零となります。誤りが $\nu$ シンボルで発生した場合

$$ e(x) = e_{j_1} x^{j_1} + e_{j_2} x^{j_2} + \cdots + e_{j_\nu} x^{j_\nu} $$

ここで $j_1, j_2, \dots, j_\nu$ は誤り位置、$e_{j_l}$ は各位置の誤り値($\mathrm{GF}(2^m)$ の非零元)です。

シンドロームの定義

シンドローム(syndrome)は、受信語を生成多項式の根で評価した値です。$2t$ 個のシンドロームを次のように定義します。

$$ S_i = r(\alpha^i) \quad (i = 1, 2, \dots, 2t) $$

符号語 $c(x)$ は $g(x)$ の根 $\alpha^1, \alpha^2, \dots, \alpha^{2t}$ を共有するため、$c(\alpha^i) = 0$ です。したがって

$$ S_i = r(\alpha^i) = c(\alpha^i) + e(\alpha^i) = e(\alpha^i) $$

シンドロームは誤り多項式のみに依存し、送信された符号語には依存しません。これが非常に重要な性質です。シンドロームが誤りの「症状」を伝えてくれるので、元の符号語を知らなくても誤りを特定できるのです。

シンドロームと誤りの関係

$\nu$ 個の誤りがある場合、シンドロームを誤りの位置と値で表すと次のようになります。

$$ S_i = \sum_{l=1}^{\nu} e_{j_l} (\alpha^{j_l})^i = \sum_{l=1}^{\nu} e_{j_l} X_l^i $$

ここで $X_l = \alpha^{j_l}$ を誤り位置子(error locator)と呼びます。シンドローム方程式を整理すると、次の連立方程式が得られます。

$$ \begin{cases} S_1 = e_{j_1} X_1 + e_{j_2} X_2 + \cdots + e_{j_\nu} X_\nu \\ S_2 = e_{j_1} X_1^2 + e_{j_2} X_2^2 + \cdots + e_{j_\nu} X_\nu^2 \\ \vdots \\ S_{2t} = e_{j_1} X_1^{2t} + e_{j_2} X_2^{2t} + \cdots + e_{j_\nu} X_\nu^{2t} \end{cases} $$

$2t$ 個の方程式に対して、未知数は $\nu$ 個の誤り位置 $X_l$ と $\nu$ 個の誤り値 $e_{j_l}$、合計 $2\nu$ 個です。したがって $2\nu \leq 2t$、つまり $\nu \leq t$ のとき(誤りの数が訂正能力以内のとき)、方程式系は一意に解けます。

全シンドロームがゼロの場合

すべてのシンドロームが $S_1 = S_2 = \cdots = S_{2t} = 0$ である場合、誤りなし($e(x) = 0$)と判定します。正確に言うと、「検出可能な誤りがない」ということです。$t$ 個を超える誤りが偶然にすべてのシンドロームをゼロにする場合は検出不可能ですが、その確率は極めて小さくなります。

シンドロームから誤りの位置と値を求めるアルゴリズムは複数ありますが、ここでは最も基本的な流れを説明します。次のセクションでは、誤り訂正能力の理論的な根拠を見ていきましょう。

誤り訂正能力 t = (n – k) / 2

なぜ t シンボルまで訂正できるのか

RS符号 $\mathrm{RS}(n, k)$ が $t = (n – k) / 2$ シンボルまでの誤りを訂正できる理由を、2つの視点から理解しましょう。

視点1: 最小距離からの議論

前提記事で述べた通り、最小距離 $d_{\min}$ の符号は $\lfloor (d_{\min} – 1) / 2 \rfloor$ 個の誤りを訂正できます。RS符号の最小距離は

$$ d_{\min} = n – k + 1 = 2t + 1 $$

なので

$$ \text{訂正能力} = \left\lfloor \frac{2t + 1 – 1}{2} \right\rfloor = t $$

視点2: シンドローム方程式からの議論

$\nu$ 個の誤りを訂正するには、$\nu$ 個の誤り位置と $\nu$ 個の誤り値、計 $2\nu$ 個の未知数を決定する必要があります。シンドロームは $2t$ 個あるので、$2\nu \leq 2t$、すなわち $\nu \leq t$ のとき一意解が存在します。

この2つの視点は、同じ事実を異なる角度から述べたものです。代数的な構造(最小距離)と、方程式の解の一意性が一致するのは、RS符号がMDS符号であることの美しい帰結です。

誤り訂正と誤り検出のトレードオフ

RS符号の能力は、訂正だけに使う必要はありません。$\rho$ シンボルの誤りを訂正しつつ、$\sigma$ シンボルの誤りを検出する(訂正はできないが存在を知る)ことができ、次の条件を満たす必要があります。

$$ 2\rho + \sigma \leq n – k = 2t $$

たとえば $\mathrm{RS}(255, 223)$($n – k = 32$)では

  • 訂正のみ: $t = 16$ シンボル
  • 検出のみ: $32$ シンボルまで検出可能
  • 混合: 例えば $\rho = 10$ シンボル訂正 + $\sigma = 12$ シンボル検出

この柔軟性は、通信システムの設計において非常に有用です。たとえば、訂正不能な誤りが検出された場合は再送を要求する、というハイブリッド方式(ARQ/FEC組み合わせ)が可能になります。

消失訂正

さらに、消失(erasure)— 誤りの位置はわかっているが値がわからない場合 — に対しては、最大 $n – k = 2t$ シンボルまで訂正できます。消失の場合は位置が既知なので、未知数は誤り値のみとなり、$2t$ 個のシンドローム方程式で $2t$ 個の未知数を解けるからです。

$$ \text{消失訂正能力} = n – k = 2t $$

これは誤り訂正能力の2倍であり、位置情報の価値がいかに大きいかを示しています。CDプレーヤーでは、傷の位置はある程度推定可能なため、消失訂正を活用して実効的な訂正能力を高めています。

誤り訂正能力の理論が分かったところで、次はRS符号の最大の売りである「バースト誤りへの強さ」について掘り下げましょう。

バースト誤りへの強さの理由

ビット単位 vs. シンボル単位

RS符号がバースト誤りに強い理由は、その動作単位がシンボルであることに尽きます。この点をハミング符号と比較して明確にしましょう。

$\mathrm{GF}(2)$ 上のハミング符号 $(7, 4)$ は、7ビットの符号語中で1ビットの誤りを訂正できます。連続する3ビットが壊れたら、もはや訂正不能です。

一方、$\mathrm{GF}(2^8)$ 上の $\mathrm{RS}(255, 223)$ は、255シンボル(= 255バイト)の符号語中で16シンボルの誤りを訂正できます。ここで決定的に重要なのは、1シンボル内のビット誤りが何ビットであろうと、RS符号から見ると1シンボルの誤りとして扱われるという点です。

バースト長の計算

1シンボルが $m$ ビットの RS 符号が $t$ シンボルの誤りを訂正できるとき、訂正可能なバースト誤りの最大長を考えてみましょう。

最良の場合: バースト誤りがシンボル境界にきれいに収まるとき、$t$ シンボル分、つまり $t \times m$ ビットのバースト誤りを訂正できます。

$$ \text{最大バースト長(最良)} = t \times m $$

最悪の場合: バースト誤りがシンボル境界をまたぐとき、先頭と末尾のシンボルは部分的にしか壊れていなくても1シンボル分の誤りとしてカウントされるため、実質的に $t + 1$ シンボルが影響を受けるケースがあります。このときの最大バースト長は

$$ \text{最大バースト長(最悪)} = (t – 1) \times m + 1 $$

$\mathrm{RS}(255, 223)$($m = 8$, $t = 16$)の場合

  • 最良: $16 \times 8 = 128$ ビット(16バイト連続)のバースト誤りを訂正
  • 最悪: $(16 – 1) \times 8 + 1 = 121$ ビットのバースト誤りを訂正

つまり、約120〜128ビット(15〜16バイト)の連続した誤りを訂正できます。ビット単位の符号では考えられない、驚異的なバースト誤り訂正能力です。

インターリーブとの組み合わせ

実際のシステムでは、RS符号をインターリーバ(interleaver)と組み合わせて使うことで、バースト誤り訂正能力をさらに拡大します。

インターリーブの基本的なアイデアは、複数の符号語のシンボルを混ぜ合わせて送信することです。$D$ 個の符号語をインターリーブする場合、送信時にはシンボルを「ラウンドロビン」方式で交互に配置します。

$$ \underbrace{c_1[0], c_2[0], \dots, c_D[0]}_{\text{各符号語の第0シンボル}},\ c_1[1], c_2[1], \dots, c_D[1],\ \dots $$

受信側でデインターリーブ(元の順序に戻す)すると、バースト誤りが $D$ 個の符号語に分散されます。元のバースト誤りの長さが $L$ ビットだとすると、各符号語には最大 $\lceil L / D \rceil$ ビット程度の誤りしか降りかかりません。

深さ $D$ のインターリーブを行った $\mathrm{RS}(255, 223)$ では、訂正可能なバースト長が $D$ 倍に拡大します。

$$ \text{最大バースト長} \approx t \times m \times D = 128D \text{ ビット} $$

CDでは $D = 4$ のインターリーブが使われており、約4,000ビット(約2.5mm分のトラック長)の連続した傷を訂正できます。

なぜバイナリ符号では不十分なのか

同じ冗長度で $\mathrm{GF}(2)$ 上のバイナリ符号を使った場合と比較してみましょう。$\mathrm{RS}(255, 223)$ の冗長ビット数は $(255 – 223) \times 8 = 256$ ビットです。

仮に256ビットの冗長度を持つ $\mathrm{GF}(2)$ 上の線形符号を構成すると、そのバースト誤り訂正能力は高々 $(256 – 1) / 2 \approx 128$ ビット程度に制限されます(Reiger限界)。しかし実際には、ランダム誤りにも対応する必要があるため、バースト誤り専用に冗長度を使い切ることはできず、実効的なバースト訂正能力はさらに低下します。

RS符号は「シンボル単位」という抽象化のおかげで、バースト誤りとランダム誤りの両方に対して均衡の取れた性能を発揮できるのです。

バースト誤りへの強さの理由が分かったところで、RS符号と密接に関連するBCH符号との関係を整理しておきましょう。

BCH符号との関係

BCH符号の概要

BCH符号(Bose-Chaudhuri-Hocquenghem code)は、RS符号と同じく代数的に構成される巡回符号の一族です。1959年にBoseとRay-Chaudhuri、そして独立にHocquenghemが発見しました。RS符号(1960年、ReedとSolomon)はBCH符号のわずか1年後に発表されています。

BCH符号は $\mathrm{GF}(2)$ 上(バイナリ)の符号であり、生成多項式は原始元の連続する冪を根に持つという点でRS符号と構造が似ています。具体的には、$\mathrm{GF}(2)$ 上のBCH符号の生成多項式は

$$ g_{\text{BCH}}(x) = \mathrm{lcm}\left(m_1(x), m_2(x), \dots, m_{2t}(x)\right) $$

ここで $m_i(x)$ は $\alpha^i$ の $\mathrm{GF}(2)$ 上の最小多項式(minimal polynomial)です。

RS符号はBCH符号の一般化

RS符号とBCH符号の関係は次のようにまとめられます。

RS符号は $\mathrm{GF}(2^m)$ 上のBCH符号である。 より正確に言うと、RS符号は符号のアルファベットと根の体が一致するBCH符号です。

BCH符号では、根 $\alpha^i$ は $\mathrm{GF}(2^m)$ の元ですが、生成多項式の係数は $\mathrm{GF}(2)$ に制限されます。そのため、根 $\alpha^i$ を追加するとその共役根($\alpha^{2i}$, $\alpha^{4i}$, … )もすべて根に加えなければならず、生成多項式の次数が大きくなりがちです。

一方、RS符号では符号シンボルが $\mathrm{GF}(2^m)$ の元なので、生成多項式の係数も $\mathrm{GF}(2^m)$ で構いません。共役根を追加する必要がなく、$2t$ 個の根に対して生成多項式の次数がちょうど $2t$ となります。

この違いを表にまとめます。

特性 BCH符号 RS符号
アルファベット $\mathrm{GF}(2)$(ビット) $\mathrm{GF}(2^m)$(シンボル)
生成多項式の係数体 $\mathrm{GF}(2)$ $\mathrm{GF}(2^m)$
生成多項式の次数 $\geq 2t$(共役根のため膨らみうる) 正確に $2t$
最小距離 $\geq 2t + 1$(設計距離) $= 2t + 1$(MDS)
バースト誤り訂正 ビット単位 シンボル単位($m$ ビット束)

設計距離と真の最小距離

BCH符号の最小距離は設計距離 $2t + 1$ 以上であり、真の最小距離がそれより大きい場合があります。これに対し、RS符号は常にシングルトン限界に一致する $d_{\min} = 2t + 1$ を達成します。

BCH符号の方が真の最小距離が大きくなりうるという事実は一見有利に見えますが、これは「生成多項式の次数が大きい = 冗長度が高い」ことの副産物であり、符号化率の点ではRS符号の方が効率的です。

BCH符号との関係を踏まえて、RS符号の理論的な位置づけが明確になりました。次は、ここまでの理論をPythonで実装して確認しましょう。

Pythonでの実装 — GF(2^8) 上のRS符号

GF(2^8) の演算の実装

まず、$\mathrm{GF}(2^8)$ の基本演算(加算、乗算、除算、冪乗)を実装します。効率のため、対数表と逆対数表を事前に計算します。

import numpy as np

class GF256:
    """GF(2^8) の演算クラス(原始多項式: x^8 + x^4 + x^3 + x^2 + 1 = 0x11D)"""

    PRIM_POLY = 0x11D  # x^8 + x^4 + x^3 + x^2 + 1

    def __init__(self):
        # 対数表と逆対数表を構築
        self.exp_table = [0] * 512  # 余裕を持って2倍のサイズ
        self.log_table = [0] * 256

        x = 1
        for i in range(255):
            self.exp_table[i] = x
            self.log_table[x] = i
            x <<= 1  # x * alpha
            if x & 0x100:  # 8ビットを超えたら原始多項式で mod
                x ^= self.PRIM_POLY

        # 周期性のため拡張(mod 255 の演算を簡略化)
        for i in range(255, 512):
            self.exp_table[i] = self.exp_table[i - 255]

    def add(self, a, b):
        """GF(2^8) の加算 = XOR"""
        return a ^ b

    def sub(self, a, b):
        """GF(2^8) の減算 = XOR(加算と同じ)"""
        return a ^ b

    def mul(self, a, b):
        """GF(2^8) の乗算"""
        if a == 0 or b == 0:
            return 0
        return self.exp_table[self.log_table[a] + self.log_table[b]]

    def div(self, a, b):
        """GF(2^8) の除算"""
        if b == 0:
            raise ZeroDivisionError("GF(2^8) でゼロ除算")
        if a == 0:
            return 0
        return self.exp_table[(self.log_table[a] - self.log_table[b]) % 255]

    def power(self, a, n):
        """GF(2^8) の冪乗"""
        if a == 0:
            return 0
        return self.exp_table[(self.log_table[a] * n) % 255]

    def inverse(self, a):
        """GF(2^8) の乗法逆元"""
        if a == 0:
            raise ZeroDivisionError("0の逆元は存在しません")
        return self.exp_table[255 - self.log_table[a]]

gf = GF256()

このコードでは、$\mathrm{GF}(2^8)$ のすべての演算を対数表/逆対数表によるテーブルルックアップで実現しています。乗算は対数の加算、除算は対数の減算に帰着するため、個々の演算が $O(1)$ で実行できます。exp_table を512要素に拡張しているのは、対数の加算結果が255を超える場合にmod演算なしでアクセスできるようにするためです。

GF(2^8) の演算を確認する

対数表と逆対数表が正しく構築されていることを、いくつかの計算で確かめましょう。

# 原始元 alpha = 2 (GF(2^8) での表現)
alpha = 2

# alpha の冪を確認
print("=== alpha の冪乗 ===")
for i in range(10):
    print(f"  alpha^{i} = {gf.exp_table[i]:3d} (0x{gf.exp_table[i]:02X})")

# alpha^8 = alpha^4 + alpha^3 + alpha^2 + 1 = 16 + 8 + 4 + 1 = 29 = 0x1D であることを確認
print(f"\nalpha^8 = {gf.exp_table[8]} (期待値: 29 = 0x1D)")

# 巡回性: alpha^255 = 1 であることを確認
print(f"alpha^255 = {gf.exp_table[255]} (期待値: 1)")

# 演算の確認: alpha^3 * alpha^5 = alpha^8
a = gf.exp_table[3]
b = gf.exp_table[5]
print(f"\nalpha^3 * alpha^5 = {gf.mul(a, b)} (期待値: {gf.exp_table[8]})")

# 逆元: alpha^3 * (alpha^3)^(-1) = 1
a = gf.exp_table[3]
a_inv = gf.inverse(a)
print(f"alpha^3 * (alpha^3)^(-1) = {gf.mul(a, a_inv)} (期待値: 1)")

# 加算: GF(2^8) での加算 = XOR
print(f"\n3 + 5 = {gf.add(3, 5)} (期待値: {3 ^ 5})")
print(f"7 + 7 = {gf.add(7, 7)} (期待値: 0, 自分自身との加算)")

上のコードを実行すると、$\alpha^8 = 29$(= 0x1D = $\alpha^4 + \alpha^3 + \alpha^2 + 1$)であることが確認でき、$\alpha^{255} = 1$ という巡回性も成り立つはずです。乗法逆元の計算も、$a \cdot a^{-1} = 1$ を満たすことが確認できます。

GF(2^8) 上の多項式演算

次に、$\mathrm{GF}(2^8)$ 係数の多項式演算を実装します。RS符号の符号化・復号は、すべて多項式演算に帰着されます。

class GF256Poly:
    """GF(2^8) 上の多項式演算"""

    def __init__(self, gf):
        self.gf = gf

    def poly_add(self, p, q):
        """多項式の加算(= XOR)"""
        r = [0] * max(len(p), len(q))
        for i in range(len(p)):
            r[i + len(r) - len(p)] ^= p[i]
        for i in range(len(q)):
            r[i + len(r) - len(q)] ^= q[i]
        return r

    def poly_mul(self, p, q):
        """多項式の乗算"""
        r = [0] * (len(p) + len(q) - 1)
        for i in range(len(p)):
            for j in range(len(q)):
                r[i + j] ^= self.gf.mul(p[i], q[j])
        return r

    def poly_div(self, dividend, divisor):
        """多項式の除算(商と余りを返す)"""
        result = list(dividend)
        for i in range(len(dividend) - len(divisor) + 1):
            coeff = result[i]
            if coeff != 0:
                for j in range(1, len(divisor)):
                    result[i + j] ^= self.gf.mul(divisor[j], coeff)

        sep = len(dividend) - len(divisor) + 1
        quotient = result[:sep]
        remainder = result[sep:]
        return quotient, remainder

    def poly_eval(self, p, x):
        """多項式を点 x で評価(ホーナー法)"""
        result = p[0]
        for i in range(1, len(p)):
            result = self.gf.add(self.gf.mul(result, x), p[i])
        return result

gfp = GF256Poly(gf)

多項式は係数のリストで表現し、先頭が最高次の係数です。poly_eval はホーナー法を使って効率的に多項式を評価します。ホーナー法は、$p(x) = a_0 x^n + a_1 x^{n-1} + \cdots + a_n$ を $(\cdots((a_0 x + a_1) x + a_2) x + \cdots) x + a_n$ として計算する方法で、乗算回数を最小化します。

RS符号の生成多項式と符号化

生成多項式の構築と符号化を実装します。

class ReedSolomon:
    """GF(2^8) 上のリードソロモン符号"""

    def __init__(self, n, k):
        self.n = n          # 符号長
        self.k = k          # 情報シンボル数
        self.nsym = n - k   # 検査シンボル数 = 2t
        self.t = self.nsym // 2  # 訂正能力
        self.gf = GF256()
        self.gfp = GF256Poly(self.gf)
        self.generator = self._build_generator()

    def _build_generator(self):
        """生成多項式 g(x) = Π(x - alpha^i) for i=1..2t を構築"""
        g = [1]
        for i in range(1, self.nsym + 1):
            # (x - alpha^i) = (x + alpha^i)(GF(2^m)では -1 = 1)
            g = self.gfp.poly_mul(g, [1, self.gf.exp_table[i]])
        return g

    def encode(self, message):
        """組織符号として符号化"""
        if len(message) != self.k:
            raise ValueError(f"情報シンボル数が不正: {len(message)} (期待: {self.k})")

        # x^(n-k) * m(x)
        msg_shifted = message + [0] * self.nsym

        # g(x) で割って余りを求める
        _, remainder = self.gfp.poly_div(msg_shifted, self.generator)

        # 符号語 = 情報シンボル || 検査シンボル
        codeword = message + remainder
        return codeword

    def calc_syndromes(self, received):
        """シンドロームの計算: S_i = r(alpha^i) for i=1..2t"""
        syndromes = []
        for i in range(1, self.nsym + 1):
            s = self.gfp.poly_eval(received, self.gf.exp_table[i])
            syndromes.append(s)
        return syndromes

rs = ReedSolomon(255, 223)
print(f"RS({rs.n}, {rs.k}): 訂正能力 t = {rs.t}")
print(f"生成多項式の次数: {len(rs.generator) - 1}")
print(f"生成多項式の先頭5係数: {rs.generator[:5]}")

_build_generator メソッドは、$g(x) = \prod_{i=1}^{2t}(x + \alpha^i)$ を逐次的に構築しています。最初に $g = [1]$(定数多項式 $1$)から始めて、$(x + \alpha^i)$ を順に掛けていくことで $2t$ 次の生成多項式を得ます。encode メソッドは組織符号の手順に忠実に実装しています。

符号化のデモンストレーション

小さな例で符号化と検証を行ってみましょう。

# RS(255, 223) の符号化デモ
np.random.seed(42)
message = list(np.random.randint(0, 256, size=rs.k))  # ランダムな223バイトのメッセージ

# 符号化
codeword = rs.encode(message)
print(f"メッセージ長: {len(message)} シンボル")
print(f"符号語長:     {len(codeword)} シンボル")
print(f"検査シンボル数: {len(codeword) - len(message)}")

# 符号語の先頭が情報シンボルと一致することを確認(組織符号)
assert codeword[:rs.k] == message, "組織符号の条件が満たされていません"
print("組織符号の確認: OK(先頭 k シンボルが情報シンボルと一致)")

# シンドロームがすべてゼロであることを確認
syndromes = rs.calc_syndromes(codeword)
print(f"シンドローム: {syndromes[:6]}... (すべてゼロなら誤りなし)")
assert all(s == 0 for s in syndromes), "シンドロームが非ゼロです"
print("シンドローム検証: OK(すべてゼロ = 誤りなし)")

このデモでは、ランダムな223バイトのメッセージを符号化して255バイトの符号語を生成しています。組織符号なので符号語の先頭223バイトがメッセージそのものであること、そして正しい符号語のシンドロームがすべてゼロになることを確認しています。

誤り訂正の実装

ここからが本丸です。シンドロームから誤りの位置と値を特定し、訂正するアルゴリズムを実装します。誤り位置多項式の求解にはBerlekamp-Masseyアルゴリズム、誤り位置の特定にはChien探索、誤り値の計算にはForney公式を使います。

def berlekamp_massey(gf, syndromes, nsym):
    """Berlekamp-Masseyアルゴリズム: シンドロームから誤り位置多項式を求める"""
    # 誤り位置多項式 sigma(x)
    sigma = [1]    # 現在の推定
    old_sigma = [1]  # 前回の推定

    for i in range(nsym):
        # 不一致量 delta を計算
        delta = syndromes[i]
        for j in range(1, len(sigma)):
            delta ^= gf.mul(sigma[j], syndromes[i - j])

        old_sigma.append(0)

        if delta != 0:
            if len(old_sigma) > len(sigma):
                # 多項式の更新
                new_sigma = [0] * len(old_sigma)
                for j in range(len(old_sigma)):
                    new_sigma[j] = gf.mul(delta, old_sigma[j])
                for j in range(len(sigma)):
                    new_sigma[j + len(old_sigma) - len(sigma)] ^= sigma[j]
                old_sigma = [gf.mul(gf.inverse(delta), s) for s in sigma]
                # old_sigma のパディング
                while len(old_sigma) < len(new_sigma):
                    old_sigma.append(0)
                sigma = new_sigma
            else:
                # sigma の更新のみ
                new_sigma = list(sigma)
                for j in range(len(old_sigma)):
                    new_sigma[j + len(sigma) - len(old_sigma)] ^= gf.mul(delta, old_sigma[j])
                sigma = new_sigma

    return sigma


def chien_search(gf, sigma, n):
    """Chien探索: 誤り位置多項式の根を全探索で見つける"""
    error_positions = []
    for i in range(n):
        # sigma(alpha^(-i)) = 0 なら位置 i に誤りがある
        val = 0
        for j in range(len(sigma)):
            val ^= gf.mul(sigma[j], gf.power(gf.exp_table[255 - i], j))
        if val == 0:
            error_positions.append(i)
    return error_positions


def forney_algorithm(gf, syndromes, sigma, error_positions, nsym):
    """Forney公式: 誤り値を計算"""
    # 誤り評価多項式 omega(x) = S(x) * sigma(x) mod x^(2t)
    # S(x) = S_1 + S_2*x + ... + S_2t*x^(2t-1)
    s_poly = list(reversed(syndromes))  # 低次から高次

    # sigma を低次から高次に
    sigma_rev = list(reversed(sigma))

    # omega = S(x) * sigma(x)
    omega = [0] * (len(s_poly) + len(sigma_rev) - 1)
    for i in range(len(s_poly)):
        for j in range(len(sigma_rev)):
            omega[i + j] ^= gf.mul(s_poly[i], sigma_rev[j])

    # mod x^nsym で切り詰め
    omega = omega[:nsym]

    # sigma の形式的微分 sigma'(x)(奇数次の項のみ残す)
    sigma_deriv = []
    for i in range(len(sigma_rev)):
        if i % 2 == 1:  # 奇数次
            sigma_deriv.append(sigma_rev[i])
        else:
            sigma_deriv.append(0)
    # 微分なので次数が1つ下がる
    sigma_deriv = sigma_deriv[1:] if len(sigma_deriv) > 1 else [0]

    error_values = []
    for pos in error_positions:
        xi_inv = gf.exp_table[255 - pos]  # X_l^(-1) = alpha^(-pos)

        # omega(X_l^(-1))
        omega_val = 0
        for i in range(len(omega)):
            omega_val ^= gf.mul(omega[i], gf.power(xi_inv, i))

        # sigma'(X_l^(-1))
        sigma_d_val = 0
        for i in range(len(sigma_deriv)):
            sigma_d_val ^= gf.mul(sigma_deriv[i], gf.power(xi_inv, i))

        if sigma_d_val == 0:
            error_values.append(0)
        else:
            # e_l = -omega(X_l^(-1)) / sigma'(X_l^(-1))
            # GF(2^m) では -1 = 1 なので符号は不要
            error_values.append(gf.div(omega_val, sigma_d_val))

    return error_values

Berlekamp-Masseyアルゴリズムは、シンドローム列から最短の線形帰還シフトレジスタ(LFSR)を求めるアルゴリズムです。このLFSRの接続多項式が誤り位置多項式 $\sigma(x)$ に対応します。Chien探索は $\sigma(x)$ の根を $\mathrm{GF}(2^m)$ の全元について試す全探索ですが、$2^m$ 回の評価で済むため実用上十分高速です。Forney公式は誤り位置がわかった後に誤り値を計算する公式で、$\sigma(x)$ の形式的微分を用います。

復号の統合と誤り訂正のシミュレーション

これらを統合して、誤りの挿入・検出・訂正をシミュレーションします。

import matplotlib.pyplot as plt

def rs_decode(rs, received):
    """RS符号の復号: 誤り検出・訂正"""
    syndromes = rs.calc_syndromes(received)

    # シンドロームがすべてゼロなら誤りなし
    if all(s == 0 for s in syndromes):
        return list(received), 0, []

    # Berlekamp-Massey で誤り位置多項式を求める
    sigma = berlekamp_massey(rs.gf, syndromes, rs.nsym)

    # Chien探索で誤り位置を特定
    error_positions = chien_search(rs.gf, sigma, rs.n)

    num_errors = len(sigma) - 1  # 誤り位置多項式の次数 = 誤り数

    if len(error_positions) != num_errors:
        raise ValueError(f"復号失敗: 誤り位置数({len(error_positions)})と"
                        f"誤り位置多項式の次数({num_errors})が不一致")

    # Forney公式で誤り値を計算
    error_values = forney_algorithm(rs.gf, syndromes, sigma, error_positions, rs.nsym)

    # 誤りを訂正
    corrected = list(received)
    for pos, val in zip(error_positions, error_values):
        corrected[pos] ^= val

    return corrected, num_errors, error_positions


# シミュレーション: 誤り数と訂正成功率の関係
rs = ReedSolomon(255, 223)
n_trials = 200
max_errors = 20

success_rates = []

for num_errors in range(max_errors + 1):
    successes = 0
    for trial in range(n_trials):
        # ランダムメッセージ
        msg = list(np.random.randint(0, 256, size=rs.k))
        codeword = rs.encode(msg)

        # ランダムな位置に誤りを挿入
        received = list(codeword)
        error_pos = np.random.choice(rs.n, size=num_errors, replace=False)
        for pos in error_pos:
            error_val = np.random.randint(1, 256)  # 非ゼロの誤り
            received[pos] ^= error_val

        # 復号を試行
        try:
            corrected, _, _ = rs_decode(rs, received)
            if corrected == codeword:
                successes += 1
        except (ValueError, ZeroDivisionError):
            pass  # 復号失敗

    success_rates.append(successes / n_trials)
    print(f"誤り数 {num_errors:2d}: 訂正成功率 {success_rates[-1]:.2%}")

# 可視化
plt.figure(figsize=(10, 6))
plt.bar(range(max_errors + 1), success_rates, color='steelblue', alpha=0.8)
plt.axvline(x=rs.t, color='red', linestyle='--', linewidth=2, label=f't = {rs.t} (理論的訂正限界)')
plt.xlabel('Number of Symbol Errors', fontsize=12)
plt.ylabel('Correction Success Rate', fontsize=12)
plt.title('RS(255, 223) Error Correction Performance', fontsize=14)
plt.legend(fontsize=12)
plt.xticks(range(0, max_errors + 1, 2))
plt.ylim(0, 1.05)
plt.grid(axis='y', alpha=0.3)
plt.tight_layout()
plt.show()

このシミュレーションの結果からは、次の重要な特徴が読み取れます。

  1. 誤り数が $t = 16$ 以下のとき、訂正成功率はほぼ100% — 理論通り、訂正能力内の誤りは確実に訂正されます。これはRS符号が代数的に「完全な」復号アルゴリズムを持つことの裏付けです。

  2. $t = 16$ を超えると急激に訂正成功率が低下する — 赤い破線の位置が理論的な訂正限界であり、この境界を超えると復号は失敗します。ただし、$t + 1$ 個の誤りで偶然正しく復号されるケースがごくまれに生じることがあります。

  3. 性能の劣化は「崖」のように急峻である — 訂正能力の範囲内では完璧に動き、限界を超えるとほぼ完全に失敗するという、非常にシャープな特性を示しています。

バースト誤りの訂正シミュレーション

次に、RS符号のバースト誤り訂正能力を可視化します。

# バースト誤りの訂正シミュレーション
rs = ReedSolomon(255, 223)
n_trials = 300

# バースト長を変化させてテスト
burst_lengths_bytes = range(1, 25)  # バイト単位のバースト長
burst_success_rates = []

for burst_len in burst_lengths_bytes:
    successes = 0
    for trial in range(n_trials):
        msg = list(np.random.randint(0, 256, size=rs.k))
        codeword = rs.encode(msg)
        received = list(codeword)

        # 連続するバースト誤りを挿入
        start = np.random.randint(0, rs.n - burst_len)
        for i in range(burst_len):
            error_val = np.random.randint(1, 256)
            received[start + i] ^= error_val

        try:
            corrected, _, _ = rs_decode(rs, received)
            if corrected == codeword:
                successes += 1
        except (ValueError, ZeroDivisionError):
            pass

    burst_success_rates.append(successes / n_trials)

# 比較: ランダム誤り(同じ数のシンボルが壊れるがバラバラ)
random_success_rates = []
for num_err in burst_lengths_bytes:
    successes = 0
    for trial in range(n_trials):
        msg = list(np.random.randint(0, 256, size=rs.k))
        codeword = rs.encode(msg)
        received = list(codeword)

        positions = np.random.choice(rs.n, size=num_err, replace=False)
        for pos in positions:
            error_val = np.random.randint(1, 256)
            received[pos] ^= error_val

        try:
            corrected, _, _ = rs_decode(rs, received)
            if corrected == codeword:
                successes += 1
        except (ValueError, ZeroDivisionError):
            pass

    random_success_rates.append(successes / n_trials)

plt.figure(figsize=(10, 6))
plt.plot(list(burst_lengths_bytes), burst_success_rates, 'o-', color='steelblue',
         linewidth=2, markersize=6, label='Burst errors (consecutive)')
plt.plot(list(burst_lengths_bytes), random_success_rates, 's--', color='coral',
         linewidth=2, markersize=6, label='Random errors (scattered)')
plt.axvline(x=rs.t, color='red', linestyle=':', linewidth=2, alpha=0.7,
            label=f't = {rs.t} (correction capability)')
plt.xlabel('Number of Corrupted Symbols', fontsize=12)
plt.ylabel('Correction Success Rate', fontsize=12)
plt.title('RS(255, 223): Burst Errors vs Random Errors', fontsize=14)
plt.legend(fontsize=11)
plt.ylim(-0.05, 1.1)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

このグラフから読み取れる重要な知見は以下の通りです。

  1. バースト誤りもランダム誤りも、$t = 16$ シンボルまでほぼ同じ成功率で訂正される — RS符号にとって、誤りが連続しているか散らばっているかは本質的に関係ありません。RS符号は「何シンボルが壊れたか」だけを見ており、誤りのパターンには依存しないのです。

  2. バースト誤りとランダム誤りの訂正成功率カーブがほぼ重なっている — これはRS符号がMDS符号であることの直接的な帰結です。シンボル数が同じなら、どのような誤りパターンでも訂正能力は変わりません。

  3. バースト誤りに「特に強い」というよりは、「バースト誤りが不利にならない」 — ビット単位の符号では長いバーストが致命的になるのに対し、RS符号では$m$ビットの連続誤りが1シンボルに収まるため、同じ長さのバーストでも影響を受けるシンボル数が劇的に少なくなります。

ビット誤り率 (BER) と符号語誤り率 (CWER) のシミュレーション

最後に、二元対称通信路(BSC)を模擬して、ビット誤り率に対する性能を評価します。

# BER vs 符号語誤り率 (CWER) のシミュレーション
rs = ReedSolomon(255, 223)
n_trials = 500

ber_values = np.logspace(-3, -0.5, 15)  # BER: 0.001 ~ 0.3
cwer_coded = []    # RS符号あり
cwer_uncoded = []  # 符号なし

for ber in ber_values:
    coded_errors = 0
    uncoded_errors = 0

    for trial in range(n_trials):
        msg = list(np.random.randint(0, 256, size=rs.k))
        codeword = rs.encode(msg)

        # ビット単位でBSCエラーを発生させる
        codeword_bits = []
        for sym in codeword:
            for bit in range(8):
                codeword_bits.append((sym >> (7 - bit)) & 1)

        # ビット反転
        noise = np.random.random(len(codeword_bits)) < ber
        received_bits = [b ^ int(n) for b, n in zip(codeword_bits, noise)]

        # ビット列をシンボルに戻す
        received = []
        for i in range(0, len(received_bits), 8):
            sym = 0
            for bit in range(8):
                sym = (sym << 1) | received_bits[i + bit]
            received.append(sym)

        # 符号化なしの場合(最初の223シンボルをそのまま比較)
        if received[:rs.k] != msg:
            uncoded_errors += 1

        # RS復号
        try:
            corrected, _, _ = rs_decode(rs, received)
            if corrected[:rs.k] != msg:
                coded_errors += 1
        except (ValueError, ZeroDivisionError):
            coded_errors += 1

    cwer_coded.append(coded_errors / n_trials)
    cwer_uncoded.append(uncoded_errors / n_trials)
    print(f"BER={ber:.4f}: 符号なしCWER={uncoded_errors/n_trials:.3f}, "
          f"RS符号ありCWER={coded_errors/n_trials:.3f}")

plt.figure(figsize=(10, 6))
plt.semilogy(ber_values, cwer_uncoded, 'o--', color='gray', linewidth=2,
             markersize=6, label='Uncoded (raw)')
plt.semilogy(ber_values, [max(c, 1e-4) for c in cwer_coded], 's-',
             color='steelblue', linewidth=2, markersize=6, label='RS(255,223)')
plt.xlabel('Bit Error Rate (BER)', fontsize=12)
plt.ylabel('Codeword Error Rate (CWER)', fontsize=12)
plt.title('RS(255, 223) Performance over BSC', fontsize=14)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3, which='both')
plt.tight_layout()
plt.show()

このBER-CWERグラフからは、RS符号の実用的な効果が明確に読み取れます。

  1. 低BER領域で劇的な改善 — ビット誤り率が $10^{-2}$ 以下の領域では、RS符号により符号語誤り率が大幅に低減されます。符号なしでは頻繁にエラーが発生する通信路でも、RS符号を適用することで高い信頼性を確保できます。

  2. 符号化利得(coding gain)の存在 — 同じCWERを達成するために必要なBERが、符号化によってどれだけ緩和されるかが「符号化利得」であり、RS符号は数dBの利得を提供します。

  3. 高BER領域では限界がある — BERが非常に高い($> 0.1$ 程度)と、16シンボル以上のエラーが頻発するため、RS符号の訂正能力を超えてしまいます。このような場合は、RS符号と畳み込み符号の連接符号や、ターボ符号、LDPC符号などのより強力な符号が必要です。

RS符号の応用事例

ここまでの理論と実装を踏まえて、RS符号が実際にどのように使われているかを具体的に見ておきましょう。

CD/DVD/Blu-ray

光ディスクではCIRC(Cross-Interleaved Reed-Solomon Code)と呼ばれる方式が採用されています。これは2段のRS符号とインターリーバを組み合わせた連接符号です。CDの場合、内符号として $\mathrm{RS}(32, 28)$、外符号として $\mathrm{RS}(28, 24)$ が使われ、深さ4のインターリーブにより約2.5mmの傷(約4,000ビット)を訂正できます。

QRコード

QRコードは $\mathrm{GF}(2^8)$ 上のRS符号を使用し、4段階の誤り訂正レベルを提供します。

  • レベルL: 約7%の損傷を訂正
  • レベルM: 約15%の損傷を訂正
  • レベルQ: 約25%の損傷を訂正
  • レベルH: 約30%の損傷を訂正

レベルHの場合、QRコードの30%が読めなくてもデータを復元できるという事実は、RS符号の訂正能力の驚異的な実例です。デザインQRコード(中央にロゴを入れたもの)が読み取れるのも、この冗長性のおかげです。

衛星通信と深宇宙通信

NASAの深宇宙ネットワークでは、ボイジャーやマーズ・パスファインダーなどの探査機との通信に $\mathrm{RS}(255, 223)$ が長年使われてきました。宇宙空間の通信は極めて低い信号対雑音比(SNR)で行われるため、強力な誤り訂正が不可欠です。現代の深宇宙通信ではターボ符号やLDPC符号に移行しつつありますが、RS符号はその信頼性の高さから今なお広く使われています。

地上デジタル放送(ISDB-T)

日本の地上デジタル放送規格ISDB-Tでは、$\mathrm{RS}(204, 188)$ が外符号として使用されています。内符号の畳み込み符号と組み合わせた連接符号構造により、マルチパスフェージングや都市部のノイズに対する耐性を確保しています。

まとめ

本記事では、リードソロモン符号の理論をガロア体の基礎から丁寧に解説しました。

  • ガロア体 $\mathrm{GF}(2^m)$ は原始多項式と原始元によって構成される有限体であり、RS符号のすべての演算はこの体の上で行われます。加算はXOR、乗算は対数表を用いた指数の加算で効率的に実装できます。

  • RS符号 $\mathrm{RS}(n, k)$ はシンボル単位で動作する線形符号で、$t = (n – k) / 2$ シンボルの誤りを訂正できます。シングルトン限界を等号で達成するMDS符号であり、冗長シンボルを最も効率的に活用しています。

  • 生成多項式 $g(x) = \prod_{i=1}^{2t}(x + \alpha^i)$ による組織符号化は、多項式の剰余演算に帰着されます。受信側ではシンドロームの計算、Berlekamp-Masseyアルゴリズム、Chien探索、Forney公式を組み合わせて誤りを訂正します。

  • バースト誤りへの強さ は、シンボル単位の処理という設計思想から自然に生まれます。$m$ ビットの連続誤りが1シンボルに収まるため、ビット単位の符号と比べて圧倒的に長いバースト誤りを訂正できます。

  • BCH符号との関係 では、RS符号は「アルファベットと根の体が一致するBCH符号」として位置づけられ、共役根を含める必要がないため最も効率的な符号構成が実現されています。

RS符号は1960年の発明以来、60年以上にわたってデジタル通信・記録の信頼性を支え続けてきた、工学史上最も成功した誤り訂正符号の一つです。その理論的な美しさ(MDS性、ガロア体上の代数構造)と実用性(バースト誤り耐性、効率的な実装)の両立が、長寿命の秘訣と言えるでしょう。

次のステップとして、以下の記事も参考にしてください。