Sage code for Shamir's Secret Sharing Algorithm and Reconstruction using Lagrange Interpolation. Go to cocalc to run the code. You can also install sage on your local machine.

reset()
        
# Random number generation in Zp
def rand_gen(p) :
    Zp = Integers(p)
    return Zp.list()[randint(0,p-1)]
        
# Return a  random polynomial in Zp[X] of degree <= (t-1)
# And through the point (0,s)
def rand_poly_share(s,t,p) :
    Zp = Integers(p)
    R.<X> = Zp[]
    coeffs = []
    for i in [0,..,t-1] :
        if (i==0) :
            coeffs.append(Zp.list()[s])
        else :
            coeffs.append(rand_gen(p))
    #print coeffs
    PX = coeffs[0]
    for i in [1,..,t-1] :
        PX = PX + coeffs[i]*(X^i)
    return PX
            
# Secret share s
# into n-parties
# t-parties can reconstruct
# over Zp
# shares[0] has the secret for sanity check purposes
def shamir_share(s,n,t,p) :
    if (s<0 or s>p-1) :
        return -1
    if (n<1 or n>p-1) :
        return -1
    # Generate a random polynomial consistent with the secret
    PX = rand_poly_share(s,t,p)
    print "...[Internal State] Secret Sharing Polynomial used:", PX
    # Evaluate the polynomial to obtain the shares
    shares = []
    for i in [0..n] :
        shares.append(PX(X=i))
    return shares
                    
                    
# Lagrange Interpolate the list of coordinates using Zp[X] polynomial
def lagrange_interpolation(coords,p) :
    t = len(coords)
    Zp = Integers(p)
    R. = Zp[]
    x_tup, y_tup = zip(*coords)
    x_list = list(x_tup)
    y_list = list(y_tup)
    polys = []
    for i in [0,..,t-1] :
        polys.append(1)
        for j in [0,..,t-1] :
            if (j == i) :
                polys[i] = polys[i]* y_tup[i]
            else :
                polys[i] = polys[i]* (X-x_tup[j]) / (x_tup[i] - x_tup[j])
    print "...[Internal State] List of P_i(X) polynomials:", polys
    PprimeX = 0
    for i in [0,..,t-1] :
        PprimeX = PprimeX + polys[i]
    print "...[Internal State] Reconstructed Polynomial:", PprimeX
    return PprimeX
                        
                        
# Reconstruction algorithm
def reconstruct_secret(coords,p) :
    PprimeX = lagrange_interpolation(coords,p)
    return PprimeX(X=0)
                        
                        
                        
s = 10 ### Input by User ###
n = 16 ### Input by User ###
t = 5  ### Input by User ###
p = 17 ### Input by User ###
                        
print "Secret to be shared:", s
print "Number of parties:", n
print "Number of parties that can reconstruct the secret:", t
print "The prime used as the field size:", p
                        
secret_shares = shamir_share(s,n,t,p)
print "\nThe Secret Shares of Parties:", secret_shares
                        
# Indices of t parties who want to reconstruct the secret
ind = [1,3,4,6,15] ### Input by User ###
print "\nList of parties who want to reconstruct the secret", ind
                        
vals = []
for i in ind :
    vals.append(secret_shares[i])
# Coordinates of reconstruction parties and their respective secret shares
coords = zip(ind,vals)
print "\nCoordinates to Interpolate:" , coords
sprime = reconstruct_secret(coords,p)
print "Reconstructed Secret:", sprime
                        
print "\nIs the reconstructed secret identical to the original secret?", s==sprime