From 69c966134cdba3c8a38038b9fd28a0c78e88cc73 Mon Sep 17 00:00:00 2001
From: Determinant <tederminant@gmail.com>
Date: Thu, 22 Aug 2024 17:39:58 -0700
Subject: support Shamir's secret sharing with minimal code (`shamir.py`)

---
 shamir.py | 113 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 113 insertions(+)
 create mode 100644 shamir.py

(limited to 'shamir.py')

diff --git a/shamir.py b/shamir.py
new file mode 100644
index 0000000..c7d606a
--- /dev/null
+++ b/shamir.py
@@ -0,0 +1,113 @@
+# Code modified from https://github.com/kurtbrose/shamir/tree/master
+
+'''
+An implementation of shamir secret sharing algorithm.
+
+To the extent possible under law,
+all copyright and related or neighboring rights
+are hereby waived.  (CC0, see LICENSE file.)
+
+All possible patents arising from this code
+are renounced under the terms of
+the Open Web Foundation CLA 1.0
+(http://www.openwebfoundation.org/legal/the-owf-1-0-agreements/owfa-1-0)
+'''
+
+from __future__ import division
+import random
+import functools
+
+# 12th Mersenne Prime is 2**127 - 1
+# use the prime to be compatible with the Shamir tool developed by Ava Labs:
+# https://github.com/ava-labs/mnemonic-shamir-secret-sharing-cli/tree/main
+# This is a 257-bit prime, so the returned points could be either 256-bit or
+# (infrequently) 257-bit.
+_PRIME = 187110422339161656731757292403725394067928975545356095774785896842956550853219
+# 13th Mersenne Prime is 2**521 - 1
+
+def _eval_at(poly, x, prime):
+    'evaluate polynomial (coefficient tuple) at x'
+    accum = 0
+    for coeff in reversed(poly):
+        accum *= x
+        accum += coeff
+        accum %= prime
+    return accum
+
+
+def split(secret, minimum, shares, prime=_PRIME):
+    '''
+    Generates a random shamir pool, returns
+    the secret and the share points.
+    '''
+    randint = functools.partial(random.SystemRandom().randint, 0)
+    if minimum > shares:
+        raise ValueError("pool secret would be irrecoverable")
+    poly = [randint(prime) for i in range(minimum - 1)]
+    poly.insert(0, secret)
+    return [_eval_at(poly, i, prime) for i in range(1, shares + 1)]
+
+
+# division in integers modulus p means finding the inverse of the denominator
+# modulo p and then multiplying the numerator by this inverse
+# (Note: inverse of A is B such that A*B % p == 1)
+# this can be computed via extended euclidean algorithm
+# http://en.wikipedia.org/wiki/Modular_multiplicative_inverse#Computation
+def _extended_gcd(a, b):
+    x = 0
+    last_x = 1
+    y = 1
+    last_y = 0
+    while b != 0:
+        quot = a // b
+        a, b = b,  a%b
+        x, last_x = last_x - quot * x, x
+        y, last_y = last_y - quot * y, y
+    return last_x, last_y
+
+
+def _divmod(num, den, p):
+    '''
+    compute num / den modulo prime p
+    To explain what this means, the return
+    value will be such that the following is true:
+    den * _divmod(num, den, p) % p == num
+    '''
+    inv, _ = _extended_gcd(den, p)
+    return num * inv
+
+
+def _lagrange_interpolate(x, x_s, y_s, p):
+    '''
+    Find the y-value for the given x, given n (x, y) points;
+    k points will define a polynomial of up to kth order
+    '''
+    k = len(x_s)
+    assert k == len(set(x_s)), "points must be distinct"
+    def PI(vals):  # upper-case PI -- product of inputs
+        accum = 1
+        for v in vals:
+            accum *= v
+        return accum
+    nums = []  # avoid inexact division
+    dens = []
+    for i in range(k):
+        others = list(x_s)
+        cur = others.pop(i)
+        nums.append(PI(x - o for o in others))
+        dens.append(PI(cur - o for o in others))
+    den = PI(dens)
+    num = sum([_divmod(nums[i] * den * y_s[i] % p, dens[i], p)
+               for i in range(k)])
+    return (_divmod(num, den, p) + p) % p
+
+
+def combine(shares, prime=_PRIME):
+    '''
+    Recover the secret from share points
+    (x,y points on the polynomial)
+    '''
+    if len(shares) < 2:
+        raise ValueError("need at least two shares")
+    x_s, y_s = zip(*shares)
+    return _lagrange_interpolate(0, x_s, y_s, prime)
-- 
cgit v1.2.3-70-g09d2