bitcoin の bech32 実装 の、わかりにくい関数についてのメモ。
前提
GF(2)
0, 1 の元で構成されるフィールド。
GF(32)
GF(2) の exntension field で、係数が GF(2) の、defining polynomial の次数未満の多項式をエンコードしたものが元になる。上位ビットが上位の次数に対応させて、係数をビットで表してエンコードする。例えば、$x^4+x^2+1$ は 0b10101。
$x^5 + x^3 + 1$ を defining polynomial として、$\mathbb{F}_2[x] / x^5 + x^3 + 1$ の quotient field から元を取っていく。quotient field は、$x^5 + x^3 + 1$ をジェネレーターとして作られるので、$x^5 + x^3 + 1$ はゼロ元に対応する。Primitive element $\alpha$ は、$x^5 + x^3 + 1 = 0$ になるような仮想的な値で、defining polynomial $x^5 + x^3 + 1 = 0$ の ルートとして扱われる。Primitive element はそれ自体を自身に掛けていくことで、GF(32)のゼロを除いた元全体を網羅できる。
元の作りかた
$1$ からはじめて、一つ前の元に $\alpha$ を掛けていく。すると、どこかで defining polynomial の次数と同じになるので、defining polynomial の最大の次数の項以外を、右辺に移項して、次数を下げるための式をつくる。
例えば、$x^5 + x^3 + 1 = 0$ なら、$x^5 = -x^3 - 1$ として、$x^5$ を $-x^3 - 1$ で置き替える。
この作業を続けていくと $1$ に戻る。そこまでに出来たものが、そのフィールドの全ての元で、primitive element を自身に掛け続けることで、全体を網羅できる。
GF(1024)
GF(1024) over GF(32) で、係数が GF(32) の多項式フィールド。defining polynomial は $x^2 + 9x + 23$ で、元の次数は 1 になる。GF32 と同じ要領で元が作られる。
Generator Polynomial
generator polynomial は特定の primitive element $e$ の root を含んでいる。bech32 の generator polynomial は、$e^{997}, e^{998}, e^{999}$ が root になるように作られている。
一般に $\alpha$ が root になる最も低い次数の多項式が minimum polynomial で、bech32 の generator polynomial は、$e^{997}, e^{998}, e^{999}$ の minimum polynomial の最小公倍数の多項式。
この generator polynomial は特定のゴールを満たすために一般的でない作りかたで作られている。
ここ に、一般的な generator polynomial の生成方法についての説明がある。
ある codeword から別の codeword に変換するために必要な変更ビット数 minimum distance が少なくとも $d$ になる BCH code を作るには、primitive element $\alpha$ として、$\alpha^i$ が root になる minimum polynomial を $m_i(x)$ とすると、generator polynomial $g(x)$ は $LCM(m_i(x), m_{i+1}, \cdots, m_{d-1})$ になる。
minium distance $d$ の BCH code は、最大 $d-1$ のエラーを検知することが可能で、最大修正出来るエラー数を $t$ とすると、$d \le 2t + 1$ になる。
また、$GF(2)$ の extension field では $f(x^2)=f(x)^2$ が成り立つので、
$m_1(x) = m_2(x) = m_4(x) \cdots$
$m_3(x) = m_6(x) = m_9(x) \cdots$
これを利用すると generator polynomial の生成計算を簡略化できる。
生成スクリプト
以下の Sage math のスクリプトは、generator polynomial 生成する。
primitive element になる元は複数存在して、スクリプトではランダムに primitive element を探しているので、実行するごと出力結果が変わるが、何回も実行していると $[1, 29, 22, 20, 21, 29, 18]$ を出力することがある。これは、bech32 の generator polynomial。
// The output is a 30-bit integer whose 5-bit groups are the coefficients of the remainder of
// v(x) mod g(x), where g(x) is the Bech32 generator,
// x^6 + {29}x^5 + {22}x^4 + {20}x^3 + {21}x^2 + {29}x + {18}.
# Define the base field GF(2)
GF2 = GF(2)
# Construct the intermediate field GF(32) using a modulus polynomial over GF(2)
R.<x> = PolynomialRing(GF2)
GF32.<a> = GF(2^5, name='a', modulus=x^5 + x^3 + 1)
# Construct the field GF(1024) as an extension of GF(32)
R32.<x> = PolynomialRing(GF32)
defining_poly_32 = x^2 + GF32(9)*x + GF32(23)
GF1024.<b> = GF32.extension(defining_poly_32, 'b')
# Find a primitive element of GF(1024)
def find_primitive_element(field, order):
while True:
elem = field.random_element()
if all(elem ** i != 1 for i in range(1, order)) and elem ** order == 1:
return elem
primitive_element = find_primitive_element(GF1024, 1023)
# Verify that the element found has order 1023
order_check = all(primitive_element ** i != 1 for i in range(1, 1023)) and primitive_element ** 1023 == 1
print(f"Primitive element found has order 1023: {order_check}")
primitive_element_wrap_around = primitive_element ** 1024
pe_check = primitive_element == primitive_element_wrap_around
print(f"Found a cycle at 1024 power: {pe_check}")
# Find the minimal polynomials
m997 = (primitive_element ** 997).minpoly()
m998 = (primitive_element ** 998).minpoly()
m999 = (primitive_element ** 999).minpoly()
generator_polynomial = lcm([m997, m998, m999])
coeffs = [c.integer_representation() for c in generator_polynomial.coefficients()]
coeffs.reverse()
print(f"{generator_polynomial}")
print(f"{coeffs}")
GenerateGFTables()
constexpr std::pair<std::array<int16_t, 1023>, std::array<int16_t, 1024>> GenerateGFTables()
{
std::array<int8_t, 31> GF32_EXP{};
std::array<int8_t, 32> GF32_LOG{};
const int fmod = 41;
GF32_EXP[0] = 1;
GF32_LOG[0] = -1;
GF32_LOG[1] = 0;
int v = 1;
for (int i = 1; i < 31; ++i) {
v = v << 1;
if (v & 32) v ^= fmod;
GF32_EXP[i] = v;
GF32_LOG[v] = i;
}
std::array<int16_t, 1023> GF1024_EXP{};
std::array<int16_t, 1024> GF1024_LOG{};
GF1024_EXP[0] = 1;
GF1024_LOG[0] = -1;
GF1024_LOG[1] = 0;
v = 1;
for (int i = 1; i < 1023; ++i) {
int v0 = v & 31;
int v1 = v >> 5;
int v0n = v1 ? GF32_EXP.at((GF32_LOG.at(v1) + GF32_LOG.at(23)) % 31) : 0;
int v1n = (v1 ? GF32_EXP.at((GF32_LOG.at(v1) + GF32_LOG.at(9)) % 31) : 0) ^ v0;
v = v1n << 5 | v0n;
GF1024_EXP[i] = v;
GF1024_LOG[v] = i;
}
return std::make_pair(GF1024_EXP, GF1024_LOG);
}
なにをやっているか
GF(1024) の Primitive Element を $e$ とすると、$k \rightarrow e^k$ のマッピング GF1024_EXP と、$e^k \rightarrow k$ のマッピング GF1024_LOG を作成している。
実装の詳細
GF32_EXP/LOG 部分
GF(32) は、GF(2) の 5th degree extension で、係数が GF(2) の多項式がその元になる。また、その元は、defining polynomial $x^5 + x^3 + 1$ が定義している。GF(32) の元の多項式は、次数 4 の多項式で、GF(2) の係数を、高い次数から低い次数へ、左から右に並べる形でエンコードされる。
例えば、$x^4 + x + 1$ の場合は、$1 \cdot x^4 + 0 \cdot x^3 + 0 \cdot x^2 + 1 \cdot x^1 + 1 \cdot x^0$ と考えて、10011 がエンコードされた値になる。
この考え方で、defining polynomial をエンコードすると 101001 で、十進数では 41 になる。これを fmod に保存する。
const int fmod = 41;
次に、GF32 の primitive element $\alpha$ が基数 の GF32_EXP/LOG を計算する。任意の数の 0 乗は 1 なので、
GF32_EXP[0] = 1;
GF32_LOG[0] = -1;
GF32_LOG[1] = 0;
で初期化する。$log_{\alpha} 0$ は未定義なので -1 をその意味で置く。
続く for ループでは、$\alpha^1$ から $\alpha^{30}$ までを計算する。$\alpha^{31} = 1$ なので、31 乗は計算しない。
int v = 1;
for (int i = 1; i < 31; ++i) {
v = v << 1;
if (v & 32) v ^= fmod;
GF32_EXP[i] = v;
GF32_LOG[v] = i;
}
v = v << 1;
で、前回の計算結果に $\alpha$ を掛けて、$\alpha$ の乗数を 1 上げる。
if (v & 32) v ^= fmod;
は、$\alpha$ を掛けた結果、元の多項式が次数 5 になった場合に、defining polynomial を足して、$x^5$ の項を相殺し、v を次数 4 以下の多項式に戻す。
そして、結果を保存。
GF32_EXP[i] = v;
GF32_LOG[v] = i;
GF1024_EXP/LOG 部分
GF(1024) は、GF(32) の 2nd-degree extension で、係数が GF(32) の 1 次の多項式が元になる。defining polynomial は、$x^2 + 9x + 23$。primitive element $e$ は、defining polynomial の root なので、$e^2 + 9x + 23 = 0$ で、ここから $x^2 = -9x - 23$ が導き出すことができる。この式を使って次数の reduction を行なう。
まず、GF32_EXP/LOG と同じ要領で、GF1024 の primitive element $e$ が基数 の GF1024_EXP/LOG を計算していく。
任意の数の 0 乗は 1 なので、
GF1024_EXP[0] = 1;
GF1024_LOG[0] = -1;
GF1024_LOG[1] = 0;
続く for ループでは、$e^1$ から $e^{1022}$ までを計算する。$e^{1023} = 1$ なので、1023 乗は計算しない。
v = 1;
for (int i = 1; i < 1023; ++i) {
int v0 = v & 31;
int v1 = v >> 5;
int v0n = v1 ? GF32_EXP.at((GF32_LOG.at(v1) + GF32_LOG.at(23)) % 31) : 0;
int v1n = (v1 ? GF32_EXP.at((GF32_LOG.at(v1) + GF32_LOG.at(9)) % 31) : 0) ^ v0;
v = v1n << 5 | v0n;
GF1024_EXP[i] = v;
GF1024_LOG[v] = i;
}
v は、1 次の多項式 $v_1x + v_0$ の GF(32) の元である 5 ビットの係数 $v_1, v_0$ を格納している。
まず、ここから $v_1, v_0$ を抽出する。
int v0 = v & 31;
int v1 = v >> 5;
次に、$v_1, v_0$ に $e$ を掛ける。
$v = v_1e + v_0$ とすると、$v \cdot e = (v_1e + v_0) \cdot e = v_1e^2 + v_0e$
GF(1024) over GF(32) の元は 1 次の多項式である必要があるので、$e^2 = -9e - 23$ を使って次数を下げると、
$v_1e^2 + v_0e = v_1(-9e - 23) + v_0e = -9v_1e -23v_1 + v_0e = (v_0 - 9v_1)e - 23v_1$
なので、$v \cdot e = v_1'x + v_0'$ とすると、
$v_0' = -23v_1$
$v_1' = v_0 - 9v_1$
GF(32) では、加算と減算はどちらも同じ XOR 演算なので、
$v_0' = 23v_1$
$v_1' = 9v_1 + v_0$
これをコードにすると、
int v0n = v1 ? GF32_EXP.at((GF32_LOG.at(v1) + GF32_LOG.at(23)) % 31) : 0;
int v1n = (v1 ? GF32_EXP.at((GF32_LOG.at(v1) + GF32_LOG.at(9)) % 31) : 0) ^ v0;
最後に、$v_1, v_0$ の 5 ビットの値を連結して、 $v$ に 10 ビットの値として保存する。
v = v1n << 5 | v0n;
GF1024_EXP[i] = v;
GF1024_LOG[v] = i;
PolyMod
uint32_t PolyMod(const data& v)
{
uint32_t c = 1;
for (const auto v_i : v) {
uint8_t c0 = c >> 25;
c = ((c & 0x1ffffff) << 5) ^ v_i;
if (c0 & 1) c ^= 0x3b6a57b2;
if (c0 & 2) c ^= 0x26508e6d;
if (c0 & 4) c ^= 0x1ea119fa;
if (c0 & 8) c ^= 0x3d4233dd;
if (c0 & 16) c ^= 0x2a1462b3;
}
return c;
}
なにをやっているか
5 ビットの値が格納された uint8_t のベクターを受け取り、
多項式 $c$ を $1$ で初期化して、ベクターの先頭から $GF(32)$ の 5 ビットの元を読み込みながら、
- $c$ に $x$ 掛けて次数を一つ上げる
- 読み込んだ元を多項式の定数項にする
- $x^6$ の項を、ジェネレーター多項式で mod したものを $c$ に足す
という形で、$c$ を更新していく。
最後に、$x^5 + v_4x^4 + v_3x^3 + v_2x^2 + v_1x + v_0 $ の、$v_n$ を $v_4 || v_3 || v_2 || v_1 || v_0$ の形で連結したものを返す。各 $v_n$ は 5 ビットの $GF(32)$ の元なので、30 ビットの値が返される。
チェックサムの生成
入力データーを GF(32) の元のベクターに変換してから、最後に 0 を 6 つ追加する。
最後の 6 つの元がゼロなので、入力データー内の元は、全て $x^6$ の位置までシフトされて、ジェネレーター多項式で mod される。
チェックサムの検証
入力データーとチェックサムを連結したものが、PolyMod の入力になる。入力データーとチェックサムが正しい場合には、結果がゼロになる。
仮に入力データーを全てジェネレーター多項式で mod し終った段階の状態を仮に c = 00001 00010 00011 00100 00101 00110 とすると、この c がチェックサムになる。
ここで、このチェックサムを処理していくと、まず、c に $x$ が掛けられて、c = 00010 00011 00100 00101 00110 00000 となり、00001 がこの c に足されて、c = 00010 00011 00100 00101 00110 00001 になる。シフトの結果発生した $x^6$ の項の係数は 00001。
もし、00001 をジェネレーター多項式で mod した結果が 00001 になるのであれば、この結果を c に足すと、c = 00010 00011 00100 00101 00110 00000 となる。これを残りのチェックサムの元について繰り返していくと、最終的に c = 0 になる。
恐らくこういうことなのかなと想像するものの、実際のところ、チェックサムの元を $x^6$ に掛けたものをジェネレーター多項式で mod すると、チェックサムの元になるのかどうかは不明。
実装の詳細
まず、GF(32) の元のベクター $v$ を受け取って、最終的に戻り値になる $c$ を定数多項式 $c1$ で初期化する。
uint32_t c = 1;
次に、$v$ に含まれている GF(32) の元 $v_i$ を使って順に以下の処理で $c$ を更新していく。
uint8_t c0 = c >> 25;
で、$x^5$ の係数を $c0$ に退避して、
c = ((c & 0x1ffffff) << 5) ^ v_i;
で、$c$ から $x^5$ を取り除いて、全ての項の次数を 1 上げる。(= $x$ を掛ける)。次に、その結果に $v_i$ を足す。
if (c0 & 1) c ^= 0x3b6a57b2;
if (c0 & 2) c ^= 0x26508e6d;
if (c0 & 4) c ^= 0x1ea119fa;
if (c0 & 8) c ^= 0x3d4233dd;
if (c0 & 16) c ^= 0x2a1462b3;
ここでは、$c0 \cdot x^6 \bmod g(x)$ を $c$ に足している。
各行では、GF(32) の 5 ビットの元のそれぞれのビットについて、$c0$ のビットがセットされている場合、そのビット掛ける $x^6$ をジェネレーター多項式で mod したものを $c$ に足している。
例えば、c0 & 2 の行は、ジェネレーター多項式を $g(x)$ とすると $2x^6 \bmod g(x)$ の値を $c$ に足している。
GenerateSyndromeConstants
constexpr std::array<uint32_t, 25> GenerateSyndromeConstants() {
std::array<uint32_t, 25> SYNDROME_CONSTS{};
for (int k = 1; k < 6; ++k) {
for (int shift = 0; shift < 5; ++shift) {
int16_t b = GF1024_LOG.at(size_t{1} << shift);
int16_t c0 = GF1024_EXP.at((997*k + b) % 1023);
int16_t c1 = GF1024_EXP.at((998*k + b) % 1023);
int16_t c2 = GF1024_EXP.at((999*k + b) % 1023);
uint32_t c = c2 << 20 | c1 << 10 | c0;
int ind = 5*(k-1) + shift;
SYNDROME_CONSTS[ind] = c;
}
}
return SYNDROME_CONSTS;
}
何をやっているか
シンドロームの計算では、受け取った 5次の多項式の変数部分を、primitive element $e$ の $997, 998, 999$ 乗で置き替えて評価する。
受けとった多項式にエラーが含まれていない場合は、その多項式は、ジェネレーター多項式の倍数で、かつジェネレーター多項式は、$e^{997}, e^{998}, e^{999}$ がルートになる多項式なので、評価結果はゼロになる。
この関数では、受け取った多項式に必要になる、$e^{997}, e^{998}, e^{999}$ の 1 乗から 5 乗までの数の 1, 2, 4, 8, 16 倍 (LSB から見て、1 ビット目から 5 ビット目までに対応) を計算する。
実装の詳細
外側のループの $k$ は 1 から 5 までの乗数を、内側のループの $shift$ は、0 から 4 までのシフト幅。
ループの内側は、
int16_t b = GF1024_LOG.at(size_t{1} << shift);
int16_t c0 = GF1024_EXP.at((997*k + b) % 1023);
int16_t c1 = GF1024_EXP.at((998*k + b) % 1023);
int16_t c2 = GF1024_EXP.at((999*k + b) % 1023);
となっていて c0, c1, c2 はそれぞれ $e^{997}, e^{998}, e^{999}$ の $k$ 乗の、LSB から見た
$b+1$ ビット目に対応する値を計算している。
GF1024_LOG/EXP は、ベースが primitive element $e$ の logarithm を事前に計算して保存したものなので、例えば c0 は、$GF(1024)$ 上での $((e^{997})^k)^{b+1}$ になる。
この計算結果は、
uint32_t c = c2 << 20 | c1 << 10 | c0;
int ind = 5*(k-1) + shift;
SYNDROME_CONSTS[ind] = c;
によって、c にパックされて、SYNDROME_CONSTS に、LSB から MSB、低い次数から高い次数、の順で格納される。
Syndrome
uint32_t Syndrome(const uint32_t residue) {
uint32_t low = residue & 0x1f;
uint32_t result = low ^ (low << 10) ^ (low << 20);
for (int i = 0; i < 25; ++i) {
result ^= ((residue >> (5+i)) & 1 ? SYNDROME_CONSTS.at(i) : 0);
}
return result;
}
何をやっているか
$R(x) = E(x) + C(x)$ の $E(x)$ が residue として渡されてくる。$R(x)$ は受信データー、$C(x)$ はエラーが含まれていない送信されたデーターで、$E(x)$ がエラー部分。
このエラー多項式を、$e^{997}, e^{998}, e^{999}$ で評価して、$e^{999}, e^{998}, e^{997}$ の順に並べて result に格納して返している。
実装の詳細
residue は 5次多項式の、5 ビットの $GF(32)$ の元の 6 つの係数を、$x^5, x^4, x^3, x^2, x^1, x^0$ の順で並べた、30 ビットを格納している。
まず、
uint32_t low = residue & 0x1f;
で、$x^0$ の係数を抽出する。
次に、
uint32_t result = low ^ (low << 10) ^ (low << 20);
result は、$s_{999}$, $s_{998}$, $s_{999}$ を、それぞれ 10 ビットを割りあてながら、その順番で保存する。
このコードは、$x^0$ の係数を、$s_{999}$, $s_{998}$, $s_{999}$ の格納場所にコピーしている。
最後に、
for (int i = 0; i < 25; ++i) {
result ^= ((residue >> (5+i)) & 1 ? SYNDROME_CONSTS.at(i) : 0);
}
が、$x^1$ の最下位 ビット から、$x^5$ の最上位ビットの順で、residue を走査しながら、そのビットに対応する値を result に足していく。
SYNDROME_CONST は、$e^{999}$, $e^{998}$, $e^{999}$ の特定乗の、特定ビットの値を、result と同じ形式でパックしたものなので、そのまま result に足していくとができる。
LocateErrors
std::pair<std::string, std::vector<int>> LocateErrors(const std::string& str) {
...
// str が正しくエンコードされた規定の長さの文字列であること確認後 values ベクターに変換して
// BECH32 と BECH32M エンコーディングの両方で、以下のコードを実行
...
std::vector<int> possible_errors;
uint32_t residue = PolyMod(Cat(ExpandHRP(hrp), values)) ^ EncodingConstant(encoding);
if (residue != 0) {
uint32_t syn = Syndrome(residue);
int s0 = syn & 0x3FF;
int s1 = (syn >> 10) & 0x3FF;
int s2 = syn >> 20;
int l_s0 = GF1024_LOG.at(s0);
int l_s1 = GF1024_LOG.at(s1);
int l_s2 = GF1024_LOG.at(s2);
if (l_s0 != -1 && l_s1 != -1 && l_s2 != -1 && (2 * l_s1 - l_s2 - l_s0 + 2046) % 1023 == 0) {
size_t p1 = (l_s1 - l_s0 + 1023) % 1023;
int l_e1 = l_s0 + (1023 - 997) * p1;
if (p1 < length && !(l_e1 % 33)) {.
possible_errors.push_back(str.size() - p1 - 1);
}
} else {
for (size_t p1 = 0; p1 < length; ++p1) {
int s2_s1p1 = s2 ^ (s1 == 0 ? 0 : GF1024_EXP.at((l_s1 + p1) % 1023));
if (s2_s1p1 == 0) continue;
int l_s2_s1p1 = GF1024_LOG.at(s2_s1p1);
int s1_s0p1 = s1 ^ (s0 == 0 ? 0 : GF1024_EXP.at((l_s0 + p1) % 1023));
if (s1_s0p1 == 0) continue;
int l_s1_s0p1 = GF1024_LOG.at(s1_s0p1);
size_t p2 = (l_s2_s1p1 - l_s1_s0p1 + 1023) % 1023;
if (p2 >= length || p1 == p2) continue;
int s1_s0p2 = s1 ^ (s0 == 0 ? 0 : GF1024_EXP.at((l_s0 + p2) % 1023));
if (s1_s0p2 == 0) continue;
int l_s1_s0p2 = GF1024_LOG.at(s1_s0p2);
int inv_p1_p2 = 1023 - GF1024_LOG.at(GF1024_EXP.at(p1) ^ GF1024_EXP.at(p2));
int l_e2 = l_s1_s0p1 + inv_p1_p2 + (1023 - 997) * p2;
if (l_e2 % 33) continue;
int l_e1 = l_s1_s0p2 + inv_p1_p2 + (1023 - 997) * p1;
if (l_e1 % 33) continue;
if (p1 > p2) {
possible_errors.push_back(str.size() - p1 - 1);
possible_errors.push_back(str.size() - p2 - 1);
} else {
possible_errors.push_back(str.size() - p2 - 1);
possible_errors.push_back(str.size() - p1 - 1);
}
break;
}
}
} else {
return std::make_pair("", std::vector<int>{});
}
}
何をやっているか
- 入力文字列が正しくエンコードされたものかをチェックし、問題があれば失敗で実行終了
- 文字列をデコードして HRP とデーターを取り出して PolyMod 関数を呼び出して $E(x)$ を取得する
- $E(x)=0$ であれば、データーとチェックサムの整合が取れているので、成功で実行終了
- $E(x) \ne 0$ でなければ、入力文字列は不正な文字を含んでいる
- $E(x)$ を入力として、3 つのシンドロームを計算する
- 1 つの文字列が不正なケースかを確認して、もしそうならその文字列の位置を記録する
- そうでなければ、2 つの文字列が不正だという前提で、不正な 2 つの文字列の位置を探す。見つかれば位置を記録する
実装の詳細
チェックサムの確認
まず $R(x) = C(x) + E(x)$ の $E(x)$ を PolyMod で取得する。$R(x)$ は受信したデコード後の文字列、$C(x)$ は送信されたデコード後の文字列、送信された文字列が正しく受信された場合は、$R(x) = C(x)$ になるので、$E(x)$ が正しいデコード後の文字列との差分（エラー部分）になる。$\text{residue} = E(x)$ で、$\text{residue} \ne 0$ で、エラーがある場合にのみ処理をしている。
uint32_t residue = PolyMod(Cat(ExpandHRP(hrp), values)) ^ EncodingConstant(encoding);
if (residue != 0) {
...
シンドロームの計算
次に、$E(x)$ に、$e^{997}, e^{998}, e^{999}$ を適用して、シンドローム $s_{997}, s_{998}, s_{999}$ を計算し、それぞれを $s_0, s_1, s_2$ として保存する。また、それらの $log_{e}$ も計算しておく。$e$ は $GF(1024)$ の primitive element。
uint32_t syn = Syndrome(residue);
int s0 = syn & 0x3FF;
int s1 = (syn >> 10) & 0x3FF;
int s2 = syn >> 20;
int l_s0 = GF1024_LOG.at(s0);
int l_s1 = GF1024_LOG.at(s1);
int l_s2 = GF1024_LOG.at(s2);
不正な文字が 1 つ含まれている場合の処理
この場合の $E(x) = e1 \cdot e^{p1}$ で、不正な文字の正しい係数との差分とその位置を、それぞれ $e1$ と $p1$ が表している。
この $E(x)$ シンドロームを計算すると、$s_{0}=e1 \cdot (e^{997})^{p1}, s_{1}=e1 \cdot (e^{998})^{p1}, s_{2}=e1 \cdot (e^{999})^{p1}$ で、
$\frac{s_1}{s_0} = \frac{e1 \cdot (e^{998})^{p1}}{e1 \cdot (e^{997})^{p1}} = e^{p1}$
$\frac{s_2}{s_1} = \frac{e1 \cdot (e^{999})^{p1}}{e1 \cdot (e^{998})^{p1}} = e^{p1}$
から
$\frac{s_1}{s_0} = \frac{s_2}{s_1}$
$s_1^2 = s_0 \cdot s_2$
$s_1^2 - s_0 \cdot s_2 = 0$
$\log_{e}$ を取ると
$\log_{e} s_1^2 - \log_{e} (s_0 \cdot s_2) = 0$
$2 \log_{e} s_1 - \log_{e} s_0 - \log_{e} s_2 = 0$
これが成り立っていれば、不正な文字が 1 つ含まれていることになるので、コードでチェックする。
if (l_s0 != -1 && l_s1 != -1 && l_s2 != -1 && (2 * l_s1 - l_s2 - l_s0 + 2046) % 1023 == 0) {
...
前半で l_s0, l_s1, l_s2 が未定義になってないこと確認して、後半で上の式が成り立つかを確認している。
$1023 \cdot 2 = 2046$ を 引いているのは、2 つの減算で結果が負の値になること防ぐため。$\bmod 1023$ なので、$2046$ を足しても結果は変わらない。
不正な文字が 1 つ含まれていることが確認できた場合は、その位置と正しい文字との差分を計算して記録する。
$\frac{s_1}{s_0} = e^{p1}$ なので $\log_e$ を取って $\log_e s_1 - \log_e s_0 = p1$ で、$p1$ を抽出できる。
コードにすると、
size_t p1 = (l_s1 - l_s0 + 1023) % 1023;
$1023$ を足しているのは、同じく - l_s0 が計算結果を負の値にするのを防ぐため。
次に
int l_e1 = l_s0 + (1023 - 997) * p1;
で、$e1$ を抽出する。このために、シンドローム $s_0$ の式 $s0 = e1 \cdot (e^{997})^p1$ を変形する。
$s0 = e1 \cdot (e^{997})^{p1}$
$e1 = \frac{s0}{(e^{997})^{p1}} = s0 \cdot (e^{-997})^{p1}$ ... (ア)
$e^{1023} = 1$、$e^{997} \cdot e^{-997} = 1$ なので、
$e^{997} \cdot e^{-997}= 1 = e^{1023}$
$e^{-997} = e^{1023} \cdot e^{-997}$
$e^{-997} = e^{1023 - 997}$
ここで (ア) の式の $e^{-997}$ を、$e^{1023 - 997}$ で置き換えると
$e1 = s0 \cdot (e^{1023 - 997})^{p1}$
$\log_e$ を取って
$\log_e e1 = \log_e s0 \cdot \log_e (e^{1023-997})^{p1} = \log_e s0 \cdot \log_e e^{(1023-997) \cdot p1}$
$\log a^c = c \log a$ なので
$\log_e e1 = \log_e s0 \cdot (1023- 997) \cdot p1 \cdot \log_e e$
$\log_e e1 = \log_e s0 \cdot (1023- 997) \cdot p1$
ここでコードと一致する。
最後に、計算した $p1$ と $e1$ がおかしな値でないかを確認して、問題なければ、位置 $p1$ のみ記録する。
if (p1 < length && !(l_e1 % 33)) {
possible_errors.push_back(str.size() - p1 - 1);
}
$GF(1024), GF(32)$ の加法単位元を除いた元の数は、それぞれ $1023, 31$ で、$GF(1024)$ のジェネレーターを $e$ とすると、$GF(32)$ の元は $e, e^{33}, e^{66}, \cdots, e^{990}$ となっていて、$33$ の間を置いて存在しているらしい。なので $log_e$ を取って $\bmod 33$ が $0$ かどうかを確認することで、正しい $GF(32)$ の元なのかどうかが分かる。
また、データー部分のインデックスに $p1$ がなっていない場合もエラーとしている。
不正な文字が 2 つ含まれている場合の処理
ここではデーター部分の全ての文字について、その文字が不正だと仮定しながら、2 つの不正な文字を見つける計算をする。ただ、どこにその 2 つの不正な文字があるのかは、わからないので、全ての文字について順に、その文字が不正な文字であると仮定して計算を進める。
途中で矛盾が生じた場合は、その文字をスキップして、最後まで辻褄が合った場合には、その 2 つの文字の位置を記録する。
不正な文字が 2 つ含まれている場合の $E(x) = e1 \cdot e^{p1} + e2 \cdot e^{p2}$ でのシンドロームは、
$s0 = E(e^{997}) = e1 \cdot (e^{997})^{p1} + e2 \cdot (e^{997})^{p2}$
$s1 = E(e^{998}) = e1 \cdot (e^{998})^{p1} + e2 \cdot (e^{998})^{p2}$
$s2 = E(e^{999}) = e1 \cdot (e^{999})^{p1} + e2 \cdot (e^{999})^{p2}$
ここで、$s2 + s1 \cdot e^{p1}$ を計算する。もしゼロになってしまった場合は以降の計算は成り立たないので、その文字をスキップする。非ゼロの値が取れた場合は、そのロガリズムも計算しておく。
int s2_s1p1 = s2 ^ (s1 == 0 ? 0 : GF1024_EXP.at((l_s1 + p1) % 1023));
if (s2_s1p1 == 0) continue;
int l_s2_s1p1 = GF1024_LOG.at(s2_s1p1);
また、今後の計算の為に、以下の変形を作っておく。
$s2 + s1 \cdot e^{p1}$
$= e1 \cdot (e^{999})^{p1} + e2 \cdot (e^{999})^{p2} + (e1 \cdot (e^{998})^{p1} + e2 \cdot (e^{998})^{p2}) \cdot e^{p1}$
$= e1 \cdot e^{999p1} + e2 \cdot e^{999p2} + (e1 \cdot e^{998p1} + e2 \cdot e^{998p2}) \cdot e^{p1}$
$= e1 \cdot e^{999p1} + e2 \cdot e^{999p2} + e1 \cdot e^{998p1 + p1} + e2 \cdot e^{998p2 + p1}$
$= e1 \cdot e^{999p1} + e2 \cdot e^{999p2} + e1 \cdot e^{(998 + 1) \cdot p1} + e2 \cdot e^{998p2 + p1}$
$= e1 \cdot e^{999p1} + e2 \cdot e^{999p2} + e1 \cdot e^{999p1} + e2 \cdot e^{998p2 + p1}$
$= 2e1 \cdot e^{999p1} + e2 \cdot e^{999p2} + e2 \cdot e^{998p2 + p1}$
$GF(32)$ では 2 を掛けると 0 になるので、
$= e2 \cdot e^{999p2} + e2 \cdot e^{998p2 + p1}$
$= e2 \cdot e^{998p2} \cdot e^{p2} + e2 \cdot e^{998p2} \cdot e^{p1}$
$= e2 \cdot e^{998p2} (e^{p2} + e^{p1})$ ... (ア)
同様の計算と変形を、$s1, s0$ の組み合わせでも行なうと、
$s1 + s0 \cdot e^{p1}$
$= e2 \cdot e^{997p2} (e^{p2} + e^{p1})$ ... (イ)
コードは、
int s1_s0p1 = s1 ^ (s0 == 0 ? 0 : GF1024_EXP.at((l_s0 + p1) % 1023));
if (s1_s0p1 == 0) continue;
int l_s1_s0p1 = GF1024_LOG.at(s1_s0p1);
もしゼロの値が出れば次の文字にスキップする。
ここで、(ア) と (イ) の比を取ると、$p2$ が抽出できる。
$\frac{s2 + s1 \cdot e^{p1}}{s1 + s0 \cdot e^{p1}} = \frac{e2 \cdot e^{998p2} (e^{p2} + e^{p1})}{e2 \cdot e^{997p2} (e^{p2} + e^{p1})} = e^{p2}$
コードでは、$log_e$ ベースで計算している。1023 を足しているのは、負の数になることを防ぐため。
size_t p2 = (l_s2_s1p1 - l_s1_s0p1 + 1023) % 1023;
ここで $p2$ がおかしな値であれば次の文字にスキップする。
if (p2 >= length || p1 == p2) continue;
次は、$e1, e2$ を求めていく。
まず $s1 + s0 \cdot e^{p2}$ を計算する。やり方は $s2 + s1 \cdot e^{p1}$ と同じ。
結果は、$e1 \cdot e^{997p1} (e^{p1} + e^{p2})$ になる。コードは、
int s1_s0p2 = s1 ^ (s0 == 0 ? 0 : GF1024_EXP.at((l_s0 + p2) % 1023));
if (s1_s0p2 == 0) continue;
int l_s1_s0p2 = GF1024_LOG.at(s1_s0p2);
次に、$s1 + s0 \cdot e^{p1}$ と $s1 + s0 \cdot e^{p2}$ 共通項の逆元 $(e^{p1} + e^{p2})^{-1}$ のロガリズムを用意する。$GF(1024)$ で $e^x$ の逆元は $e^{1023-x}$ なので、
int inv_p1_p2 = 1023 - GF1024_LOG.at(GF1024_EXP.at(p1) ^ GF1024_EXP.at(p2));
ここで、$s1 + s0 \cdot e^{p1}$ に 計算した逆元を掛けると
$(s1 + s0 \cdot e^{p1}) \cdot (e^{p1} + e^{p2})^{-1} = e2 \cdot e^{997p2} (e^{p2} + e^{p1}) \cdot (e^{p1} + e^{p2})^{-1} = e2 \cdot e^{997p2}$
となるので、これを $ e^{997p2}$ で割ることで、$e2$ が抽出できる。
$log_e$ ベースでコードにすると
int l_e2 = l_s1_s0p1 + inv_p1_p2 + (1023 - 997) * p2;
計算した結果が $GF(32)$ の元でない場合は、この文字をスキップする。
if (l_e2 % 33) continue;
同様の処理を $s1 + s0 \cdot e^{p2}$ について行なう。
$(s1 + s0 \cdot e^{p2}) \cdot (e^{p1} + e^{p2})^{-1} = e1 \cdot e^{997p1} (e^{p1} + e^{p2}) \cdot (e^{p1} + e^{p2})^{-1} = e1 \cdot e^{997p1}$
なので、これを $ e^{997p1}$ で割り、$e1$ が抽出する。コードは、
int l_e1 = l_s1_s0p2 + inv_p1_p2 + (1023 - 997) * p1;
同じく、計算した結果が $GF(32)$ の元でなければは、この文字をスキップ。
if (l_e1 % 33) continue;
最後に計算結果を記録する。
if (p1 > p2) {
possible_errors.push_back(str.size() - p1 - 1);
possible_errors.push_back(str.size() - p2 - 1);
} else {
possible_errors.push_back(str.size() - p2 - 1);
possible_errors.push_back(str.size() - p1 - 1);
}