Inside Bitcoin

251    return out
252
253def poly_frobeniusmod(poly, gf):
254    """Compute x^(2^field_size) mod poly."""
255    out = [0, 1]
256    for _ in range(gf.field_size):
257        _, out = poly_divmod(poly_sqr(out, gf), poly, gf)
258    return out
259
260def poly_find_roots(poly, gf):
261    """Find the roots of poly if fully factorizable with unique roots, [] otherwise."""
262    assert len(poly) > 0
263    # If the polynomial is constant (and nonzero), it has no roots.
264    if len(poly) == 1:
265        return []
266    # Make the polynomial monic (which doesn't change its roots).
267    poly = poly_monic(poly, gf)
268    # If the polynomial is of the form x+a, return a.
269    if len(poly) == 2:
270        return [poly[0]]
271    # Otherwise, first test that poly can be completely factored into unique roots. The polynomial
272    # x^(2^fieldsize)-x has every field element once as root. Thus we want to know that that is a
273    # multiple of poly. Compute x^(field_size) mod poly, which needs to equal x if that is the case
274    # (unless poly has degree <= 1, but that case is handled above).
275    if poly_frobeniusmod(poly, gf) != [0, 1]:
276        return []
277
278    def rec_split(poly, randv):
279        """Recursively split poly using the Berlekamp trace algorithm."""
280        # See https://hal.archives-ouvertes.fr/hal-00626997/document.
281        assert len(poly) > 1 and poly[-1] == 1 # Require a monic poly.
282        # If poly is of the form x+a, its root is a.
283        if len(poly) == 2:
284            return [poly[0]]
285        # Try consecutive randomization factors randv, until one is found that factors poly.
286        while True:
287            # Compute the trace of (randv*x) mod poly. This is a polynomial that maps half of the
288            # domain to 0, and the other half to 1. Which half that is is controlled by randv.
289            # By taking it modulo poly, we only add a multiple of poly. Thus the result has at least
290            # the shared roots of the trace polynomial and poly still, but may have others.
291            trace = poly_tracemod(poly, randv, gf)
292            # Using the set {2^i*a for i=0..fieldsize-1} gives optimally independent randv values
293            # (no more than fieldsize are ever needed).
294            randv = gf.mul2(randv)
295            # Now take the GCD of this trace polynomial with poly. The result is a polynomial
296            # that only has the shared roots of the trace polynomial and poly as roots.
297            gcd = poly_gcd(trace, poly, gf)
298            # If the result has a degree higher than 1, and lower than that of poly, we found a
299            # useful factorization.
300            if len(gcd) != len(poly) and len(gcd) > 1:
301                break
302            # Otherwise, continue with another randv.
303        # Find the actual factors: the monic version of the GCD above, and poly divided by it.
304        factor1 = poly_monic(gcd, gf)
305        factor2, _ = poly_divmod(poly, gcd, gf)
306        # Recurse.
307        return rec_split(factor1, randv) + rec_split(factor2, randv)
308
309    # Invoke the recursive splitting with a random initial factor, and sort the results.
310    return sorted(rec_split(poly, random.randrange(1, 1 << gf.field_size)))
311
312class TestPolyFindRoots(unittest.TestCase):
313    """Test class for poly_find_roots."""
314
315    def field_size_test(self, field_size):
316        """Run tests for given field_size."""
317        gf = GF2Ops(field_size)
318        for test_size in [0, 1, 2, 3, 10]:
319            roots = [random.randrange(1 << field_size) for _ in range(test_size)]
320            roots_set = set(roots)
321            # Construct a polynomial with all elements of roots as roots (with multiplicity).
322            poly = [1]
323            for root in roots:
324                new_poly = [0] + poly
325                for n, c in enumerate(poly):
326                    new_poly[n] ^= gf.mul(c, root)
327                poly = new_poly
328            # Invoke the root finding algorithm.
329            found_roots = poly_find_roots(poly, gf)
330            # The result must match the input, unless any roots were repeated.
331            if len(roots) == len(roots_set):
332                self.assertEqual(found_roots, sorted(roots))
333            else:
334                self.assertEqual(found_roots, [])
335
336    def test(self):
337        """Run tests."""
338        for field_size in range(2, 65):
339            self.field_size_test(field_size)
340
341def berlekamp_massey(syndromes, gf):
342    """Implement the Berlekamp-Massey algorithm.
343
344    Takes as input a sequence of GF(2^field_size) elements, and returns the shortest LSFR
345    that generates it, represented as a polynomial.
346    """
347    # See https://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm.
348    current = [1]
349    prev = [1]
350    b_inv = 1
351    for n, discrepancy in enumerate(syndromes):
352        # Compute discrepancy
353        for i in range(1, len(current)):
354            discrepancy ^= gf.mul(syndromes[n - i], current[i])
355
356        # Correct if discrepancy is nonzero.
357        if discrepancy:
358            x = n + 1 - (len(current) - 1) - (len(prev) - 1)
359            if 2 * (len(current) - 1) <= n:
360                tmp = list(current)
361                current.extend(0 for _ in range(len(prev) + x - len(current)))
362                mul = gf.mul(discrepancy, b_inv)
363                for i, v in enumerate(prev):
364                    current[i + x] ^= gf.mul(mul, v)
365                prev = tmp
366                b_inv = gf.inv(discrepancy)
367            else:
368                mul = gf.mul(discrepancy, b_inv)
369                for i, v in enumerate(prev):
370                    current[i + x] ^= gf.mul(mul, v)
371    return current
372
373class Minisketch:
374    """A Minisketch sketch.
375
376    This represents a sketch of a certain capacity, with elements of a certain bit size.
377    """
378
379    def __init__(self, field_size, capacity):
380        """Initialize an empty sketch with the specified field_size size and capacity."""
381        self.field_size = field_size
382        self.capacity = capacity
383        self.odd_syndromes = [0] * capacity
384        self.gf = GF2Ops(field_size)
385
386    def add(self, element):
387        """Add an element to this sketch. 1 <= element < 2**field_size."""
388        sqr = self.gf.sqr(element)
389        for pos in range(self.capacity):
390            self.odd_syndromes[pos] ^= element
391            element = self.gf.mul(sqr, element)
392
393    def serialized_size(self):
394        """Compute how many bytes a serialization of this sketch will be in size."""
395        return (self.capacity * self.field_size + 7) // 8
396
397    def serialize(self):
398        """Serialize this sketch to bytes."""
399        val = 0
400        for i in range(self.capacity):
401            val |= self.odd_syndromes[i] << (self.field_size * i)
402        return val.to_bytes(self.serialized_size(), 'little')
403
404    def deserialize(self, byte_data):
405        """Deserialize a byte array into this sketch, overwriting its contents."""
406        assert len(byte_data) == self.serialized_size()
407        val = int.from_bytes(byte_data, 'little')
408        for i in range(self.capacity):
409            self.odd_syndromes[i] = (val >> (self.field_size * i)) & ((1 << self.field_size) - 1)
410
411    def clone(self):
412        """Return a clone of this sketch."""
413        ret = Minisketch(self.field_size, self.capacity)
414        ret.odd_syndromes = list(self.odd_syndromes)
415        ret.gf = self.gf
416        return ret
417
418    def merge(self, other):
419        """Merge a sketch with another sketch. Corresponds to XOR'ing their serializations."""
420        assert self.capacity == other.capacity
421        assert self.field_size == other.field_size
422        for i in range(self.capacity):
423            self.odd_syndromes[i] ^= other.odd_syndromes[i]
424
425    def decode(self, max_count=None):
426        """Decode the contents of this sketch.
427
428        Returns either a list of elements or None if undecodable.
429        """
430        # We know the odd syndromes s1=x+y+..., s3=x^3+y^3+..., s5=..., and reconstruct the even
431        # syndromes from this:
432        #  * s2 = x^2+y^2+.... = (x+y+...)^2 = s1^2
433        #  * s4 = x^4+y^4+.... = (x^2+y^2+...)^2 = s2^2
434        #  * s6 = x^6+y^6+.... = (x^3+y^3+...)^2 = s3^2
435        all_syndromes = [0 for _ in range(2 * len(self.odd_syndromes))]
436        for i in range(len(self.odd_syndromes)):
437            all_syndromes[i * 2] = self.odd_syndromes[i]
438            all_syndromes[i * 2 + 1] = self.gf.sqr(all_syndromes[i])
439        # Given the syndromes, find the polynomial that generates them.
440        poly = berlekamp_massey(all_syndromes, self.gf)
441        # Deal with failure and trivial cases.
442        if len(poly) == 0:
443            return None
444        if len(poly) == 1:
445            return []
446        if max_count is not None and len(poly) > 1 + max_count:
447            return None
448        # If the polynomial can be factored into (1-m1*x)*(1-m2*x)*...*(1-mn*x), then {m1,m2,...,mn}
449        # is our set. As each factor (1-m*x) has 1/m as root, we're really just looking for the
450        # inverses of the roots. We find these by reversing the order of the coefficients, and
451        # finding the roots.
452        roots = poly_find_roots(list(reversed(poly)), self.gf)
453        if len(roots) == 0:
454            return None
455        return roots
456
457class TestMinisketch(unittest.TestCase):
458    """Test class for Minisketch."""
459
460    @classmethod
461    def construct_data(cls, field_size, num_a_only, num_b_only, num_both):
462        """Construct two random lists of elements in [1..2**field_size-1].
463
464        Each list will have unique elements that don't appear in the other (num_a_only in the first
465        and num_b_only in the second), and num_both elements will appear in both."""
466        sample = []
467        # Simulate random.sample here (which doesn't work with ranges over 2**63).
468        for _ in range(num_a_only + num_b_only + num_both):
469            while True:
470                r = random.randrange(1, 1 << field_size)
471                if r not in sample:
472                    sample.append(r)
473                    break
474        full_a = sample[:num_a_only + num_both]
475        full_b = sample[num_a_only:]
476        random.shuffle(full_a)
477        random.shuffle(full_b)
478        return full_a, full_b
479
480    def field_size_capacity_test(self, field_size, capacity):
481        """Test Minisketch methods for a specific field and capacity."""
482        used_capacity = random.randrange(capacity + 1)
483        num_a = random.randrange(used_capacity + 1)
484        num_both = random.randrange(min(2 * capacity, (1 << field_size) - 1 - used_capacity) + 1)
485        full_a, full_b = self.construct_data(field_size, num_a, used_capacity - num_a, num_both)
486        sketch_a = Minisketch(field_size, capacity)
487        sketch_b = Minisketch(field_size, capacity)
488        for v in full_a:
489            sketch_a.add(v)
490        for v in full_b:
491            sketch_b.add(v)
492        sketch_combined = sketch_a.clone()
493        sketch_b_ser = sketch_b.serialize()
494        sketch_b_received = Minisketch(field_size, capacity)
495        sketch_b_received.deserialize(sketch_b_ser)
496        sketch_combined.merge(sketch_b_received)
497        decode = sketch_combined.decode()
498        self.assertEqual(decode, sorted(set(full_a) ^ set(full_b)))
499
500    def test(self):