Add references for the generator/constant used in Bech32(m)

src/bech32.cpp

@@ -66,6 +66,26 @@ uint32_t PolyMod(const data& v)
     // the above example, `c` initially corresponds to 1 mod g(x), and after processing 2 inputs of
     // v, it corresponds to x^2 + v0*x + v1 mod g(x). As 1 mod g(x) = 1, that is the starting value
     // for `c`.
+
+    // The following Sage code constructs the generator used:
+    //
+    // B = GF(2) # Binary field
+    // BP.<b> = B[] # Polynomials over the binary field
+    // F_mod = b**5 + b**3 + 1
+    // F.<f> = GF(32, modulus=F_mod, repr='int') # GF(32) definition
+    // FP.<x> = F[] # Polynomials over GF(32)
+    // E_mod = x**2 + F.fetch_int(9)*x + F.fetch_int(23)
+    // E.<e> = F.extension(E_mod) # GF(1024) extension field definition
+    // for p in divisors(E.order() - 1): # Verify e has order 1023.
+    //    assert((e**p == 1) == (p % 1023 == 0))
+    // G = lcm([(e**i).minpoly() for i in range(997,1000)])
+    // print(G) # Print out the generator
+    //
+    // It demonstrates that g(x) is the least common multiple of the minimal polynomials
+    // of 3 consecutive powers (997,998,999) of a primitive element (e) of GF(1024).
+    // That guarantees it is, in fact, the generator of a primitive BCH code with cycle
+    // length 1023 and distance 4. See https://en.wikipedia.org/wiki/BCH_code for more details.
+
     uint32_t c = 1;
     for (const auto v_i : v) {
         // We want to update `c` to correspond to a polynomial with one extra term. If the initial
@@ -88,12 +108,21 @@ uint32_t PolyMod(const data& v)
         // Then compute c1*x^5 + c2*x^4 + c3*x^3 + c4*x^2 + c5*x + v_i:
         c = ((c & 0x1ffffff) << 5) ^ v_i;
 
-        // Finally, for each set bit n in c0, conditionally add {2^n}k(x):
+        // Finally, for each set bit n in c0, conditionally add {2^n}k(x). These constants can be
+        // computed using the following Sage code (continuing the code above):
+        //
+        // for i in [1,2,4,8,16]: # Print out {1,2,4,8,16}*(g(x) mod x^6), packed in hex integers.
+        //     v = 0
+        //     for coef in reversed((F.fetch_int(i)*(G % x**6)).coefficients(sparse=True)):
+        //         v = v*32 + coef.integer_representation()
+        //     print("0x%x" % v)
+        //
         if (c0 & 1)  c ^= 0x3b6a57b2; //     k(x) = {29}x^5 + {22}x^4 + {20}x^3 + {21}x^2 + {29}x + {18}
         if (c0 & 2)  c ^= 0x26508e6d; //  {2}k(x) = {19}x^5 +  {5}x^4 +     x^3 +  {3}x^2 + {19}x + {13}
         if (c0 & 4)  c ^= 0x1ea119fa; //  {4}k(x) = {15}x^5 + {10}x^4 +  {2}x^3 +  {6}x^2 + {15}x + {26}
         if (c0 & 8)  c ^= 0x3d4233dd; //  {8}k(x) = {30}x^5 + {20}x^4 +  {4}x^3 + {12}x^2 + {30}x + {29}
         if (c0 & 16) c ^= 0x2a1462b3; // {16}k(x) = {21}x^5 +     x^4 +  {8}x^3 + {24}x^2 + {21}x + {19}
+
     }
     return c;
 }
@@ -125,7 +154,8 @@ Encoding VerifyChecksum(const std::string& hrp, const data& values)
     // PolyMod computes what value to xor into the final values to make the checksum 0. However,
     // if we required that the checksum was 0, it would be the case that appending a 0 to a valid
     // list of values would result in a new valid list. For that reason, Bech32 requires the
-    // resulting checksum to be 1 instead. In Bech32m, this constant was amended.
+    // resulting checksum to be 1 instead. In Bech32m, this constant was amended. See
+    // https://gist.github.com/sipa/14c248c288c3880a3b191f978a34508e for details.
     const uint32_t check = PolyMod(Cat(ExpandHRP(hrp), values));
     if (check == EncodingConstant(Encoding::BECH32)) return Encoding::BECH32;
     if (check == EncodingConstant(Encoding::BECH32M)) return Encoding::BECH32M;