Add lots of comments to Bech32

src/bech32.cpp

@@ -30,6 +30,33 @@ const int8_t CHARSET_REV[128] = {
      1,  0,  3, 16, 11, 28, 12, 14,  6,  4,  2, -1, -1, -1, -1, -1
 };
 
+// We work with the finite field GF(1024) defined as a degree 2 extension of the base field GF(32)
+// The defining polynomial of the extension is x^2 + 9x + 23
+// Let (e) be a primitive element of GF(1024), that is, a generator of the field.
+// Every non-zero element of the field can then be represented as (e)^k for some power k.
+// The array GF1024_EXP contains all these powers of (e) - GF1024_EXP[k] = (e)^k in GF(1024).
+// Conversely, GF1024_LOG contains the discrete logarithms of these powers, so
+// GF1024_LOG[GF1024_EXP[k]] == k
+// Each element v of GF(1024) is encoded as a 10 bit integer in the following way:
+// v = v1 || v0 where v0, v1 are 5-bit integers (elements of GF(32)).
+//
+// The element (e) is encoded as 9 || 15. Given (v), we compute (e)*(v) by multiplying in the following way:
+// v0' = 27*v1 + 15*v0
+// v1' = 6*v1 + 9*v0
+// e*v = v1' || v0'
+//
+// The following sage code can be used to reproduce both _EXP and _LOG arrays
+// GF1024_LOG = [-1] + [0] * 1023
+// GF1024_EXP = [1] * 1024
+// v = 1
+// for i in range(1, 1023):
+//     v0 = v & 31
+//     v1 = v >> 5
+//     v0n = F.fetch_int(27)*F.fetch_int(v1) + F.fetch_int(15)*F.fetch_int(v0)
+//     v1n = F.fetch_int(6)*F.fetch_int(v1) + F.fetch_int(9)*F.fetch_int(v0)
+//     v = v1n.integer_representation() << 5 | v0n.integer_representation()
+//     GF1024_EXP[i] = v
+//     GF1024_LOG[v] = i
 
 const int16_t GF1024_EXP[] = {
     1, 303, 635, 446, 997, 640, 121, 142, 959, 420, 350, 438, 166, 39, 543,
@@ -103,6 +130,7 @@ const int16_t GF1024_EXP[] = {
     433, 610, 116, 855, 180, 479, 910, 1014, 599, 915, 905, 306, 516, 731,
     626, 978, 825, 344, 605, 654, 209
 };
+// As above, GF1024_EXP contains all elements of GF(1024) except 0
 static_assert(std::size(GF1024_EXP) == 1023, "GF1024_EXP length should be 1023");
 
 const int16_t GF1024_LOG[] = {
@@ -276,8 +304,51 @@ uint32_t PolyMod(const data& v)
     return c;
 }
 
+/** Syndrome computes the values s_j = R(e^j) for j in [997, 998, 999]. As described above, the
+ * generator polynomial G is the LCM of the minimal polynomials of (e)^997, (e)^998, and (e)^999.
+ *
+ * Consider a codeword with errors, of the form R(x) = C(x) + E(x). The residue is the bit-packed
+ * result of computing R(x) mod G(X), where G is the generator of the code. Because C(x) is a valid
+ * codeword, it is a multiple of G(X), so the residue is in fact just E(x) mod G(x). Note that all
+ * of the (e)^j are roots of G(x) by definition, so R((e)^j) = E((e)^j).
+ *
+ * Syndrome returns the three values packed into a 30-bit integer, where each 10 bits is one value.
+ */
 uint32_t Syndrome(const uint32_t residue) {
+    // Let R(x) = r1*x^5 + r2*x^4 + r3*x^3 + r4*x^2 + r5*x + r6
+    // low is the first 5 bits, corresponding to the r6 in the residue
+    // (the constant term of the polynomial).
+
     uint32_t low = residue & 0x1f;
+
+    // Recall that XOR corresponds to addition in a characteristic 2 field.
+    //
+    // To compute R((e)^j), we are really computing:
+    // r1*(e)^(j*5) + r2*(e)^(j*4) + r3*(e)^(j*3) + r4*(e)^(j*2) + r5*(e)^j + r6
+    // Now note that all of the (e)^(j*i) for i in [5..0] are constants and can be precomputed
+    // for efficiency. But even more than that, we can consider each coefficient as a bit-string.
+    // For example, take r5 = (b_5, b_4, b_3, b_2, b_1) written out as 5 bits. Then:
+    // r5*(e)^j = b_1*(e)^j + b_2*(2*(e)^j) + b_3*(4*(e)^j) + b_4*(8*(e)^j) + b_5*(16*(e)^j)
+    // where all the (2^i*(e)^j) are constants and can be precomputed. Then we just add each
+    // of these corresponding constants to our final value based on the bit values b_i.
+    // This is exactly what is done below. Note that all three values of s_j for j in (997, 998,
+    // 999) are computed simultaneously.
+    //
+    // We begin by setting s_j = low = r6 for all three values of j, because these are unconditional.
+    // Then for each following bit, we add the corresponding precomputed constant if the bit is 1.
+    // For example, 0x31edd3c4 is 1100011110 1101110100 1111000100 when unpacked in groups of 10
+    // bits, corresponding exactly to a^999 || a^998 || a^997 (matching the corresponding values in
+    // GF1024_EXP above).
+    //
+    // The following sage code reproduces these constants:
+    // for k in range(1, 6):
+    //     for b in [1,2,4,8,16]:
+    //         c0 = GF1024_EXP[(997*k + GF1024_LOG[b]) % 1023]
+    //         c1 = GF1024_EXP[(998*k + GF1024_LOG[b]) % 1023]
+    //         c2 = GF1024_EXP[(999*k + GF1024_LOG[b]) % 1023]
+    //         c = c2 << 20 | c1 << 10 | c0
+    //         print("0x%x" % c)
+
     return low ^ (low << 10) ^ (low << 20) ^
         ((residue >> 5) & 1 ? 0x31edd3c4 : 0) ^
         ((residue >> 6) & 1 ? 0x335f86a8 : 0) ^
@@ -469,44 +540,104 @@ std::string LocateErrors(const std::string& str, std::vector<int>& error_locatio
     // We can't simply use the segwit version, because that may be one of the errors
     for (Encoding encoding : {Encoding::BECH32, Encoding::BECH32M}) {
         std::vector<int> possible_errors;
+        // Recall that (ExpandHRP(hrp) ++ values) is interpreted as a list of coefficients of a polynomial
+        // over GF(32). PolyMod computes the "remainder" of this polynomial modulo the generator G(x).
         uint32_t residue = PolyMod(Cat(ExpandHRP(hrp), values)) ^ EncodingConstant(encoding);
+
+        // All valid codewords should be multiples of G(x), so this remainder (after XORing with the encoding
+        // constant) should be 0 - hence 0 indicates there are no errors present.
         if (residue != 0) {
+            // If errors are present, our polynomial must be of the form C(x) + E(x) where C is the valid
+            // codeword (a multiple of G(x)), and E encodes the errors.
             uint32_t syn = Syndrome(residue);
+
+            // Unpack the three 10-bit syndrome values
             int s0 = syn & 0x3FF;
             int s1 = (syn >> 10) & 0x3FF;
             int s2 = syn >> 20;
+
+            // Get the discrete logs of these values in GF1024 for more efficient computation
             int l_s0 = GF1024_LOG[s0];
             int l_s1 = GF1024_LOG[s1];
             int l_s2 = GF1024_LOG[s2];
 
+            // First, suppose there is only a single error. Then E(x) = e1*x^p1 for some position p1
+            // Then s0 = E((e)^997) = e1*(e)^(997*p1) and s1 = E((e)^998) = e1*(e)^(998*p1)
+            // Therefore s1/s0 = (e)^p1, and by the same logic, s2/s1 = (e)^p1 too.
+            // Hence, s1^2 == s0*s2, which is exactly the condition we check first:
             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;
+                // Compute the error position p1 as l_s1 - l_s0 = p1 (mod 1023)
+                size_t p1 = (l_s1 - l_s0 + 1023) % 1023; // the +1023 ensures it is positive
+                // Now because s0 = e1*(e)^(997*p1), we get e1 = s0/((e)^(997*p1)). Remember that (e)^1023 = 1,
+                // so 1/((e)^997) = (e)^(1023-997).
                 int l_e1 = l_s0 + (1023 - 997) * p1;
+                // Finally, some sanity checks on the result:
+                // - The error position should be within the length of the data
+                // - e1 should be in GF(32), which implies that e1 = (e)^(33k) for some k (the 31 non-zero elements
+                // of GF(32) form an index 33 subgroup of the 1023 non-zero elements of GF(1024)).
                 if (p1 < length && !(l_e1 % 33)) {
+                    // Polynomials run from highest power to lowest, so the index p1 is from the right.
+                    // We don't return e1 because it is dangerous to suggest corrections to the user,
+                    // the user should check the address themselves.
                     possible_errors.push_back(str.size() - p1 - 1);
                 }
+            // Otherwise, suppose there are two errors. Then E(x) = e1*x^p1 + e2*x^p2.
             } else {
+                // For all possible first error positions p1
                 for (size_t p1 = 0; p1 < length; ++p1) {
+                    // We have guessed p1, and want to solve for p2. Recall that E(x) = e1*x^p1 + e2*x^p2, so
+                    // s0 = E((e)^997) = e1*(e)^(997^p1) + e2*(e)^(997*p2), and similar for s1 and s2.
+                    //
+                    // Consider s2 + s1*(e)^p1
+                    //          = 2e1*(e)^(999^p1) + e2*(e)^(999*p2) + e2*(e)^(998*p2)*(e)^p1
+                    //          = e2*(e)^(999*p2) + e2*(e)^(998*p2)*(e)^p1
+                    //    (Because we are working in characteristic 2.)
+                    //          = e2*(e)^(998*p2) ((e)^p2 + (e)^p1)
+                    //
                     int s2_s1p1 = s2 ^ (s1 == 0 ? 0 : GF1024_EXP[(l_s1 + p1) % 1023]);
                     if (s2_s1p1 == 0) continue;
+                    int l_s2_s1p1 = GF1024_LOG[s2_s1p1];
 
+                    // Similarly, s1 + s0*(e)^p1
+                    //          = e2*(e)^(997*p2) ((e)^p2 + (e)^p1)
                     int s1_s0p1 = s1 ^ (s0 == 0 ? 0 : GF1024_EXP[(l_s0 + p1) % 1023]);
                     if (s1_s0p1 == 0) continue;
-
                     int l_s1_s0p1 = GF1024_LOG[s1_s0p1];
-                    size_t p2 = (GF1024_LOG[s2_s1p1] - l_s1_s0p1 + 1023) % 1023;
+
+                    // So, putting these together, we can compute the second error position as
+                    // (e)^p2 = (s2 + s1^p1)/(s1 + s0^p1)
+                    // p2 = log((e)^p2)
+                    size_t p2 = (l_s2_s1p1 - l_s1_s0p1 + 1023) % 1023;
+
+                    // Sanity checks that p2 is a valid position and not the same as p1
                     if (p2 >= length || p1 == p2) continue;
 
+                    // Now we want to compute the error values e1 and e2.
+                    // Similar to above, we compute s1 + s0*(e)^p2
+                    //          = e1*(e)^(997*p1) ((e)^p1 + (e)^p2)
                     int s1_s0p2 = s1 ^ (s0 == 0 ? 0 : GF1024_EXP[(l_s0 + p2) % 1023]);
                     if (s1_s0p2 == 0) continue;
+                    int l_s1_s0p2 = GF1024_LOG[s1_s0p2];
 
+                    // And compute (the log of) 1/((e)^p1 + (e)^p2))
                     int inv_p1_p2 = 1023 - GF1024_LOG[GF1024_EXP[p1] ^ GF1024_EXP[p2]];
+
+                    // Then (s1 + s0*(e)^p1) * (1/((e)^p1 + (e)^p2)))
+                    //         = e2*(e)^(997*p2)
+                    // Then recover e2 by dividing by (e)^(997*p2)
                     int l_e2 = l_s1_s0p1 + inv_p1_p2 + (1023 - 997) * p2;
+                    // Check that e2 is in GF(32)
                     if (l_e2 % 33) continue;
 
-                    int l_e1 = GF1024_LOG[s1_s0p2] + inv_p1_p2 + (1023 - 997) * p1;
+                    // In the same way, (s1 + s0*(e)^p2) * (1/((e)^p1 + (e)^p2)))
+                    //         = e1*(e)^(997*p1)
+                    // So recover e1 by dividing by (e)^(997*p1)
+                    int l_e1 = l_s1_s0p2 + inv_p1_p2 + (1023 - 997) * p1;
+                    // Check that e1 is in GF(32)
                     if (l_e1 % 33) continue;
 
+                    // Again, we do not return e1 or e2 for safety.
+                    // Order the error positions from the left of the string and return them
                     if (p1 > p2) {
                         possible_errors.push_back(str.size() - p1 - 1);
                         possible_errors.push_back(str.size() - p2 - 1);