+#include "../include/minisketch.h"
+
+int main(void) {
+
+ minisketch *sketch_a = minisketch_create(12, 0, 4);
+
+ for (int i = 3000; i < 3010; ++i) {
+ minisketch_add_uint64(sketch_a, i);
+ }
+
+ size_t sersize = minisketch_serialized_size(sketch_a);
+ assert(sersize == 12 * 4 / 8); // 4 12-bit values is 6 bytes.
+ unsigned char *buffer_a = malloc(sersize);
+ minisketch_serialize(sketch_a, buffer_a);
+ minisketch_destroy(sketch_a);
+
+ minisketch *sketch_b = minisketch_create(12, 0, 4); // Bob's own sketch
+ for (int i = 3002; i < 3012; ++i) {
+ minisketch_add_uint64(sketch_b, i);
+ }
+
+ sketch_a = minisketch_create(12, 0, 4); // Alice's sketch
+ minisketch_deserialize(sketch_a, buffer_a); // Load Alice's sketch
+ free(buffer_a);
+
+ // Merge the elements from sketch_a into sketch_b. The result is a sketch_b
+ // which contains all elements that occurred in Alice's or Bob's sets, but not
+ // in both.
+ minisketch_merge(sketch_b, sketch_a);
+
+ uint64_t differences[4];
+ ssize_t num_differences = minisketch_decode(sketch_b, 4, differences);
+ minisketch_destroy(sketch_a);
+ minisketch_destroy(sketch_b);
+ if (num_differences < 0) {
+ printf("More than 4 differences!\n");
+ } else {
+ ssize_t i;
+ for (i = 0; i < num_differences; ++i) {
+ printf("%u is in only one of the two sets\n", (unsigned)differences[i]);
+ }
+ }
+}
diff --git a/doc/gen_basefpbits.sage b/doc/gen_basefpbits.sage
new file mode 100644
index 0000000000..d1e75a6e29
--- /dev/null
+++ b/doc/gen_basefpbits.sage
@@ -0,0 +1,78 @@
+# Require exact values up to
+FPBITS = 256
+
+# Overkill accuracy
+F = RealField(400)
+
+def BaseFPBits(bits, capacity):
+ return bits * capacity - int(ceil(F(log(sum(binomial(2**bits - 1, i) for i in range(capacity+1)), 2))))
+
+def Log2Factorial(capacity):
+ return int(floor(log(factorial(capacity), 2)))
+
+print("uint64_t BaseFPBits(uint32_t bits, uint32_t capacity) {")
+print(" // Correction table for low bits/capacities")
+TBLS={}
+FARS={}
+SKIPS={}
+for bits in range(1, 32):
+ TBL = []
+ for capacity in range(1, min(2**bits, FPBITS)):
+ exact = BaseFPBits(bits, capacity)
+ approx = Log2Factorial(capacity)
+ TBL.append((exact, approx))
+ MIN = 10000000000
+ while len(TBL) and ((TBL[-1][0] == TBL[-1][1]) or (TBL[-1][0] >= FPBITS and TBL[-1][1] >= FPBITS)):
+ MIN = min(MIN, TBL[-1][0] - TBL[-1][1])
+ TBL.pop()
+ while len(TBL) and (TBL[-1][0] - TBL[-1][1] == MIN):
+ TBL.pop()
+ SKIP = 0
+ while SKIP < len(TBL) and TBL[SKIP][0] == TBL[SKIP][1]:
+ SKIP += 1
+ DIFFS = [TBL[i][0] - TBL[i][1] for i in range(SKIP, len(TBL))]
+ if len(DIFFS) > 0 and len(DIFFS) * Integer(max(DIFFS)).nbits() > 64:
+ print(" static constexpr uint8_t ADD%i[] = {%s};" % (bits, ", ".join(("%i" % (TBL[i][0] - TBL[i][1])) for i in range(SKIP, len(TBL)))))
+ TBLS[bits] = DIFFS
+ FARS[bits] = MIN
+ SKIPS[bits] = SKIP
+print("")
+print(" if (capacity == 0) return 0;")
+print(" uint64_t ret = 0;")
+print(" if (bits < 32 && capacity >= (1U << bits)) {")
+print(" ret = uint64_t{bits} * (capacity - (1U << bits) + 1);")
+print(" capacity = (1U << bits) - 1;")
+print(" }")
+print(" ret += Log2Factorial(capacity);")
+print(" switch (bits) {")
+for bits in sorted(TBLS.keys()):
+ if len(TBLS[bits]) == 0:
+ continue
+ width = Integer(max(TBLS[bits])).nbits()
+ if len(TBLS[bits]) == 1:
+ add = "%i" % TBLS[bits][0]
+ elif len(TBLS[bits]) * width <= 64:
+ code = sum((2**(width*i) * TBLS[bits][i]) for i in range(len(TBLS[bits])))
+ if width == 1:
+ add = "(0x%x >> (capacity - %i)) & 1" % (code, 1 + SKIPS[bits])
+ else:
+ add = "(0x%x >> %i * (capacity - %i)) & %i" % (code, width, 1 + SKIPS[bits], 2**width - 1)
+ else:
+ add = "ADD%i[capacity - %i]" % (bits, 1 + SKIPS[bits])
+ if len(TBLS[bits]) + SKIPS[bits] == 2**bits - 1:
+ print(" case %i: return ret + (capacity <= %i ? 0 : %s);" % (bits, SKIPS[bits], add))
+ else:
+ print(" case %i: return ret + (capacity <= %i ? 0 : capacity > %i ? %i : %s);" % (bits, SKIPS[bits], len(TBLS[bits]) + SKIPS[bits], FARS[bits], add))
+print(" default: return ret;")
+print(" }")
+print("}")
+
+print("void TestBaseFPBits() {")
+print(" static constexpr uint16_t TBL[20][100] = {%s};" % (", ".join("{" + ", ".join(("%i" % BaseFPBits(bits, capacity)) for capacity in range(0, 100)) + "}" for bits in range(1, 21))))
+print(" for (int bits = 1; bits <= 20; ++bits) {")
+print(" for (int capacity = 0; capacity < 100; ++capacity) {")
+print(" uint64_t computed = BaseFPBits(bits, capacity), exact = TBL[bits - 1][capacity];")
+print(" CHECK(exact == computed || (exact >= 256 && computed >= 256));")
+print(" }")
+print(" }")
+print("}")
diff --git a/doc/gen_params.sage b/doc/gen_params.sage
new file mode 100755
index 0000000000..1cf036adb4
--- /dev/null
+++ b/doc/gen_params.sage
@@ -0,0 +1,333 @@
+#!/usr/bin/env sage
+r"""
+Generate finite field parameters for minisketch.
+
+This script selects the finite fields used by minisketch
+ for various sizes and generates the required tables for
+ the implementation.
+
+The output (after formatting) can be found in src/fields/*.cpp.
+
+"""
+B. = GF(2)
+P. = B[]
+
+def apply_map(m, v):
+ r = 0
+ i = 0
+ while v != 0:
+ if (v & 1):
+ r ^^= m[i]
+ i += 1
+ v >>= 1
+ return r
+
+def recurse_moduli(acc, maxweight, maxdegree):
+ for pos in range(maxweight, maxdegree + 1, 1):
+ poly = acc + p^pos
+ if maxweight == 1:
+ if poly.is_irreducible():
+ return (pos, poly)
+ else:
+ (deg, ret) = recurse_moduli(poly, maxweight - 1, pos - 1)
+ if ret is not None:
+ return (pos, ret)
+ return (None, None)
+
+def compute_moduli(bits):
+ # Return all optimal irreducible polynomials for GF(2^bits)
+ # The result is a list of tuples (weight, degree of second-highest nonzero coefficient, polynomial)
+ maxdegree = bits - 1
+ result = []
+ for weight in range(1, bits, 2):
+ deg, res = None, None
+ while True:
+ ret = recurse_moduli(p^bits + 1, weight, maxdegree)
+ if ret[0] is not None:
+ (deg, res) = ret
+ maxdegree = deg - 1
+ else:
+ break
+ if res is not None:
+ result.append((weight + 2, deg, res))
+ return result
+
+def bits_to_int(vals):
+ ret = 0
+ base = 1
+ for val in vals:
+ ret += Integer(val) * base
+ base *= 2
+ return ret
+
+def sqr_table(f, bits, n=1):
+ ret = []
+ for i in range(bits):
+ ret.append((f^(2^n*i)).integer_representation())
+ return ret
+
+# Compute x**(2**n)
+def pow2(x, n):
+ for i in range(n):
+ x = x**2
+ return x
+
+def qrt_table(F, f, bits):
+ # Table for solving x2 + x = a
+ # This implements the technique from https://www.raco.cat/index.php/PublicacionsMatematiques/article/viewFile/37927/40412, Lemma 1
+ for i in range(bits):
+ if (f**i).trace() != 0:
+ u = f**i
+ ret = []
+ for i in range(0, bits):
+ d = f^i
+ y = sum(pow2(d, j) * sum(pow2(u, k) for k in range(j)) for j in range(1, bits))
+ ret.append(y.integer_representation() ^^ (y.integer_representation() & 1))
+ return ret
+
+def conv_tables(F, NF, bits):
+ # Generate a F(2) linear projection that maps elements from one field
+ # to an isomorphic field with a different modulus.
+ f = F.gen()
+ fp = f.minimal_polynomial()
+ assert(fp == F.modulus())
+ nfp = fp.change_ring(NF)
+ nf = sorted(nfp.roots(multiplicities=False))[0]
+ ret = []
+ matrepr = [[B(0) for x in range(bits)] for y in range(bits)]
+ for i in range(bits):
+ val = (nf**i).integer_representation()
+ ret.append(val)
+ for j in range(bits):
+ matrepr[j][i] = B((val >> j) & 1)
+ mat = Matrix(matrepr).inverse().transpose()
+ ret2 = []
+ for i in range(bits):
+ ret2.append(bits_to_int(mat[i]))
+
+ for t in range(100):
+ f1a = F.random_element()
+ f1b = F.random_element()
+ f1r = f1a * f1b
+ f2a = NF.fetch_int(apply_map(ret, f1a.integer_representation()))
+ f2b = NF.fetch_int(apply_map(ret, f1b.integer_representation()))
+ f2r = NF.fetch_int(apply_map(ret, f1r.integer_representation()))
+ f2s = f2a * f2b
+ assert(f2r == f2s)
+
+ for t in range(100):
+ f2a = NF.random_element()
+ f2b = NF.random_element()
+ f2r = f2a * f2b
+ f1a = F.fetch_int(apply_map(ret2, f2a.integer_representation()))
+ f1b = F.fetch_int(apply_map(ret2, f2b.integer_representation()))
+ f1r = F.fetch_int(apply_map(ret2, f2r.integer_representation()))
+ f1s = f1a * f1b
+ assert(f1r == f1s)
+
+ return (ret, ret2)
+
+def fmt(i,typ):
+ if i == 0:
+ return "0"
+ else:
+ return "0x%x" % i
+
+def lintranstype(typ, bits, maxtbl):
+ gsize = min(maxtbl, bits)
+ array_size = (bits + gsize - 1) // gsize
+ bits_list = []
+ total = 0
+ for i in range(array_size):
+ rsize = (bits - total + array_size - i - 1) // (array_size - i)
+ total += rsize
+ bits_list.append(rsize)
+ return "RecLinTrans<%s, %s>" % (typ, ", ".join("%i" % x for x in bits_list))
+
+INT=0
+CLMUL=1
+CLMUL_TRI=2
+MD=3
+
+def print_modulus_md(mod):
+ ret = ""
+ pos = mod.degree()
+ for c in reversed(list(mod)):
+ if c:
+ if ret:
+ ret += " + "
+ if pos == 0:
+ ret += "1"
+ elif pos == 1:
+ ret += "x"
+ else:
+ ret += "x%i" % pos
+ pos -= 1
+ return ret
+
+def pick_modulus(bits, style):
+ # Choose the lexicographicly-first lowest-weight modulus
+ # optionally subject to implementation specific constraints.
+ moduli = compute_moduli(bits)
+ if style == INT or style == MD:
+ multi_sqr = False
+ need_trans = False
+ elif style == CLMUL:
+ # Fast CLMUL reduction requires that bits + the highest
+ # set bit are less than 66.
+ moduli = list(filter((lambda x: bits+x[1] <= 66), moduli)) + moduli
+ multi_sqr = True
+ need_trans = True
+ if not moduli or moduli[0][2].change_ring(ZZ)(2) == 3 + 2**bits:
+ # For modulus 3, CLMUL_TRI is obviously better.
+ return None
+ elif style == CLMUL_TRI:
+ moduli = list(filter(lambda x: bits+x[1] <= 66, moduli)) + moduli
+ moduli = list(filter(lambda x: x[0] == 3, moduli))
+ multi_sqr = True
+ need_trans = True
+ else:
+ assert(False)
+ if not moduli:
+ return None
+ return moduli[0][2]
+
+def print_result(bits, style):
+ if style == INT:
+ multi_sqr = False
+ need_trans = False
+ table_id = "%i" % bits
+ elif style == MD:
+ pass
+ elif style == CLMUL:
+ multi_sqr = True
+ need_trans = True
+ table_id = "%i" % bits
+ elif style == CLMUL_TRI:
+ multi_sqr = True
+ need_trans = True
+ table_id = "TRI%i" % bits
+ else:
+ assert(False)
+
+ nmodulus = pick_modulus(bits, INT)
+ modulus = pick_modulus(bits, style)
+ if modulus is None:
+ return
+
+ if style == MD:
+ print("* *%s*" % print_modulus_md(modulus))
+ return
+
+ if bits > 32:
+ typ = "uint64_t"
+ elif bits > 16:
+ typ = "uint32_t"
+ elif bits > 8:
+ typ = "uint16_t"
+ else:
+ typ = "uint8_t"
+
+ ttyp = lintranstype(typ, bits, 4)
+ rtyp = lintranstype(typ, bits, 6)
+
+ F. = GF(2**bits, modulus=modulus)
+
+ include_table = True
+ if style != INT and style != CLMUL:
+ cmodulus = pick_modulus(bits, CLMUL)
+ if cmodulus == modulus:
+ include_table = False
+ table_id = "%i" % bits
+
+ if include_table:
+ print("typedef %s StatTable%s;" % (rtyp, table_id))
+ rtyp = "StatTable%s" % table_id
+ if (style == INT):
+ print("typedef %s DynTable%s;" % (ttyp, table_id))
+ ttyp = "DynTable%s" % table_id
+
+ if need_trans:
+ if modulus != nmodulus:
+ # If the bitstream modulus is not the best modulus for
+ # this implementation a conversion table will be needed.
+ ctyp = rtyp
+ NF. = GF(2**bits, modulus=nmodulus)
+ ctables = conv_tables(NF, F, bits)
+ loadtbl = "&LOAD_TABLE_%s" % table_id
+ savetbl = "&SAVE_TABLE_%s" % table_id
+ if include_table:
+ print("constexpr %s LOAD_TABLE_%s({%s});" % (ctyp, table_id, ", ".join([fmt(x,typ) for x in ctables[0]])))
+ print("constexpr %s SAVE_TABLE_%s({%s});" % (ctyp, table_id, ", ".join([fmt(x,typ) for x in ctables[1]])))
+ else:
+ ctyp = "IdTrans"
+ loadtbl = "&ID_TRANS"
+ savetbl = "&ID_TRANS"
+ else:
+ assert(modulus == nmodulus)
+
+ if include_table:
+ print("constexpr %s SQR_TABLE_%s({%s});" % (rtyp, table_id, ", ".join([fmt(x,typ) for x in sqr_table(f, bits, 1)])))
+ if multi_sqr:
+ # Repeated squaring is a linearised polynomial so in F(2^n) it is
+ # F(2) linear and can be computed by a simple bit-matrix.
+ # Repeated squaring is especially useful in powering ladders such as
+ # for inversion.
+ # When certain repeated squaring tables are not in use, use the QRT
+ # table instead to make the C++ compiler happy (it always has the
+ # same type).
+ sqr2 = "&QRT_TABLE_%s" % table_id
+ sqr4 = "&QRT_TABLE_%s" % table_id
+ sqr8 = "&QRT_TABLE_%s" % table_id
+ sqr16 = "&QRT_TABLE_%s" % table_id
+ if ((bits - 1) >= 4):
+ if include_table:
+ print("constexpr %s SQR2_TABLE_%s({%s});" % (rtyp, table_id, ", ".join([fmt(x,typ) for x in sqr_table(f, bits, 2)])))
+ sqr2 = "&SQR2_TABLE_%s" % table_id
+ if ((bits - 1) >= 8):
+ if include_table:
+ print("constexpr %s SQR4_TABLE_%s({%s});" % (rtyp, table_id, ", ".join([fmt(x,typ) for x in sqr_table(f, bits, 4)])))
+ sqr4 = "&SQR4_TABLE_%s" % table_id
+ if ((bits - 1) >= 16):
+ if include_table:
+ print("constexpr %s SQR8_TABLE_%s({%s});" % (rtyp, table_id, ", ".join([fmt(x,typ) for x in sqr_table(f, bits, 8)])))
+ sqr8 = "&SQR8_TABLE_%s" % table_id
+ if ((bits - 1) >= 32):
+ if include_table:
+ print("constexpr %s SQR16_TABLE_%s({%s});" % (rtyp, table_id, ", ".join([fmt(x,typ) for x in sqr_table(f, bits, 16)])))
+ sqr16 = "&SQR16_TABLE_%s" % table_id
+ if include_table:
+ print("constexpr %s QRT_TABLE_%s({%s});" % (rtyp, table_id, ", ".join([fmt(x,typ) for x in qrt_table(F, f, bits)])))
+
+ modulus_weight = modulus.hamming_weight()
+ modulus_degree = (modulus - p**bits).degree()
+ modulus_int = (modulus - p**bits).change_ring(ZZ)(2)
+
+ lfsr = ""
+
+ if style == INT:
+ print("typedef Field<%s, %i, %i, %s, %s, &SQR_TABLE_%s, &QRT_TABLE_%s%s> Field%i;" % (typ, bits, modulus_int, rtyp, ttyp, table_id, table_id, lfsr, bits))
+ elif style == CLMUL:
+ print("typedef Field<%s, %i, %i, %s, &SQR_TABLE_%s, %s, %s, %s, %s, &QRT_TABLE_%s, %s, %s, %s%s> Field%i;" % (typ, bits, modulus_int, rtyp, table_id, sqr2, sqr4, sqr8, sqr16, table_id, ctyp, loadtbl, savetbl, lfsr, bits))
+ elif style == CLMUL_TRI:
+ print("typedef FieldTri<%s, %i, %i, %s, &SQR_TABLE_%s, %s, %s, %s, %s, &QRT_TABLE_%s, %s, %s, %s> FieldTri%i;" % (typ, bits, modulus_degree, rtyp, table_id, sqr2, sqr4, sqr8, sqr16, table_id, ctyp, loadtbl, savetbl, bits))
+ else:
+ assert(False)
+
+for bits in range(2, 65):
+ print("#ifdef ENABLE_FIELD_INT_%i" % bits)
+ print("// %i bit field" % bits)
+ print_result(bits, INT)
+ print("#endif")
+ print("")
+
+for bits in range(2, 65):
+ print("#ifdef ENABLE_FIELD_INT_%i" % bits)
+ print("// %i bit field" % bits)
+ print_result(bits, CLMUL)
+ print_result(bits, CLMUL_TRI)
+ print("#endif")
+ print("")
+
+for bits in range(2, 65):
+ print_result(bits, MD)
diff --git a/doc/log2_factorial.sage b/doc/log2_factorial.sage
new file mode 100644
index 0000000000..afc6d66c57
--- /dev/null
+++ b/doc/log2_factorial.sage
@@ -0,0 +1,85 @@
+import bisect
+
+INPUT_BITS = 32
+TABLE_BITS = 5
+INT_BITS = 64
+EXACT_FPBITS = 256
+
+F = RealField(100) # overkill
+
+def BestOverApproxInvLog2(mulof, maxd):
+ """
+ Compute denominator of an approximation of 1/log(2).
+
+ Specifically, find the value of d (<= maxd, and a multiple of mulof)
+ such that ceil(d/log(2))/d is the best approximation of 1/log(2).
+ """
+ dist=1
+ best=0
+ # Precomputed denominators that lead to good approximations of 1/log(2)
+ for d in [1, 2, 9, 70, 131, 192, 445, 1588, 4319, 11369, 18419, 25469, 287209, 836158, 3057423, 8336111, 21950910, 35565709, 49180508, 161156323, 273132138, 385107953, 882191721]:
+ kd = lcm(mulof, d)
+ if kd <= maxd:
+ n = ceil(kd / log(2))
+ dis = F((n / kd) - 1 / log(2))
+ if dis < dist:
+ dist = dis
+ best = kd
+ return best
+
+
+LOG2_TABLE = []
+A = 0
+B = 0
+C = 0
+D = 0
+K = 0
+
+def Setup(k):
+ global LOG2_TABLE, A, B, C, D, K
+ K = k
+ LOG2_TABLE = []
+ for i in range(2 ** TABLE_BITS):
+ LOG2_TABLE.append(int(floor(F(K * log(1 + i / 2**TABLE_BITS, 2)))))
+
+ # Maximum for (2*x+1)*LogK2(x)
+ max_T = (2^(INPUT_BITS + 1) - 1) * (INPUT_BITS*K - 1)
+ # Maximum for A
+ max_A = (2^INT_BITS - 1) // max_T
+ D = BestOverApproxInvLog2(2 * K, max_A * 2 * K)
+ A = D // (2 * K)
+ B = int(ceil(F(D/log(2))))
+ C = int(floor(F(D*log(2*pi,2)/2)))
+
+def LogK2(n):
+ assert(n >= 1 and n < (1 << INPUT_BITS))
+ bits = Integer(n).nbits()
+ return K * (bits - 1) + LOG2_TABLE[((n << (INPUT_BITS - bits)) >> (INPUT_BITS - TABLE_BITS - 1)) - 2**TABLE_BITS]
+
+def Log2Fact(n):
+ # Use formula (A*(2*x+1)*LogK2(x) - B*x + C) / D
+ return (A*(2*n+1)*LogK2(n) - B*n + C) // D + (n < 3)
+
+RES = [int(F(log(factorial(i),2))) for i in range(EXACT_FPBITS * 10)]
+
+best_worst_ratio = 0
+
+for K in range(1, 10000):
+ Setup(K)
+ assert(LogK2(1) == 0)
+ assert(LogK2(2) == K)
+ assert(LogK2(4) == 2 * K)
+ good = True
+ worst_ratio = 1
+ for i in range(1, EXACT_FPBITS * 10):
+ exact = RES[i]
+ approx = Log2Fact(i)
+ if not (approx <= exact and ((approx == exact) or (approx >= EXACT_FPBITS and exact >= EXACT_FPBITS))):
+ good = False
+ break
+ if worst_ratio * exact > approx:
+ worst_ratio = approx / exact
+ if good and worst_ratio > best_worst_ratio:
+ best_worst_ratio = worst_ratio
+ print("Formula: (%i*(2*x+1)*floor(%i*log2(x)) - %i*x + %i) / %i; log(max_ratio)=%f" % (A, K, B, C, D, RR(-log(worst_ratio))))
+ print("LOG2K_TABLE: %r" % LOG2_TABLE)
diff --git a/doc/math.md b/doc/math.md
new file mode 100644
index 0000000000..cf46f193ab
--- /dev/null
+++ b/doc/math.md
@@ -0,0 +1,117 @@
+# The mathematics of Minisketch sketches
+
+This is an unconventional mathematical overview of the PinSketch algorithm without references to coding theory[[1]](#myfootnote1).
+
+## Set sketches
+
+A sketch, for the purpose of this description, can be seen as a "set checksum" with two peculiar properties:
+
+* Sketches have a predetermined capacity, and when the number of elements in the set is not higher than the capacity, minisketch will always recover the entire set from the sketch. A sketch of *b*-bit elements with capacity *c* can be stored in *bc* bits.
+* The sketches of two sets can be combined by adding them (XOR) to obtain a sketch of the [symmetric difference](https://en.wikipedia.org/wiki/Symmetric_difference) between the two sets (*i.e.*, all elements that occur in one but not both input sets).
+
+This overview explains how sets can be converted into a sketch and how a set can be recovered from a sketch.
+
+## From field elements to sketches
+
+**Data entries as field elements**
+
+Every integer in the range *[1...2b-1]* (the acceptable data elements for a Minisketch sketch with field size *b*) can be mapped to a nonzero field element of *GF(2b)*. In this [finite field](https://en.wikipedia.org/wiki/Finite_field), we can add and multiply elements together, with many of the expected properties for those operations. Addition (and subtraction!) of field elements corresponds to bitwise XOR of the integers they correspond to, though multiplication is more involved.
+
+**Sets as power series**
+
+We define a function *S* which maps field elements *m* to the following [formal power series](https://en.wikipedia.org/wiki/Formal_power_series) (similar to a polynomial, except there can be an infinite number of terms, and we don't care about concepts like convergence as we're never going to actually evaluate it for a specific value of *x*):
+
+* *S(m) = 1 + mx + m2x2 + m3x3 + ...*.
+
+We then extend this function to operate on sets of field elements, by adding together the images of every set element. If *M = {m1, m2, ... }*:
+
+* *S(M) = S({m1,m2,...}) = S(m1) + S(m2) + ... = (1 + 1 + ...) + (m1 + m2 + ...)x + (m12 + m22 + ...)x2 + (m13 + ...*
+
+Because in our field addition corresponds to XOR of integers, it holds for every *a* that *a + a = 0*. This carries over to the *S* function, meaning that *S(a) + S(a) = 0* for every *a*. This means that the coefficients of these power series have the second of the properties we
+desire from a sketch, namely that an efficient operation exists to
+combine two sketches such that the result is a sketch of the symmetric
+difference of the sets. It holds that
+*S({m1,m2}) + S({m2,m3}) = S(m1) + (S(m2) + S(m2)) + S(m3) = S(m1) + S(m3) = S({m1,m3})*. The question is whether we can also efficiently recover the elements from their power series' coefficients.
+
+**An infinity of coefficients is hard**
+
+To make reasoning about these power series easier, notice that the series for a single element is in fact a [geometric series](https://en.wikipedia.org/wiki/Geometric_series). If we were working over real numbers rather than a finite field and *|mx| < 1*, it would converge to *(1 - mx)-1*. Convergence has no meaning in formal power series, however it is still the case that:
+
+* *(1 - mx) S(m) = 1*
+
+You can verify this by seeing that every coefficient except the constant one gets cancelled out by the multiplication. This can be generalized to the series for multiple set elements. For two elements we have:
+
+* *(1 - m1x) (1 - m2x) S({m1,m2}) = (1 - m1x) (1 - m2x) (S(m1) + S(m2)) = (1 - m2x) + (1 - m1x)*
+
+And for three:
+
+* *(1 - m1x) (1 - m2x) (1 - m3x) S({m1,m2,m3}) = (1 - m1x) (1 - m2x) (1 - m3x) (S(m1) + S(m2) + S(m3)) = (1 - m2x)(1 - m3x) + (1 - m1x)(1 - m3x) + (1 - m1x)(1 - m2x)*
+
+In each case, we notice that multiplying *S(M)* with *(1 - mix)* for each element *mi ∈ M* results in a polynomial of degree *n-1*.
+
+**Solving for the set elements**
+
+The above insight lets us build a solver that extracts the set elements from the coefficients of a power series. If we can find a polynomial *L* that is the product of *n* different *(1 - mix)* factors for various values of *mi*, such that *P = S(M)L* is an *n-1* degree polynomial, then those values *mi* are the elements of *M*.
+
+The coefficients of *P* are nontrivial expressions of the set elements themselves. However, we can just focus on the coefficients of degree *n* and higher in *P*, as those are all 0. Let *si* be the coefficients of *S(M)*, and *li* the coefficients of L. In other words, *S(M) = s0 + s1x + s2x2 + s3x3 + ...* and *L = l0 + l1x + l2x2 + l3x3 + ... + lnxn*. Note that *l0 = 1*, as it is the product of all the *1* terms in the *(1 - mix)* factors.
+
+Here are the equations for the coefficients of *S(M)L* of degree *n+1* through *2n*:
+* *sn+1 + sn+0l1 + sn-1l2 + sn-2l3 + ... + s1ln = 0*
+* *sn+2 + sn+1l1 + sn+0l2 + sn-1l3 + ... + s2ln = 0*
+* *sn+3 + sn+2l1 + sn+1l2 + sn+0l3 + ... + s3ln = 0*
+* ...
+* *s2n + s2n-1l1 + s2n-2l2 + s2n-3l3 + ... + snln = 0*
+
+These are *n* linear equations with *n* unknowns (the *li*
+values, for *i=1..n*), which can be solved using [Gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination). After doing so,
+we have the coefficients of *L*, which can then be [factored](https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields)
+into first degree factors of the form *(1 - mix)*. The resulting *m* values are our set elements.
+
+**Putting it all together**
+
+Interestingly, only *2n* coefficients of *S(M)* were needed for solving
+the set of equations above. This means we have our answer: the
+coefficients *1* through *2n* of *S(M)*, or the list
+*[m1 + m2 + ..., m12 + m22 + ..., ..., m12n + m22n + ...]*
+functions as a sketch, satisfying the two properties we want:
+
+* Sketches can be combined to form the sketch of their symmetric difference, by simply pairwise adding the list elements together.
+* With *2n* list elements we can efficiently recover *n* elements from a sketch.
+
+**Capacity and difference**
+
+The approach above only works when the number of elements *n* in the sketch is known. Of course we want to support cases where only an upper bound on the number of elements in the sketch is known, the capacity *c*. Given that we can reconstruct a set of size *c* from a sketch with *2c* terms, we should be able to reconstruct a set of size *n* too as long as *n ≤ c*. This is simply a matter of trying to solve the above set of equations assuming values of *n* that count down from *c* until a solution is found for one. This is known as the [Peterson-Gorenstein-Zierler algorithm](https://en.wikipedia.org/wiki/BCH_code#Peterson%E2%80%93Gorenstein%E2%80%93Zierler_algorithm).
+
+## Optimizations
+
+**Halving the sketch size**
+
+We can in fact only include the odd terms in the sketch, and reconstruct the even ones before solving the equation to find *L*. This means the size of a sketch becomes just *c* field elements, the same size as would be needed to send its contents naively.
+
+To see how this is possible, we need the [Frobenius endomorphism](https://en.wikipedia.org/wiki/Frobenius_endomorphism), which in short states that in fields where *x + x = 0* it holds that *(x + y)2 = x2 + y2* for every *x* and *y* (the dream of every high school math student!). This means that:
+
+* *s2 = m12 + m22 + ... = (m1 + m2 + ...)2 = s12*.
+* *s4 = m14 + m24 + ... = (m12 + m22 + ...)2 = s22*.
+* *s6 = m16 + m26 + ... = (m13 + m23 + ...)2 = s32*.
+* ...
+
+In other words, we only need to send *s1, s3, s5, ..., s2n-1* to recover all *2n* *si* values, and proceed with reconstruction.
+
+**Quadratic performance rather than cubic**
+
+Using Gaussian elimination to solve the set of equations above for the *li* values requires *O(n3)* field operations. However, due to the special structure in the equations (look at the repeated *si* values), it can be solved in *O(n2)* time using a number of techniques, including the [Berlekamp-Massey algorithm](https://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm) (BM).
+
+**Roots instead of factorization**
+
+As explained above, the polynomial *L* can be factored into *(1 - mix)* factors, where the values *mi* are the set elements. However, since we know that a decodable sketch must result in a polynomial that is fully factorizable into degree-*1* factors, we can instead use a more efficient root-finding algorithm rather than a factorization algorithm. As the root of each *(1 - mix)* factor is *mi-1*, we conclude that the set elements are in fact the inverses of the roots of *L*.
+
+**Avoiding inversions**
+
+As inversions are a relatively expensive operation, it would be useful to avoid them.
+
+Say that we're trying to find the inverses of the roots of *L = 1 + l1x + l2x2 + ... + lnxn*, then we're really interested in the solutions *y* for *1 + l1y-1 + l2y-2 + ... + lny-n = 0*. By multiplying both sides in the equations with *yn*, we find *ln + ln-1y + ln-2y2 + ... + yn = 0*.
+
+In other words, we can find the inverses of the roots of *L* by instead factoring the polynomial with the coefficients of *L* in reverse order.
+
+* [1] For those familiar with coding theory: PinSketch communicates a set difference by encoding the set members as errors in a binary [BCH](https://en.wikipedia.org/wiki/BCH_code) codeword 2bits in size and sends the syndromes.
+ The linearity of the syndromes provides all the properties needed for a sketch. Sketch decoding is simply finding the error locations. Decode is much faster than an ordinary BCH decoder for such a large codeword because the need to take a discrete log is avoided by storing the set in the roots directly instead of in an exponent (logically permuting the bits of the codeword).
diff --git a/doc/minisketch-vs.png b/doc/minisketch-vs.png
new file mode 100644
index 0000000000000000000000000000000000000000..aed810de8a2512a2e5bfb33fc9e9a99c150b98b0
GIT binary patch
literal 14212
zcmZ{L19WB05^g-PZB1|{wr$(CCOYxNP9~n%wr$&-Wa3P0TQBop{qMf}&RM;?c2#{<
z-PKjSd#!y=n1Y-*0xT{p2nYy*q=bmlN1yo7yr988u2Y)`Hy<4+P)S?}q;eAf@Z+Y}
zR728CRu%;CK|_Ol0!0M@|D*Ep1p&ncf&7aG0Wku_`!{U`O7%~j51XKjAb%{Qd^AT;
zroU+P4>~Z!5(MI(*yNA?4@vyz@~>XO(#_Q7tD%#fBQY}_Cmj_1$&lFXtBoy?%Z-=xFAc5_`VX6)
zl=v?dCu?3(4Os1E+rAMf5JcRcu6gsob0*i>0Mo2
z>0FuV>>SPM896yQ=^2>lnV4ulG-!eDwoZm_w6;L9e;E0f9T8KYv7@EElck+4@gKW}
zMt06lyriUm9R2(Ak8zq9|J#qfv!l&lUQCSXO>IoSn%X)6=^5!5>Hin*W48YW+5-O>
zjgPp{yBXSlfEfNK$kfg9KhyfZP=6%I
zf8&10{)(}Xp*j8E*gW+AQ=8{weq6G4CYEOIB8E<;e1DRPmVt$qiBaWWNyS6|UpoJZ
z?q6tON7FxPsbXhm!}s4=`VZ{iHU0+v6AP}t5E>$uAFBZTSVlfZRt6sW|0(&ONX@??
znOOgh{F~&zk-~O1c8<#ShQ_9RjQ@oEjrBL}KVjz*w)<-DXbSvmbRT2;e|Z0`uc2fL
zw6k&klazc+931~d{GI#X2zf`#j~q1o3-K|azdZb%_uu;eaiD5z>BPtMPso3<{~M+G
zANk+L^>^++5j^yNYUyva^^d~+i~CWze6Sy@`mbGp57w+i%jIJ;qLdU7RB;16&2X7aLY|GqwF`*od>EzMBl6E?J=v|x`(GOmHpTQY93f*~!j03-lZ2w|K^Mz$s&
zN{QM4S{WHkM#Fl^AGB^6o4(dWyOx%><8{q7X}#^5(<{HMh*Env&aF}U{y_Ek=J0IJ
zYpTiT`8=Jiv^1_&fe~JKQG7*w26~Cdx@z%xSe(ysXE;`xu7!N8wwA{AWGS}raA(H^
ze@GD7-hbac?Z_zAEOTF{p}ARNdwzcY?sP@F;^)J|!-e_O_&5y>&GaNjJB|*a=uCfq
zU}s`kuz3_9pUJoY*?A34Q9r%#Gxm0nj51mNe
z6J~X5w!m6jCta+3{+BHzeUHNef!msBuahDEH%7bey@@o=ijDKR{i7qm8W_o`x~?-u
zq!TwAo0;a%lM@=dk=-PHxvA{8QT5I^G2gD&yOr~neZtq)=`P(QvAi17G0e&Rf*fsI
zab9@gOuub=lrxFl8;3Bpq^OB;4ZUjXsm{4*zEbV^svrx(_lKdb*Hb(GKKhyY5p2fh
zsBI7^O6Jb`fYG*%S6hR^-F^P|$FEImcOT(wilCtsE@?*q0)Zz&bM+ulc2F57D@{l`
zBx-ull?9V!-v{}m>grfF4Cz1Gy`Rv(uXVZN-^26Y;UkJV^jBH1
zrS99Gy6fczWtf^y=4YEWHHlxch&pceO*PrlGr5QtzpW2H3w;s
z;Nl>>)hv1c&VI^mklD85_9s@yf#~)1wE>jZvA4lM*V}DX&-ukxp8`)iR{GMy1Ic@~
zbALWx3fg8z>s#LW`^8<`IRb8^lvH0eZlbBSoN6NbOV|1PHNXCCkshSq7n97bko7<-pJ!vzzU(#Rx!*YxuhyAItLft?2fv+E>AF(1!n9&U68EoU
zmt4$cz1`@#t`#iEX$d6Jw#5)Q!_rUFHmzRuAbH1(wP2n?B^*uX2+wCdzhFf)430_R
zRyrTk5Sx|)ksuJoEiG{o8N$u1ADdAZN7HzKt_vF-J;Sv4Cf)RSTq-qw?CR
z&qK7`)6@H_))+=voWmm1xZ8WhbJY@CVrtLLPuky0m#b@KO+H`RZ+^KLk6lO-ief;L
zs0_zsAYlmcYGKaeVOxafzco~Ooo@C7fZTnQ3~cA8RxPxd#tgp4jB`&)EqG^GY!HSE
zU)Xxl92E8L=Z*I_B+pH7e4`|BNLql6V3aQW%h#30jzd{HBWZ4ah}nI;dFqu$yLb6>
z?|XCpJAJwhd^J%f)^*!pOM?Wz*at}SVCL^@H1BJ706-Jp`q~wT*Fv5}8FAp-QtmYW
zCa^E#Pnk0v$h(g!XWryf`rRKEG5$fp|bM{i74>{&k0vKTKo?)s2j0mp9VveOBP>bHlXUn&A#3l?k!@u^j78aDzZb?<;xmp3Y9fx)%3~@c
zjO7&+n-9Q3fdq~&hlqXY>AF3&^Lg%kKI?l9YIyqy*Zz0?m*2pw=hc0lea_>HU?L%N
zoz8Ffwp1WH(~f@O%roSIVCaUv$OsUlk=@ZhN`RqwUT46pwzsEaPQNir2=WI?28H10
z{8#!Lknk9$^G^c33|V*eMKc`voRdXBZ^Xe7W1?pbFZVXyw|1<~lTZ5ZLqU4YW0Bfr
z1Yefa%-p#`0piPopJVYpiCQa-ANbN>GP!qXAf=V?y^NS^Dk76QW+Tx`+?-)`UIcxt
zPn#RrN6{6%-|D}!3X5@HokVe;bz6z`3Vx^eI+V_l=5>ik{Yj3ZK2IVt`6*7l>uJHF
z)CrOy8`=%sRs_^CIFnX4hqsr!g}^D`8+;50s1r)a@`yjpOIz17_XB9iL85eXVw*iK
zSb`j}AP@Jq`?K9nbUlvdrpn_}X1%LYJfEA1ZbqBB+HG4^8-nWyKb69-DFNOr5?8sf
zc!Dd8WkK5T)FwxVMlXcX3(IdrgRTC^vVMNkS6MkdeUGo_KJOhLHFcf!{^BejtdR@y
zXwLtv=~-&u743bK*A4aSPy#Z3hzgw2X{4ds7aQDCA!OUIW9vPTgV^W+U~?g@|1^u_
zeF07-Y0X%FA$LJ;lOvH#f_&vh^fVoEfY$BCe~=YoQnrBs4`UmiMws%%?sLta%)Xk$
zyvX}_m(iB>mdoF^t?y**Fb;eMKfPjZDiVOx#dEeDp_Dl8P
zO|A$aAv0!N$83u4+ZZ<1&w;J$zg;D5+%+~iz)av!C!4%t3iHt~1i
z+9GWBF+gq5wvMHsaCM$)dgHjN?R_69%%Z2`xXY;JqU3%PvDpptBFBOxAscSqoO!tB
ze+Rl4TG+XdX?LF|b-jGX#}K@|O6tOw0&AD|I#Tu(e|4_KydJ~T8FZU9u)Cj*z_atjU}d4*z~w|6DKhjk=(XFnBd0^XLs*^-haSW7
z;1s_0Q}V@1XNy8Z?;L#6+KFQEUscX}QqdzgO5J2JyD5L(LM{Gb?NE92490mUlStt(
zZ;r>3+VKWp@g6|SY4ONSrR5j;ET83i!%xcFDCp=);uI1OoeZ~t$vf}i~$6W6Kr76k5t-yc+foE
zbEc#)Fv2#m5e&jjoss(3TU<
zzS>dY@C#Wi_rIp*0i1stXgrw(;E%B(lnpbDxm3PxC%tY?S2&)7HoDrg+#VNtluI-KmeDL>m7Jt{$awr*_H
zFQv$j_!64XV8txQ0gR-wTb0lE=lh-Y!ji=qxI$uCn6{un2$2vjip_lQk5k5{4!E_6
z?7nHemW-6dnonDL!iZC2B3aYvWM;9u;d#FscUe*MhQr6+4@B7~2&wV}&2UkwVs@rA
zXSE%8He7PorK%ViGo6z1gqwm^?r|rAMR+CE^IJvzkNe&FM?@h>=-;@MU!7tmSHa9K
z#O#eswnDf(20klbhz-ztUw+?sJ#)p?qG3CMsK@IbWPJT|K;KKVGdo}$yrmV86ApF#
z_5jF!Yv1QR7tG56R{oY*hiY7n^1d1KqZX2ZeGYWF$m=kJ2~u~n$jA8|$B*FzSi8Hq
z>7=|ldkI@{0@R5-puoDs%&oos)(sh?YnvY@)8ODz^0+y8s62fkKcxLWhvnL}q#nvs
zYoRDO%nb(vMt%$Y(rnj^pMXJ*ttMf$gcK0v}#PK!W
zQ0`2-609B;sj*5M?&RcRt~IONn$p~eQ(9Bhsra7Gv1ygh4*483#<79wTQEP7?j)7~
zi3KIdolcWT5JUkpmuGX}LI)(5I}HyGoT<3R3{g`ylsG!S++b)m#ow)i*Y)RBwPWr%8K{#ZyO@7iDW
z-h}LudNqJOf9A59J~eL<6$TZW_nHhTMEPd!zI9mC19gk%>-H0hXn%9Vy$^34M0d6~
z4-Cmyu>3^mD+tv2E|o3BLG-oden%jYLWp!~l(hz(fJf!=5@sr@drje@#R+K}VuK3j
zC;#eQh|CQuH+|wrC#FxPOXxPQ!8D4uKoIfvS$bF^
z_YLfp!!E+6mcow0XNY!KhW=Cs0xS~CB#cqn=IdprOZ#7K?2J6zzIY{6QfioymxT!h
z*yI=BbrO<3(rs&QpFcY<(DX0uh0eq*L}52clA%~5<$K%m-}&n};|qB+n>bDKGAS{z
zN8)dXw%w+;`RCV<;W?T>Eb?KTK24vWrY}MlhC`JJsOFm_V@6;VS|!3LDJz;KPy{cn
zPY&&daEVF=#|uW~Di`R4{veGu9;2xMMCXCh{Iq5>7BNdt*{$6VVLA;%lj4utk%Rd(
z>o(pO`r83I>mrbHQ#2rs;oLxW8N7F8f?ICOG(^xa#t*@cJwO^d3@Z^UohSkm+`KPy
zt_Cp}NjZb6=59|wsh;^}>3JIAB$c+?m#gs;1){&V7Vb6-bDz5UY*W;j^eorCs+p^WsSQWyy*TVgJ?s*>5=
z!(Dqvn+6|;1Zj-**VINumr62K?Ssb
z!s7=MGk`EBG1o@03d0LBbb^#SBhd9g=X{B}lAE5apw153PFr^ayyOZ^tLmnZQO4c
zruhEU=p(?hTe&2TCZ%_h?7J$3$53nUz7g7B%1yv1<6U3w)Xt5Hki|}OrFMkZxorCG
zS`b%`rza=Wyw&9Gi{~G%MlH`jNV3jm|9d(hqmQH-LU)~0FQXwN9{ctn#}ZYP5z`p=
z^V))3x`<^Wx+USn>goYOuj9o;Derow(;h~=XEofszLRxyG}d>ACBFNL;iF4C3|RtB_+8o`z0_|;p5|rZK~cW
zm*m%-=1x>+Hj9gBU9og;%DZH_12x8T385u+GQobT*1CZsPkJ#Nf6(hU4NI*}2y8pZ
zxS>};B|vJK3G&EQdv@a;;@BqU{MJ^Vjm&N{&1!NQZ-lGjjv0GCBonAFrUiTimwZ{b{c(t?CP?j0<+dUDBcV0z
z#jH?+d;5`A62EK9$zTwmjoHI**pDyc8JGmsr$9t7_En$?v=u$cyzoD>+w?YW#9;aB|fek1lZ6JYWfWu;a31_Qs(3
zl2|Avt^rV}y5qJDt*mME{5jSnf^kJ?yN8lodTf?RqMd4nR*6ua#I2hsl=-(OPpb8E
zdor=gGID;nU!h^1gUV4_(BKpJf#;hgL3+OUGt4?b4|
zR}&N3k%Gg;!az$r3A%pq+3#KpJLy@n-dn(?x@c$$Thyn!jIjFhcWh8!Rx~g;1>=D2
z!HBqYm0G!q>e1ZjV&$}Q4tFYPxR5da$$uHg&3GHL65q`}sEUr(ej4GaWAbRJzj+S2
zZ9GSt2fsk_+lm}}1$-rNt7%+(ATLwWLm_?5C_OHT#WNmGsRf!O5S^B%pjz;^x)q=C
zUZEOTX}H)|`pc%GnQj35BEuLB?bURDep3h2$go+x*P?_%z_?hRq|~P~JJ!Sbva2BD
z8cWW$dz=Y{#~z2|eqN30WZp})0{0e;Hacj7B3DO6A&qP0+$*b#3YI2(n_yUvaV1j4
z8c~%d0`w02+d7G9gWB>^R5jt@*aD3ZZ)BnB7!!Y&c+5_7xq4g(VfzEznNk~EN)CUi
zvIs8$%CbHda;{-`amIQ0pKy9~a?pyVF}Cn6hFB$OZQ1+eh=;xQg7B~X8>#ytYqPY5
z8PE#9$Cr={P}y{3zL!D^?fRJGl7Z4XX{X}DEez9fU5BSq7WNM@f8U;04}{JYiKyTn
zkRIrFC>~FP)1aep-s`Vk(Bl&Ilj;ymaZnJKz51oV_#i)!G%lt-*l$W08m-Lnxd*(Z
zu{Koi2%>TeAWmf#m=_T{3d`f_sXWpyY7}FU-vOD@AumDb8%qiY^CLfOR&l!$9#Kjx
zK6aAMq=lm-2Bm!a4N!!jJ5Ry9F|Mq0CBg1lu0a(0OYT~9$~nAThFRt@+mWD_@*5+=
zt{SeiX7+-Xm>Rr%xM8oRD^Zwti5n-f+0=3jFED0nNlH`Vh9&T!8M$MpC
zNV1r&miY9O_iEhayTmAhXp~X?UXro?_A^2wLbxMRviHByGBCgSs&g;;Eunq(df-28Uw;*@5sOZD%mqB3
zn$65e@aL%ZP`{|CeH}B$mWhqhRtR)p2}MeIf84f$Zme(D@y`~eDP-Gk4At~j&^7D6
zi1^lThZ??$&4zoTyh^BWi14KhfCPd1NPtT_YMQXjkpG#HX$+BL*A7bvRfOq%hAATg;30KRT+z6f|ey~P;CK0VtR9Xa=DL`f~
z+;A?+U~$js+?OUKj|(FX1tdf(ayr|7l1o^A5b=4ukbHJ<{^mFvIV2*$oOS%GUrFW<
zd>sy12t0peo!D>5pB01&?=?(W`(O?~40I>!MKo~8)llG>$_h18Gi-5H
zBB}X%k7d!NAU_Z5J1_KKS}hlM%qM@qmUb?a>ZnlxTnhuR=mT$eU?~)e%dvYyixVJzT)_uhGNw>u
zVIQYLF~Vrja4JcoM&SjKS(MU^S|y0Pji6+zrf6x}#zoCpB+y;oriH32seV~4bw9);
z5BQN5Z~{?Cs;k?3SDyxYSoRxj#UGYA%sPvsj5#^Vgk8E>{bkIX{<31=;dc9o3R17trD`pcz>)o2mLs+&;nni3Y@^th8}_XzfO{v($Ed3sTH*^Y*#SQFl&Si
z20}Sm&nr~4ATfhapQ3VT;E{nBn;LAPi4ICz2Q{AD{D4nyqKThI3&pE`0TXmy!Em}%
zf>`aB9hf{_&e$K|^C4+7RdED4@F~FWaeGmOgd!by8}>Nrs9#9C%|hl7VeLrZnmE8O
zc8Q>L5Wgmb$vU#T2xS^n<>jH~z#dRd8N~I!%O;xIgCqBrSc!URn$X6Bc7mF!Y?k1x
zl4o%K@-C=%8Qc3}YSTt4f=K@+HJ$uY2P+pz%mYKBdBV^UH=oSp=YcAx$^
z6W1U=lUPp(zLtB8&jq+AP<`Cwz8}SyN|;y---tKSj9Md3daO!M@F;CqJLMX1Hl+;O8D*l)uzD3s{R(dk7egsk%Et+8
z_8fN`I2gxxH5HNzmMZM|9b%%{ozt1egl+8zq4;KD;#xz=?j!GgZ#tr8VpX>Tti4%z
z>MlWYKU)9rHxy9AqLUi{`orL13SnW#V9{K-dJ$=L_G(|_1r-hiIPDVaAv0@EjoD1*
z&UT^J5TR&_@hZpTO5P*I$V(w099&KyE+dVbHqaEE2W|GE=!;x|aR#d^maqy;!+$Hd
zifbCm_h~9$)smico^p3|h1eRY;m}j|^KAO5fZ+);3>)Sd2URF#k{my`+n;%=+3&|l
zF(X(8fd~Ex;hWk3m}qe7qnZh$B0RQTob=Cg(+*X~$Ar?^<|#^
zF3g-4WGNky4c?|kGHSLzdbZ{=uT1)&dw&?$6b6R{BjINZ_KFLGjQA(S^_>SH@->d~
zG%DRnCiKSF@rJ?`3|@d)Hp&>-vTPNntdtHmCLo@H?Rru>CcuR>koK6Rw5kn)iQ40e
zp25;0x|Rscz=`ixW(~)3$ICZJ?Y!YZoc9KxTezvn;rV&+G0>or)zm^$_cAlzDAnGN
zyi%lPH3X5B2=#^(Hp|HMaZf_n&KXzQ7lTFX!4qHxX){ZC7A{4NrscI41hcxgIR2nJ
z!>T95wEc`Eht(aBy|i7GHAcG;OV%U-Y;)X2W+}W3=l|BNSxg2k-zOjWaU`pch9!X#
zAK(~H;Em=E;o$$(t;$~(u~GHHLsIo1-Z-F9q;7;XJTsD}l9@2X3Zr#Qw>wk8ORZ-6_yEN>YBOI4muI!tFxHH`M)wO%+%Xh_g+1h8wr&N?c${udLK4O!n
zmuiMFAqW13JJyk9>avEizRZUnuGn9<}1fAkzr5sH}~5Xrm=MSQ5|d
znCdJ!RwJ6<=Im=BK?;4Co$DD(;-otqE|bRxA(~WK5sYc!7Lfgk)Gv#k=P_dBF!N8s
zFrtbt1Y@~xeVL#xrZ^}zbjUJT%sX5>Ej*D)Zqnmf(l~SvHPgTYO*Po?>wc697Qw?<
zNN5X4cpr_Y*ZXL5g0giob_%DWNww3i0>0FtZevw2Zfnshhv?2s?s~8Tb~>{RyFPrV
zFr3%Hcj)lHELZ%%I@$M?R
z-%r$d>swXfToZNnS>`IpQ<+B2G)kfjpb9llj`+T=H%lYj
z)E_cuoWh(2HvMNzO}K_-_QN-sjL%*D(z_Dfc4Fkkl*VQZ#W9ZaL;!2#YC3InrGYRP867
zgQ`|S#e3c!)5j?4OBNTw+LMBepb`-Xu#^_WA@+zw;nGy8p(Q|y^|q}U#ocmofe9UD
zW#ak8r(|+rxk$XYjXktV^HU{j9n&vS1Qt5(r
z*^Gvt7HHqF7J2#r6dV-$h-_O|(+a)%;Ow{0WyUCd+MtA->
z%i(mmAjF{9uozl5M8uXBTj?(dQ~4<2^ZQ6m`40z(;wwO0P}9RYjOa^=@jI?36<8!F@P>%IQ&L$@o+6kXpdtd4S7P5=m}jsw?>4QwR>C4d45U`
zl(B@#YU-1glVMT)fmoE)?7W83;^}StYskxd9-JHjXoHG(BoXsBO-%~c6!p{)!u%BW
z2nu~)09gRE566!)`1bW=d?h$HKfmC>PLKd58_
zpMWZR6#FV=VciWXsZQavIs+CzFQ)a{Kvf3A##f^kzHO$eh}Way-#=Q@5y_W8u&J(w
z6es6VV_jtZRX9XICu7sJk+$9Ktc^EH6#KO^uqgcc=gmK}Q7bAvB3q8AXQklq#Ed+T
zdVL$q)XlA`%;8tqLh4>pDz>o7^M`_VsTDrZv)ph{sh?+T(6<
zL5`=Gvs7A{YNr{!_`L%;(w*UXY=iQ(%$6P0DKup;dT1JB`uB7UW=J9_ZEGwuYh;wk
z7G(FC36|%LKdT=x`jVNUFY$3wV)?mjuMrArmeKskBi(C@6=mr};Z^pi>ny}jV(xJ&
z@-Hi5O9xs4h!oeapq
zL5#{LFr#`JYgHly`wHD=(3xUROT%Z=*%0b-dAor+_<|Zp8;~o=ez$?D_sOf;7^r=U
zVQ1_qrNw->Drz)u6QtNyUI(?#`Rz{8DBKJFV8X_1BDN1aT5osa?~ZhSyEXAOH$&B
z$^va*OfMwamRXC*G~IISxf_Uz-nhyO2?{NN)l8ZL@zBB98BC^W1;6Hd4WtLyy6|O^
zT<5>&WE$kBS8^UzD=Xxqy;&@lB80IG4?H1~0q~>df+~;|$`^4>`Lsb!k^Ptyzs}dz
z$~gAFq76584AZ~S9^f6Nt!{)@%m?e~;ViKrj{WXCN*756(BZbgYCVI7<1QOm
zI^yi0!Ov*JXbf7UN{u0fYq=CGCZ}FzJ|kp5cV{Yd!t|XNAoRGX~Eh#vk#
z&0=$$DUnxQw)2OU#$Cae!hjlcXwE(WbgyQxqjQ9EoS5JN6c}ga#|E!>YQW
z^+&=(bQyp_o@-h=GYmEeQ-jL6byn!(hZa?z_K;djRPj!FWIpAC6=bIdYhiH7meCki
zcsZ$bQu0u8r@-*2t`)ddLCNzciN#{V_C_ZOVl4|g#B&3Uc-~^(w;{~$8dH_CV5H$=
z^}`z!iE-<)O$4bD>$d`?X>_f|aoR4HS!(Vr2{F6pLak_6!)-8Rzu+17ON0pDhLj=W
zRGxa(E0@FUitldob27_6T96WHn`1zs<0a-xhL0HkTGw6ZP8ir0a-elKDz
z8IwuX=kb+3s@Q%);oYz-```;+tm-9bj$RqG-Y9M%MwW&Go0y(yxd#p%E;p_!U&Zg|
z20QJ(#VzTLQ2V!;IjCk1gd{9G$cN+u+Zy6s`GKa0c(-ifrA-vcwVRwh2
ziMxLdUz0-*S}0fg%3>poUUHDQEY6(U8wqOXmb45LRiH7KKC3tugF$yXYv5wkQ(H%F
z6%($-1FQ3_U{byg?VH7f5ZVEirPoh7C#~+lmBZ)F``5x1>Q?%gm3!o+_
zKL5~oT*tj8ynJ=5Z7z!3H^y@b+IPQVUF>aTil38G*xie$F=c=+DlYGuBO-6u0W
zhNcOcc^kv%gkK#H{BlMSM652Hs`D|%L#A7yutZSpgcpebTv4#n-|2wzVfM`S@p7=P
zpOz|T_6a$^e9x^o9qN}Fio)?bH+95g4?hmcMD_~o&9UTLzSTj$W1QC+K&h2N9J*Re
zg0>@tc`M*qMak$fuPkW>2o^wpRFP2~a{Ai`^mt6K#0vPizy>Oi+!2wBd5aws<2*rb
zV3P(4rUv&lp^OdIIHswKrKRkMEp_4xK4crW_I#R_vHwMkf(Hl%os&620tp$Pq6}7q
z^yFcLTi^=w9BhmE<=B=*aeBmelY&CFW>7q7&z(p{+N_UvN&%$bjME60P0O~|!U0{5
zghj8rL8oyskeQJ5^FP7>l0i=~6yo_5XU1(&Wx1ZB2UT4?j?P+Q1~WYXk*Xq}e9@T;
zMhoIhSGP|k+5GLd9k{Lujx9I1x`K&=Q)gvQ!XTP^FT#vjVxnw$gZL+BfVuR
YFwV018>eah{FzKrR8FK)$iV;q0L2
+
+`libminisketch` is an optimized standalone MIT-licensed library with C API for constructing and decoding *set sketches*, which can be used for compact set reconciliation and other applications.
+It is an implementation of the PinSketch[[1]](#myfootnote1) algorithm. An explanation of the algorithm can be found [here](doc/math.md).
+
+## Sketches for set reconciliation
+
+Sketches, as produced by this library, can be seen as "set checksums" with two peculiar properties:
+* Sketches have a predetermined capacity, and when the number of elements in the set is not higher than the capacity, `libminisketch` will always recover the entire set from the sketch. A sketch of *b*-bit elements with capacity *c* can be stored in *bc* bits.
+* The sketches of two sets can be combined by adding them (XOR) to obtain a sketch of the [symmetric difference](https://en.wikipedia.org/wiki/Symmetric_difference) between the two sets (*i.e.*, all elements that occur in one but not both input sets).
+
+This makes them appropriate for a very bandwidth-efficient set reconciliation protocol. If Alice and Bob each have a set of elements, and they suspect that the sets largely but not entirely overlap,
+they can use the following protocol to let both parties learn all the elements:
+* Alice and Bob both compute a sketch of their set elements.
+* Alice sends her sketch to Bob.
+* Bob combines the two sketches, and obtains a sketch of the symmetric difference.
+* Bob tries to recover the elements from the difference sketch.
+* Bob sends every element in the difference that he has to Alice.
+
+This will always succeed when the size of the difference (elements that Alice has but Bob doesn't plus elements that Bob has but Alice doesn't) does not exceed the
+capacity of the sketch that Alice sent. The interesting part is that this works regardless of the actual set sizes—only the difference matters.
+
+If the elements are large, it may be preferable to compute the sketches over *hashes* of the set elements. In that case an additional step is added to the protocol, where Bob also sends the hash
+of every element he does not have to Alice, who responds with the requested elements.
+
+The doc/ directory has additional [tips for designing reconciliation protocols using libminisketch](doc/protocoltips.md).
+
+## Evaluation
+
+
+
+
+
+**The first graph** above shows a benchmark of `libminisketch` against three other set reconciliation algorithms/implementations. The benchmarks were performed using a single core on a system with an Intel Core i7-7820HQ CPU with clock speed locked at 2.4 GHz. The diagram shows the time needed for merging of two sketches and decoding the result. The creation of a sketch on the same machine takes around 5 ns per capacity and per set element. The other implementations are:
+* [`pinsketch`](https://www.cs.bu.edu/~reyzin/code/fuzzy.html), the original PinSketch implementation.
+* [`cpisync`](https://github.com/trachten/cpisync), a software project which implements a number of set reconciliation algorithms and protocols. The included benchmark analyzes the non-probabilistic version of the original CPISync algorithm[[5]](#myfootnote5) only.
+* A high-performance custom IBLT implementation using 4 hash functions and 32-bit checksums.
+
+For the largest sizes currently of interest to the authors, such as a set of capacity 4096 with 1024 differences, `libminisketch` is forty-nine times faster than `pinsketch` and over eight thousand times faster than `cpisync`. `libminisketch` is fast enough on realistic set sizes for use on high-traffic network servers where computational resources are limited.
+
+Even where performance is latency-limited, small minisketches can be fast enough to improve performance. On the above i7-7820HQ, a set of 2500 30-bit entries with a difference of 20 elements can be communicated in less time with a minisketch than sending the raw set so long as the communications bandwidth is 1 gigabit per second or less; an eight-element difference can be communicated in better than one-fifth the time on a gigabit link.
+
+**The second graph** above shows the performance of the same algorithms on the same system, but this time keeping the capacity constant at 128, while varying the number of differences to reconcile between 1 and 128. It shows how `cpisync`'s reconciliation speed is mostly dependent on capacity, while `pinsketch`/`libminisketch` are more dependent on number of differences.
+
+**The third graph** above shows the size overhead of a typical IBLT scheme over the other algorithms (which are near-optimal bandwidth), for various levels of failure probability. IBLT takes tens of times the bandwidth of `libminisketch` sketches when the set difference size is small and the required failure rate is low.
+
+**The fourth graph** above shows the effect of the field size on speed in `libminisketch`. The three lines correspond to:
+* CLMUL 64-bit: Intel Core i7-7820HQ system at 2.4 GHz
+* Generic 64-bit: POWER9 CP9M06 system at 2.8 GHz (Talos II)
+* Generic 32-bit: Cortex-A53 at 1.2 GHz (Raspberry Pi 3B)
+
+It shows how CLMUL implementations are faster for certain fields (specifically, field sizes for which an irreducible polynomial of the form *xb + x + 1* over *GF(2)* exists, and to a lesser extent, fields which are a multiple of 8 bits). It also shows how (for now) a significant performance drop exists for fields larger than 32 bits on 32-bit platforms. Note that the three lines are not at the same scale (the Raspberry Pi 3B is around 10x slower for 32-bit fields than the Core i7; the POWER9 is around 1.3x slower).
+
+Below we compare the PinSketch algorithm (which `libminisketch` is an implementation of) with other set reconciliation algorithms:
+
+| Algorithm | Sketch size | Decode success | Decoding complexity | Difference type | Secure sketch |
+| ----------------------------------------------------- | ------------------------- | ---------------| ------------------- | --------------- | ------------- |
+| CPISync[[2]](#myfootnote2) | *(b+1)c* | Always | *O(n3)* | Both | Yes |
+| PinSketch[[1]](#myfootnote1) | *bc* | Always | *O(n2)* | Symmetric only | Yes |
+| IBLT[[6]](#myfootnote1)[[7]](#myfootnote1) | *αbc* (see graph 3) | Probabilistic | *O(n)* | Depends | No |
+
+* **Sketch size:** This column shows the size in bits of a sketch designed for reconciling *c* different *b*-bit elements. PinSketch and CPISync have a near-optimal[[11]](#myfootnote11) communication overhead, which in practice means the sketch size is very close (or equal to) *bc* bits. That is the same size as would be needed to transfer the elements of the difference naively (which is remarkable, as the difference isn't even known by the sender). For IBLT there is an overhead factor *α*, which depends on various design parameters, but is often between *2* and *10*.
+* **Decode success:** Whenever a sketch is designed with a capacity not lower than the actual difference size, CPISync and PinSketch guarantee that decoding of the difference will always succeed. IBLT always has a chance of failure, though that chance can be made arbitrarily small by increasing the communication overhead.
+* **Decoding complexity:** The space savings achieved by near-optimal algorithms come at a cost in performance, as their asymptotic decode complexity is quadratic or cubic, while IBLT is linear. This means that using near-optimal algorithms can be too expensive for applications where the difference is sufficiently large.
+* **Difference type:** PinSketch can only compute the symmetric difference from a merged sketch, while CPISync and IBLT can distinguish which side certain elements were missing on. When the decoder has access to one of the sets, this generally doesn't matter, as he can look up each of the elements in the symmetric difference with one of the sets.
+* **Secure sketch:** Whether the sketch satisfies the definition of a secure sketch[[1]](#myfootnote1), which implies a minimal amount about a set can be extracted from a sketch by anyone who does not know most of the elements already. This makes the algorithm appropriate for applications like fingerprint authentication.
+
+## Building
+
+The build system is very rudimentary for now, and [improvements](https://github.com/sipa/minisketch/pulls) are welcome.
+
+The following may work and produce a `libminisketch.a` file you can link against:
+
+```bash
+git clone https://github.com/sipa/minisketch
+cd minisketch
+./autogen.sh && ./configure && make
+```
+
+## Usage
+
+In this section Alice and Bob are trying to find the difference between their sets.
+Alice has the set *[3000 ... 3009]*, while Bob has *[3002 ... 3011]*.
+
+First, Alice creates a sketch:
+
+```c
+#include