## Adi Shamir’s Threshold Scheme

### June 17, 2011

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 *i*th *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.

I also wrote a version in Common Lisp, since I’m trying to learn a new language. If there are any CL gurus in the audience, I’d appreciate any feedback!

For good measure, here it is in C and Haskell.

Here’s my solution in C:

Err, that should rather be

I changed github accounts, so the gists I linked to above are dead (apologies). Here are new ones:

Common Lisp

Haskell

C, version 1

C, version 2 (improved with a little help from Razvan’s solution).

= [22 · 2 · 11−1 · 21 · 7−1] + [8 · 14 · 12−1 · 21(mark) · 19−1] + [5 · 14 · 16−1 · 21 · 4−1] (mod 23)

= [22 · 2 · 21 · 21 · 10] + [8 · 14 · 2 · 23(mark) · 17] + [5 · 14 · 13 · 21 · 6] (mod 23)

these lines has been copied from the solution example given above.

see the mark written portion in both the lines. how did 21 got transformed to 23?

please explain