简介

Shamir 密钥分享算法最早是由 Shamir 和 Blackly 在 1970 年基于 Lagrange 插值和矢量方法提出的。

算法有 2 个重要参数：kk 和 nn。nn 表示将明文加密为 nn 个 ShadowShadow，kk 表示 至少需要 kk 个 ShadowShadow 才可以恢复出明文。

加密

若明文为 ss，s∈ Zps\in\ Z_p，pp 为一个大素数。在 GF(p)GF(p) 任取 k−1k-1 个随机数： a1a_1、a2a_2、a3a_3、…、ak−1a_{k-1}，构造如下多项式：

f(x)=s+a1x+a2x2+a3x3+...+ak−1xk−1 mod pf(x) = s + a_1x + a_2x^2 + a_3x^3 + ... + a_{k-1}x^{k-1} \bmod p

任取 nn 个不同的数：x1x_1、x2x_2、x3x_3、…、xnx_n，分别代入多项式得到 nn 个密钥对：

{x1、f(x1)}{x2、f(x2)}{x3、f(x3)}...{xn、f(xn)}\{x_1、f(x_1)\} \\ \{x_2、f(x_2)\} \\ \{x_3、f(x_3)\} \\ ... \\ \{x_n、f(x_n)\}

将这 nn 个密钥对分发给 nn 个持有者。

解密

得到了 nn个密钥对 {x1、f(x1)}\{x_1、f(x_1)\}、{x2、f(x2)}\{x_2、f(x_2)\}、{x3、f(x3)}\{x_3、f(x_3)\}、…、{xk、f(xn)}\{x_k、f(x_n)\} 后，可以列出以下在 GF(p)GF(p) 上的方程组：

