# -*- coding: utf-8 -*-
# ENGG 2440A / ESTR 2004 Lecture 11 -- python code

from random import *
from decimal import *

# a prime number
q = 5915587277

# output d + 1 shares for secret s modulo q   
def share(d, q, s):
    c = coeffs(d, q, s)
    shares = []
    for i in range(1, d + 2):
        shares.append(eval_poly(c, i, q))
    return shares

# recover the secret from the shares using the formula
# secret = sum shares[i] (d + 1 choose i) (-1)^{d - i + 1} for i = 1 to d + 1
def recover(shares, d, q):
    secret = 0
    D = d + 1
    binom = D 
    for i in range(1, D + 1):
        secret += binom * shares[i - 1]
        binom = -binom * (D - i) / (i + 1)
    return secret % q

# statistics on pairs of shares for three parties
def stats(q = 5, s = 1, trials = 100000):
    # number of pairs observed for every two parties
    AB = [[0] * q for x in range(q)]
    AC = [[0] * q for x in range(q)]
    BC = [[0] * q for x in range(q)]

    # collect statistics    
    for trial in range(trials):
        shares = share(2, q, s)
        AB[shares[0]][shares[1]] += 1
        AC[shares[0]][shares[2]] += 1
        BC[shares[1]][shares[2]] += 1
            
    return AB, BC, AC

# coefficients of a secret sharing polynomial
def coeffs(d, q, s):
    c = []
    for i in range(d):
        c.append(randint(0, q - 1))
    c.append(s)
    return c

# calculate the value of a polynomial with the given vector of coefficients at x modulo q
def eval_poly(coeffs, x, q):
    y = 0
    d = len(coeffs) - 1
    for i in range(d):
        y = (y + coeffs[i]) * x
    y = y + coeffs[d]
    return y % q