Programming Praxis

We first import two functions from the random module:

from random import randrange, sample

These will furnish us with the random integers we require:randrange(1,p) returns a uniformly distributed random integer in the half-open interval [1, p), while sample(xrange(1, p), n) will select n of the integers in [1, p) with no repeats.

We build our polynomial via a modified version of Horner’s Method:

def horner_mod(coeffs, mod):
    return lambda x: (x,
        reduce(lambda a, b: a * x + b % mod, reversed(coeffs), 0) % mod)

This function, given a list of coefficients and a modulus, returns the function that maps x to the pair (x, y) where y is the result of running $x$ through the polynomial modulo mod whose coefficients are listed in coeffs.

From there, Shamir’s Threshold scheme involves building our list of coefficients (the first — really zeroth — of which is the secret S), then evaluating the polynomial we’ve built n times and returning the (x, y) pairs:

def shamir_threshold(S, k, n, p):
    coeffs = [S]
    coeffs.extend(randrange(1, p) for _ in xrange(k - 1))
    return map(horner_mod(coeffs, p), sample(xrange(1, p), n))

Finally, we us Lagrange Interpolation to find the constant term S given a list of k or more (x, y) pairs. For good measure, we first check that we have enough pairs, then use just the first k to find S:

def interp_const(xy_pairs, k, p):
    assert len(xy_pairs) >= k, "Not enough points for interpolation"
    x = lambda i: xy_pairs[i][0]
    y = lambda i: xy_pairs[i][1]
    return sum(y(i) * prod(x(j) * mod_inv(x(j) - x(i), p) % p for j in xrange(k)
        if j != i) for i in xrange(k)) % p

We first define inline functions that allow us to pick out the ith x and y values from our list of pairs, then use them to write the equation for S given on the previous page.

As for testing, let’s encode the words “PRAXIS” as an integer by considering all digits and letters as integers in base 36. We’ll take p to be the next prime after S — this choice of p is arbitrary, it just needs to be bigger than S. To keep things manageable, let’s take n to be 20 and k to be 5, so we’re constructing a degree 4 polynomial and we’ll hand out 20 (x, y) pairs.

if __name__ == "__main__":
    from pprint import pprint                   # Pretty printing
    S = int("PRAXIS", 36); print S              # Prints 1557514036
    n, k, p = 20, 5, 1557514061                 # p is the next prime after S
    xy_pairs = shamir_threshold(S, k, n, p)
    pprint(xy_pairs)                            # Prints all 20 (x, y) pairs
    print interp_const(pairs, k, p)             # Should print 1557514036


which gives the output (on one running; pairs depends on the random polynomial created)

1557514036
[(697286162, 445615394L),
 (471866046, 757728985L),
 (112045393, 1132162792L),
 (397324764, 486286231L),
 (135120894, 1142009194L),
 (508637994, 1556915744L),
 (488738532, 834401917L),
 (1369874096, 1345716686L),
 (91597754, 487556032L),
 (970187759, 341284274L),
 (1102805729, 224871713L),
 (245100902, 1306749801L),
 (413372256, 568733054L),
 (1218343037, 63534734L),
 (442535975, 1060000953L),
 (1173207231, 400308586L),
 (515043844, 141960722L),
 (1162691976, 374990038L),
 (73252341, 785232686L),
 (934671161, 486917357L)]
1557514036

(Ignore the L on the end of the numbers on the right; this just signifies that Python used long integer arithmetic to find them. Python 3.x, the newer iteration of Python, removes the distinction between integers and long integers.)

Our solution used mod_inv from a previous exercise; we also defined the function prod, analogous to sum from the Standard Prelude. You can run the full program at http://codepad.org/NvYVdZFT.

Pages: 1 2