Inside Bitcoin

1#!/usr/bin/env python3
2# Copyright (c) 2020 Pieter Wuille
3# Distributed under the MIT software license, see the accompanying
4# file LICENSE or http://www.opensource.org/licenses/mit-license.php.
5
6"""Native Python (slow) reimplementation of libminisketch' algorithms."""
7
8import random
9import unittest
10
11# Irreducible polynomials over GF(2) to use (represented as integers).
12#
13# Most fields can be defined by multiple such polynomials. Minisketch uses the one with the minimal
14# number of nonzero coefficients, and tie-breaking by picking the lexicographically first among
15# those.
16#
17# All polynomials for degrees 2 through 64 (inclusive) are given.
18GF2_MODULI = [
19    None, None,
20    2**2 + 2**1 + 1,
21    2**3 + 2**1 + 1,
22    2**4 + 2**1 + 1,
23    2**5 + 2**2 + 1,
24    2**6 + 2**1 + 1,
25    2**7 + 2**1 + 1,
26    2**8 + 2**4 + 2**3 + 2**1 + 1,
27    2**9 + 2**1 + 1,
28    2**10 + 2**3 + 1,
29    2**11 + 2**2 + 1,
30    2**12 + 2**3 + 1,
31    2**13 + 2**4 + 2**3 + 2**1 + 1,
32    2**14 + 2**5 + 1,
33    2**15 + 2**1 + 1,
34    2**16 + 2**5 + 2**3 + 2**1 + 1,
35    2**17 + 2**3 + 1,
36    2**18 + 2**3 + 1,
37    2**19 + 2**5 + 2**2 + 2**1 + 1,
38    2**20 + 2**3 + 1,
39    2**21 + 2**2 + 1,
40    2**22 + 2**1 + 1,
41    2**23 + 2**5 + 1,
42    2**24 + 2**4 + 2**3 + 2**1 + 1,
43    2**25 + 2**3 + 1,
44    2**26 + 2**4 + 2**3 + 2**1 + 1,
45    2**27 + 2**5 + 2**2 + 2**1 + 1,
46    2**28 + 2**1 + 1,
47    2**29 + 2**2 + 1,
48    2**30 + 2**1 + 1,
49    2**31 + 2**3 + 1,
50    2**32 + 2**7 + 2**3 + 2**2 + 1,
51    2**33 + 2**10 + 1,
52    2**34 + 2**7 + 1,
53    2**35 + 2**2 + 1,
54    2**36 + 2**9 + 1,
55    2**37 + 2**6 + 2**4 + 2**1 + 1,
56    2**38 + 2**6 + 2**5 + 2**1 + 1,
57    2**39 + 2**4 + 1,
58    2**40 + 2**5 + 2**4 + 2**3 + 1,
59    2**41 + 2**3 + 1,
60    2**42 + 2**7 + 1,
61    2**43 + 2**6 + 2**4 + 2**3 + 1,
62    2**44 + 2**5 + 1,
63    2**45 + 2**4 + 2**3 + 2**1 + 1,
64    2**46 + 2**1 + 1,
65    2**47 + 2**5 + 1,
66    2**48 + 2**5 + 2**3 + 2**2 + 1,
67    2**49 + 2**9 + 1,
68    2**50 + 2**4 + 2**3 + 2**2 + 1,
69    2**51 + 2**6 + 2**3 + 2**1 + 1,
70    2**52 + 2**3 + 1,
71    2**53 + 2**6 + 2**2 + 2**1 + 1,
72    2**54 + 2**9 + 1,
73    2**55 + 2**7 + 1,
74    2**56 + 2**7 + 2**4 + 2**2 + 1,
75    2**57 + 2**4 + 1,
76    2**58 + 2**19 + 1,
77    2**59 + 2**7 + 2**4 + 2**2 + 1,
78    2**60 + 2**1 + 1,
79    2**61 + 2**5 + 2**2 + 2**1 + 1,
80    2**62 + 2**29 + 1,
81    2**63 + 2**1 + 1,
82    2**64 + 2**4 + 2**3 + 2**1 + 1
83]
84
85class GF2Ops:
86    """Class to perform GF(2^field_size) operations on elements represented as integers.
87
88    Given that elements are represented as integers, addition is simply xor, and not
89    exposed here.
90    """
91
92    def __init__(self, field_size):
93        """Construct a GF2Ops object for the specified field size."""
94        self.field_size = field_size
95        self._modulus = GF2_MODULI[field_size]
96        assert self._modulus is not None
97
98    def mul2(self, x):
99        """Multiply x by 2 in GF(2^field_size)."""
100        x <<= 1
101        if x >> self.field_size:
102            x ^= self._modulus
103        return x
104
105    def mul(self, x, y):
106        """Multiply x by y in GF(2^field_size)."""
107        ret = 0
108        while y:
109            if y & 1:
110                ret ^= x
111            y >>= 1
112            x = self.mul2(x)
113        return ret
114
115    def sqr(self, x):
116        """Square x in GF(2^field_size)."""
117        return self.mul(x, x)
118
119    def inv(self, x):
120        """Compute the inverse of x in GF(2^field_size)."""
121        assert x != 0
122        # Use the extended polynomial Euclidean GCD algorithm on (modulus, x), over GF(2).
123        # See https://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor.
124        t1, t2 = 0, 1
125        r1, r2 = self._modulus, x
126        r1l, r2l = self.field_size + 1, r2.bit_length()
127        while r2:
128            q = r1l - r2l
129            r1 ^= r2 << q
130            t1 ^= t2 << q
131            r1l = r1.bit_length()
132            if r1 < r2:
133                t1, t2 = t2, t1
134                r1, r2 = r2, r1
135                r1l, r2l = r2l, r1l
136        assert r1 == 1
137        return t1
138
139class TestGF2Ops(unittest.TestCase):
140    """Test class for basic arithmetic properties of GF2Ops."""
141
142    def field_size_test(self, field_size):
143        """Test operations for given field_size."""
144
145        gf = GF2Ops(field_size)
146        for i in range(100):
147            x = random.randrange(1 << field_size)
148            y = random.randrange(1 << field_size)
149            x2 = gf.mul2(x)
150            xy = gf.mul(x, y)
151            self.assertEqual(x2, gf.mul(x, 2)) # mul2(x) == x*2
152            self.assertEqual(x2, gf.mul(2, x)) # mul2(x) == 2*x
153            self.assertEqual(xy == 0, x == 0 or y == 0)
154            self.assertEqual(xy == x, y == 1 or x == 0)
155            self.assertEqual(xy == y, x == 1 or y == 0)
156            self.assertEqual(xy, gf.mul(y, x)) # x*y == y*x
157            if i < 10:
158                xp = x
159                for _ in range(field_size):
160                    xp = gf.sqr(xp)
161                self.assertEqual(xp, x) # x^(2^field_size) == x
162            if y != 0:
163                yi = gf.inv(y)
164                self.assertEqual(y == yi, y == 1) # y==1/x iff y==1
165                self.assertEqual(gf.mul(y, yi), 1) # y*(1/y) == 1
166                yii = gf.inv(yi)
167                self.assertEqual(y, yii) # 1/(1/y) == y
168                if x != 0:
169                    xi = gf.inv(x)
170                    xyi = gf.inv(xy)
171                    self.assertEqual(xyi, gf.mul(xi, yi)) # (1/x)*(1/y) == 1/(x*y)
172
173    def test(self):
174        """Run tests."""
175        for field_size in range(2, 65):
176            self.field_size_test(field_size)
177
178# The operations below operate on polynomials over GF(2^field_size), represented as lists of
179# integers:
180#
181#   [a, b, c, ...] = a + b*x + c*x^2 + ...
182#
183# As an invariant, there are never any trailing zeroes in the list representation.
184#
185# Examples:
186# * [] = 0
187# * [3] = 3
188# * [0, 1] = x
189# * [2, 0, 5] = 5*x^2 + 2
190
191def poly_monic(poly, gf):
192    """Return a monic version of the polynomial poly."""
193    # Multiply every coefficient with the inverse of the top coefficient.
194    inv = gf.inv(poly[-1])
195    return [gf.mul(inv, v) for v in poly]
196
197def poly_divmod(poly, mod, gf):
198    """Return the polynomial (quotient, remainder) of poly divided by mod."""
199    assert len(mod) > 0 and mod[-1] == 1 # Require monic mod.
200    if len(poly) < len(mod):
201        return ([], poly)
202    val = list(poly)
203    div = [0 for _ in range(len(val) - len(mod) + 1)]
204    while len(val) >= len(mod):
205        term = val[-1]
206        div[len(val) - len(mod)] = term
207        # If the highest coefficient in val is nonzero, subtract a multiple of mod from it.
208        val.pop()
209        if term != 0:
210            for x in range(len(mod) - 1):
211                val[1 + x - len(mod)] ^= gf.mul(term, mod[x])
212    # Prune trailing zero coefficients.
213    while len(val) > 0 and val[-1] == 0:
214        val.pop()
215    return div, val
216
217def poly_gcd(a, b, gf):
218    """Return the polynomial GCD of a and b."""
219    if len(a) < len(b):
220        a, b = b, a
221    # Use Euclid's algorithm to find the GCD of a and b.
222    # see https://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#Euclid's_algorithm.
223    while len(b) > 0:
224        b = poly_monic(b, gf)
225        (_, b), a = poly_divmod(a, b, gf), b
226    return a
227
228def poly_sqr(poly, gf):
229    """Return the square of polynomial poly."""
230    if len(poly) == 0:
231        return []
232    # In characteristic-2 fields, thanks to Frobenius' endomorphism ((a + b)^2 = a^2 + b^2),
233    # squaring a polynomial is easy: square all the coefficients and interleave with zeroes.
234    # E.g., (3 + 5*x + 17*x^2)^2 = 3^2 + (5*x)^2 + (17*x^2)^2.
235    # See https://en.wikipedia.org/wiki/Frobenius_endomorphism.
236    return [0 if i & 1 else gf.sqr(poly[i // 2]) for i in range(2 * len(poly) - 1)]
237
238def poly_tracemod(poly, param, gf):
239    """Compute y + y^2 + y^4 + ... + y^(2^(field_size-1)) mod poly, where y = param*x."""
240    out = [0, param]
241    for _ in range(gf.field_size - 1):
242        # In each loop iteration i, we start with out = y + y^2 + ... + y^(2^i). By squaring that we
243        # transform it into out = y^2 + y^4 + ... + y^(2^(i+1)).
244        out = poly_sqr(out, gf)
245        # Thus, we just need to add y again to it to get out = y + ... + y^(2^(i+1)).
246        while len(out) < 2:
247            out.append(0)
248        out[1] = param
249        # Finally take a modulus to keep the intermediary polynomials small.
250        _, out = poly_divmod(out, poly, gf)