{s+a1x1+a2x12+...+ak−1x1k−1=f(x1)s+a1x2+a2x22+...+ak−1x2k−1=f(x2)s+a1x3+a2x32+...+ak−1x3k−1=f(x3)...s+a1xn+a2xn2+...+ak−1xnk−1=f(xn)\left\{ \begin{array}{c} s + a_1x_1 + a_2{x_1}^2 + ... + a_{k-1}{x_1}^{k-1} = f(x_1) \\ s + a_1x_2 + a_2{x_2}^2 + ... + a_{k-1}{x_2}^{k-1} = f(x_2) \\ s + a_1x_3 + a_2{x_3}^2 + ... + a_{k-1}{x_3}^{k-1} = f(x_3) \\ ... \\ s + a_1x_n + a_2{x_n}^2 + ... + a_{k-1}{x_n}^{k-1} = f(x_n) \end{array} \right.

然后就可以用 Lagrange 插值算法求出 ss 了。

python 的加解密代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
from Crypto.Util.number import *
import gmpy2


def product(vals, p):
    return reduce(lambda x, y: x * y % p, vals)


def lagrange_interpolate(x, x_s, y_s, p):
    n = len(x_s)
    assert n == len(set(x_s))  # x_s must be distinct
    num = []
    den = []
    for i in range(n):
        others = x_s[:i] + x_s[i+1:]
        num.append(product((x - o for o in others), p))
        den.append(product((x_s[i] - o for o in others), p))
    dens = product(den, p)
    res = 0
    for i in range(n):
        tmp1 = ((num[i] * dens % p) * y_s[i]) % p
        tmp2 = tmp1 * gmpy2.invert(den[i], p) % p
        res = (res + tmp2) % p 
    res = res * gmpy2.invert(dens, p) % p
    return res


def shamir_sharing_encode(s, k, n, p):
    a = [getRandomInteger(Nbits) for _ in range(k-1)]
    res = []
    for _ in range(n):
        x = getRandomInteger(Nbits)
        fx = s
        for i in range(k-1):
            fx = (fx + a[i] * pow(x, i+1, p)) % p
        res.append((x, fx))
    return res


def shamir_sharing_decode(shares, p):
    x_s, y_s = zip(*shares)
    return lagrange_interpolate(0, x_s, y_s, p)


Nbits = 1024
secret = bytes_to_long('it_is_the_top_secret')
p = getPrime(Nbits)
k = 513
n = 513

shares = shamir_sharing_encode(secret, k, n, p)

s = shamir_sharing_decode(shares, p)
print long_to_bytes(s)
# 2s

当然解密算法也可以用 sage，更方便：

1
2
3
4
5
6
7
8
9
10
P = PolynomialRing(GF(p), 'x')
ret = P(0)
for x, y in shares:
    r = P(1) * y
    for xx, yy in shares:
        if x != xx:
            r = r * P((0 - xx) / (x - xx))
    ret = ret + r
print ret
# 19s

也可以用 sage 算出所有的系数，但相应的解密时间也会很长：

1
2
3
4
5
6
7
8
9
10
P = PolynomialRing(GF(p), 'x')
ret = P(0)
for x, y in shares:
    r = P(1) * y
    for xx, yy in shares:
        if x != xx:
            r = r * P('(x - %d) / (%d - %d)' % (xx, x, xx))
    ret = ret + r
print ret[0]
# 154s

除了 Lagrange 插值，也可以用矩阵。

可以将上面的方程组化为以下在 GF(p)GF(p) 上的矩阵乘法：

(1x1x12⋯x1k−11x2x22⋯x2k−1⋮⋮⋮⋱⋮1xnxn2⋯xnk−1)×(sa1⋮ak−1)=(f(x1)f(x2)⋮f(xn))\begin{pmatrix} 1 & x_1 & {x_1}^2 & \cdots & {x_1}^{k-1} \\ 1 & x_2 & {x_2}^2 & \cdots & {x_2}^{k-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & {x_n}^2 & \cdots & {x_n}^{k-1} \\ \end{pmatrix} \times \begin{pmatrix} s \\ a_1 \\ \vdots \\ a_{k-1} \end{pmatrix} = \begin{pmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_n) \end{pmatrix}

sage 矩阵解法，时间很长：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
k = 513
n = 513

x, y = zip(*shares)

array_1 = []
for i in range(n):
    for j in range(k):
        array_1.append(pow(x[i], j, p))
A = matrix(GF(p), n, k, array_1)

array_2 = []
for i in range(n):
    array_2.append(y[i])
B = matrix(GF(p), n, 1, array_2)

a = A.solve_right(B)
print a[0][0]
# 266s

也可以将之视为多元同余方程组：

{s+a1x1+a2x12+...+ak−1x1k−1 mod p≡f(x1)s+a1x2+a2x22+...+ak−1x2k−1 mod p≡f(x2)s+a1x3+a2x32+...+ak−1x3k−1 mod p≡f(x3)...s+a1xn+a2xn2+...+ak−1xnk−1 mod p≡f(xn)\left\{ \begin{array}{c} s + a_1x_1 + a_2{x_1}^2 + ... + a_{k-1}{x_1}^{k-1} \bmod p \equiv f(x_1) \\ s + a_1x_2 + a_2{x_2}^2 + ... + a_{k-1}{x_2}^{k-1} \bmod p \equiv f(x_2) \\ s + a_1x_3 + a_2{x_3}^2 + ... + a_{k-1}{x_3}^{k-1} \bmod p \equiv f(x_3) \\ ... \\ s + a_1x_n + a_2{x_n}^2 + ... + a_{k-1}{x_n}^{k-1} \bmod p \equiv f(x_n) \end{array} \right.

记系数矩阵的行列式为：

Δ=∣1x1x12⋯x1k−11x2x22⋯x2k−1⋮⋮⋮⋱⋮1xnxn2⋯xnk−1∣≠0\Delta = \begin{vmatrix} 1 & x_1 & {x_1}^2 & \cdots & {x_1}^{k-1} \\ 1 & x_2 & {x_2}^2 & \cdots & {x_2}^{k-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & {x_n}^2 & \cdots & {x_n}^{k-1} \\ \end{vmatrix} \neq 0

要求的是 ss，因此将行列式第一列换为同余方程的值：

Δ0=∣f(x1)x1x12⋯x1k−1f(x2)x2x22⋯x2k−1⋮⋮⋮⋱⋮f(xn)xnxn2⋯xnk−1∣\Delta_0 = \begin{vmatrix} f(x_1) & x_1 & {x_1}^2 & \cdots & {x_1}^{k-1} \\ f(x_2) & x_2 & {x_2}^2 & \cdots & {x_2}^{k-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ f(x_n) & x_n & {x_n}^2 & \cdots & {x_n}^{k-1} \\ \end{vmatrix}

由克拉默法则得：

s≡Δ−1Δ0 mod ps \equiv \Delta^{-1} \Delta_0 \bmod p

sage 脚本如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
fx_num = 513

x, y = zip(*shares)

array_1 = []
for i in range(fx_num):
    for j in range(fx_num):
        array_1.append(pow(x[i], j, p))
delta = matrix(GF(p), fx_num, fx_num, array_1).determinant()

array_2 = [_ for _ in array_1]
for i in range(fx_num):
    array_2[i*fx_num] = y[i]
delta_0 = matrix(GF(p), fx_num, fx_num, array_2).determinant() % p

delta_inverse = inverse_mod(int(delta), p)
res = delta_inverse * delta_0 % p
print res
# 1200s

但是如果维数很大，直接解行列式会很慢。观察发现 Δ\Delta 其实是旋转后的范德蒙行列式。Δ0\Delta_0 是旋转后的范德蒙行列式的变形，第一行不是全 1，而是 f(x)f(x)。可以按第一行展开，余子式中将每一列除以该列元的一次方，可化为范德蒙行列式。计算速度会快非常多。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
fx_num = 513

x, y = zip(*shares)

delta = 1
for i in range(1, fx_num):
    for j in range(0, i):
        t = x[i] - x[j]
        delta = (delta * t) % p

delta_0 = 0
for i in range(fx_num):
    tmp_x = [_ for _ in x]
    tmp_x.pop(i)
    yuzishi = -y[i] if (i+1+1) % 2 else y[i]
    for j in range(fx_num-1):
        yuzishi = (yuzishi * tmp_x[j]) % p
    for j in range(1, fx_num-1):
        for k in range(0, j):
            t = tmp_x[j] - tmp_x[k]
            yuzishi = (yuzishi * t) % p
    delta_0 = (delta_0 + yuzishi) % p

delta_inverse = gmpy2.invert(delta, p)
res = delta_inverse * delta_0 % p
print res
# 224s

Problem

如果在加密 2 段明文 s1s1、s2s2 过程中，使用相同的随机数 aa 以及大素数 pp，还能得到同一个 xx 对应的两个 f(x)f(x)，如果知道一个明文 s1s1 的情况下，可以算出另一个明文。

在 GF(p)GF(p) 上：

s1+a1x+a2x2+a3x3+...+ak−1xk−1=f(x)1s_1 + a_1x + a_2x^2 + a_3x^3 + ... + a_{k-1}x^{k-1} = f(x)_1

设：

A=a1x+a2x2+a3x3+...+ak−1xk−1A = a_1x + a_2x^2 + a_3x^3 + ... + a_{k-1}x^{k-1}

即：

s1+A=f(x)1s2+A=f(x)2s_1 + A = f(x)_1 \\ s_2 + A = f(x)_2

因此：

s2=f(x)2−f(x)1+s1s_2 = f(x)_2 - f(x)_1 + s_1

如果不知道两个明文的情况下，但两明文有一定的线性关系，也可以算出两明文。

若 s2=as1+bs_2 = as_1 + b，则有同余方程组：

{s1+A−f(x)1 mod p≡0s2+A−f(x)2 mod p≡0as1+b−s2 mod p≡0\left\{ \begin{array}{c} s_1 + A - f(x)_1 \bmod p \equiv 0 \\ s_2 + A - f(x)_2 \bmod p \equiv 0 \\ as_1 + b - s_2 \bmod p \equiv 0 \end{array} \right.

可以用 sage 计算：

1
2
3
4
5
6
7
8
9
10
11
12
a = 99999999
b = 123456

PR.<s1,s2,A> = PolynomialRing(Zmod(p))
f1 = s1 + A - y1
f2 = s2 + A - y2
f3 = a * s1 + b - s2
Fs = [f1, f2, f3]
I = Ideal(Fs)
B = I.groebner_basis()
print 's1 =', ZZ(-B[0](0, 0, 0))
print 's2 =', ZZ(-B[1](0, 0, 0))