# Shamir's Secret Sharing

# cat shamir_secret_sharing.py
import random
import sys

secret = int(sys.argv[1])
components = int(sys.argv[2])
at_least = int(sys.argv[3])

def is_prime(number):
        return all(number % i for i in xrange(2, int(number ** 0.5) + 1))

def get_next_prime(number):
        next_number = number + 1
        while not is_prime(next_number):
                next_number += 1
        return next_number

def egcd(a, b):
        if a == 0:
                return (b, 0, 1)
        else:
                g, y, x = egcd(b % a, a)
                return (g, x - (b // a) * y, y)

def modinv(a, m):
        g, x, y = egcd(a, m)
        if g != 1:
                return None  # modular inverse does not exist
        else:
                return x % m

def split(prime, number, available, needed):
        coef = [number]
        shares = []
        for i in xrange(1, needed):
                coef.append(random.randrange(0, prime))
        for x in xrange(1, available + 1):
                # coef[0]x^0 + coef[1]x^1 + ... + coef[n]x^n
                accum = coef[0]
                for e in xrange(1, needed):
                        accum += ((coef[e] * pow(x, e, prime)) % prime)
                accum = accum % prime
                shares.append((x, accum))
        return coef, shares

def get_random_shares(available, number, shares):
        newshares = []
        selected = []
        for i in xrange(0, number):
                r = random.randrange(0, available)
                while r in selected:
                        r = random.randrange(0, available)
                selected.append(r)
                newshares.append(shares[r])
        return newshares

def merge(prime, shares):
        accum = 0
        # Lagrange's interpolation
        for formula in xrange(0, len(shares)):
                numerator = denominator = 1
                for count in xrange(0, len(shares)):
                        if formula == count:
                                continue
                        startpos = shares[formula][0]
                        nextpos = shares[count][0]
                        numerator = (numerator * -1 * nextpos) % prime
                        denominator = (denominator * (startpos - nextpos)) % prime
                value = shares[formula][1]
                accum = (prime + accum + (value * numerator * modinv(denominator, prime))) % prime
        return accum

print 'Secret =', secret
prime = get_next_prime(secret)
print 'Prime =', prime
coef, shares = split(prime, secret, components, at_least)
print 'Coefficients =', coef
print 'Subsecrets =', shares
newshares = get_random_shares(components, at_least, shares)
print 'Random subsecrets =', newshares
secret = merge(prime, newshares)
print 'Recovered secret =', secret

# python shamir_secret_sharing.py 1337 7 3
Secret = 1337
Prime = 1361
Coefficients = [1337, 255, 752]
Subsecrets = [(1, 983), (2, 772), (3, 704), (4, 779), (5, 997), (6, 1358), (7, 501)]
Random subsecrets = [(6, 1358), (4, 779), (7, 501)]
Recovered secret = 1337

Reference

https://en.wikipedia.org/wiki/Shamir%27s_Secret_Sharing