|
;; 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] |
