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);
}