A sample code for "Shamir's secret sharing" (http://goo.gl/7PjHF) with the Lagrange interpolation. Written in Clojure.

shared-secret.clj

;; to represent a polynomial, a vector [a b c d] means a + bx + cx^2 + dx^3
;; gen-poly makes a polynomial with n random coefficients (max value: max-coeff)
(defn gen-poly [n max-coeff]
  (vec (repeatedly n #(rand-int max-coeff))))

;; if p is negative, returns n
;; else returns n mod p
(defn modp [n p]
  (if (<= p 0)
    n
    (mod n p)))

;; calculate a polynomial with x
;; (cal-poly [4 2 1] 2) => 12
;; f(x) = x^2 + 2x + 4
;; if x equals 2, 4 + 4 + 4 = 12
(defn cal-poly
  ([poly x] (cal-poly poly x -1))
  ([poly x m]
   (loop [p poly
          z 1
          cal 0]
     (if-not (seq p)
       cal
       (recur (rest p) (modp (* z x) m) (modp (+ cal (* (first p) z)) m))))))

;; add n degrees with x
;; (fill-with-x [3 4 2] 2 5) => [3 4 2 5 5]
(defn fill-with-x [poly n x]
  (if (<= n 0)
    poly
    (recur (conj poly x) (dec n) x)))

;; add two polynomials
(defn add-poly
  ([p1 p2] (add-poly p1 p2 -1))
  ([p1 p2 m]
   (if (> (count p1) (count p2))
     (add-poly p2 p1 m)
     (let [d (- (count p2) (count p1))
           q (fill-with-x p1 d 0)]
       (->> (mapv + q p2)
            (mapv #(modp % m)))))))

;; multiply two items
;; (mul-two-item [[1 2] [3 4]]) => [4 8]
;; this represents 2x * 4x^3 = 8x^4
(defn mul-two-item [xs]
  (let [a (first xs)
        b (second xs)]
    [(+ (first a) (first b)) (* (second a) (second b))]))

;; multiply two polymomials
(defn mul-poly
  ([p1 p2] (mul-poly p1 p2 -1))
  ([p1 p2 m]
   (->> (for [x (map-indexed #(vector %1 %2) p1)
              y (map-indexed #(vector %1 %2) p2)]
          [x y])
        (map mul-two-item)
        (map #(hash-map (first %) (second %)))
        (apply merge-with +)
        (mapv val)
        (mapv #(modp % m)))))

;; drop n-th item and returns the vector
(defn dropv-nth [n xs]
  (vec (concat (subvec xs 0 n) (subvec xs (inc n) (count xs)))))

;; Euler's modular multiplicative inverse
;; a^phi(m) = 1 mod m
;; a^(phi(m) - 1) = a^-1 mod m
;; phi(m) = m - 1 if m is prime
;; a^(m - 1 - 1) = a^(m-2) = a^-1 mod m
(defn multiplicative-inverse-euler [n p]
  (bigint (.modPow (biginteger n) (biginteger (- p 2)) (biginteger p))))

;; calculate an item for Lagrange interpolation, L^n_t(t)
;; xs - a vector. x values
;; idx - represents t in L^n_t(t)
;; y - p(x)
;; m - modulus
(defn cal-lt
  ([xs idx y] (cal-lt xs idx y -1))
  ([xs idx y m]
   (let [upper (->> xs
                    (mapv #(vec (list (- %) 1)))
                    (dropv-nth idx)
                    (reduce #(mul-poly %1 %2 m)))
         w (nth xs idx)
         lower (modp (->> xs
                         (mapv #(- w %))
                         (dropv-nth idx)
                         (reduce *))
                     m)]
     (if (pos? m)
       (mapv #(modp (* (modp (* % y) m) (multiplicative-inverse-euler lower m)) m) upper)
       (mapv #(/ (* % y) lower) upper)))))

;; Lagrange interpolation
(defn lagrange-interpolation
  ([ps] (lagrange-interpolation ps -1))
  ([ps m]
   (let [xs (mapv first ps)]
     (loop [res [0]
            lt ps
            idx 0]
       (if-not (seq lt)
         res
         (recur (add-poly res (cal-lt xs idx (second (first lt)) m) m) (rest lt) (inc idx)))))))

(defn gen-x-fx-pair
  ([poly x] (gen-x-fx-pair poly x -1))
  ([poly x m]
   [x (cal-poly poly x m)]))

;; shared keys MUST have different values
(defn gen-n-distinct-vec [n limit]
  (loop [xs []]
    (if (= (count xs) n)
      xs
      (recur (vec (set (conj xs (inc (rand-int (dec limit))))))))))

(defn create-shared-secret [k n max-coeff m]
  (let [poly (gen-poly k m)]
    (list poly (map #(gen-x-fx-pair poly % m) (gen-n-distinct-vec n 100)))))

;; (create-shared-secret 6 10 1000 4523)
;; ([1804 3539 3114 2303 1721 3165] ([32 3771] [70 1282] [80 271] [50 564] [19 2022] [21 1447] [59 2831] [60 134] [29 2374] [95 2622]))

;; (lagrange-interpolation [[32 3771] [80 271] [19 2022] [59 2831] [29 2374] [95 2622]] 4523)
;; [1804N 3539N 3114N 2303N 1721N 3165N]