Python 与 RSA 算法

本周，Computer Security 这门课的作业是用代码实现 RSA 算法，要求加密 + 解密都要实现，语言不限。

这种东西网上一搜一大把，但因为 RSA 太重要了，老师额外要求要 “Defend your code!” (相当于当着他面，一对一做一次 Code Review)。没办法，老老实实的把相关课件重新学了一遍，然后动笔实现。害怕到时候讲不清楚，所以特地写一篇博客理一遍思路。

虽说语言不限，但能方便快捷的操作这么大的数 (老师要求使用至少 1024 bit 的质数！) 的语言还真不多，对比了一下，最终选择了 Python。

快速模幂运算

快速模幂运算，是与 RSA 相关的算法中最基础也是最核心的一个算法了。毕竟 RSA 加密和解密的本质就是做一个模幂运算；而求我们生成一个大质数的过程中也需要用到这个函数 (费马定理，验证这个数是否为质数)。

这个算法主要是为了减少计算机的计算量，毕竟要硬算 $ 1.047247e152^{5.80901e149} \pmod{65537} $ 这样的大数的模幂可不容易。

我们都知道，关于取模，有这么一个公式：

$$ c \bmod m = (a \cdot b) \bmod m = [(a \bmod m) \cdot (b \bmod m)] \bmod m $$

再结合 $ a^{b + c} = a^{b} \cdot a^{c} $，我们可以很容易将我们的模幂运算分解：

$$ a^e \bmod m = a^{b + c} \bmod m = [(a^{b} \bmod m) \cdot (a^{c} \bmod m)] \bmod m $$

但是，如何分解 $ e = b + c $ 最高效又成了问题。这里我们使用平方求幂原理。

首先将 $ e $ 表示成二进制，即：

$$ e = \sum_{i = 0}^{n - 1} a_i 2^i, (a_{n - 1} = 1) $$

$ b^{e} $ 可以表现为：

$$ b^{e} = b^{(\sum_{i = 0}^{n - 1} a_i 2^i)} = \prod_{i = 0}^{n - 1} (b^{2^{i}})^{a_i} $$

因此 $ e $ 为：

$$ e \equiv \prod_{i = 0}^{n - 1} (b^{2^{i}})^{a_i} \pmod m $$

具体实现的代码如下 (递归)：

1
2
3
4
def powmod(p, q, m):
    if not q: return 1
    if q & 1: return (p * pow_mod(p, q - 1, m)) % m
    return pow_mod((p ** 2) % n, q >> 1, m)

虽说递归的实现最为直观，但我在后面的实际使用过程中遇到了性能问题，不得已重新用 while 重构了一遍：

1
2
3
4
5
6
7
8
def powmod(p, q, m):
    r = 1
    while q:
        if q & 1:
            r = (r * p) % m
        q >>= 1
        p = (p * p) % m
    return r

生成一个大大的 (伪) 质数

接着，我们需要一个函数用来生成一个固定位数的大数，这个数可能是质数也可能是合数。

我们直接用 Python 的 random 模块来生成一个 w 位的二进制字符串 (仅有 0 和 1 组成的字符串)，然后将字符串转化为数字 (十进制)：

1
2
3
4
5
6
7
8
9
import random

def fake_prime(w):
    list = []
    list.append('1')
    for i in range(w - 2):
        list.append(random.choice(['0', '1']))
    list.append('1')
    return int(''.join(list), 2)

需要注意的是，二进制字符串的首位必须位 '1'，因为如果首位为 0，就不能保证转换后的数字是 w 位了。同时，末尾也必须为 '1'，因为我们的这个函数的目的是生成 (伪) 质数，而所有的偶数都是合数。

米勒-拉宾素性检验

接下来我们需要一个函数来对生成的 (伪) 质数进行验证，这里我们需要用的米勒-拉宾素性检验 (Miller-Rabin prime test)。

这种算法基于大名鼎鼎的黎曼猜想，其利用了费马小定理的结论：

对于素数 $ n $ 和任意整数 $ a $，有

$$ a^{n − 1} \equiv 1 \pmod n $$

米勒-拉宾素性检验在费马定理的基础上增加了推论：

对于任何大于 $ 2 $ 的素数 $ p $，不存在 $ 1 \bmod p $ 的“非平凡平方根”¹。

让我们假设 $ a $ 是 $ 1 \bmod p $ 的平方根之一 ($ 0 < a < p $)，于是有：

$$ a^{2} \equiv 1 \pmod p \newline (a - 1)(a + 1) \equiv 0 \pmod p $$

因为 $ p $ 是质数，$ (a - 1) \bmod p $ 和 $ (a + 1) \bmod p $ 均不可能为 $ 0 $。

那么唯一的可能性只能是：

$$ (a - 1)(a + 1) = 0 \space or \space (a - 1)(a + 1) = p $$

所以：

$$ a = 1 \space or \space a = p - 1 $$

我们获得推论：如果 $ p $ 是素数，则 $ a^{2} \equiv 1 \pmod p $ 的解为 $ a = 1 $ 或 $ a = p−1 $。

但是对于所有数都使用米勒-拉宾素性检验有点不划算，因为这种算法的开销其实还是挺大的。

我们可以先将 $ 1 $ 至 $ 100 $ 以内的所有素数列出来，用这些素数对需要测试的数取模，其结果如果为 $ 0 $ ，则必并是合数 (素数的倍数必定是合数)；然后再对非零的情况使用米勒-拉宾素性检验：

 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
def prime_miller_rabin(a, p):
    prime_numbers = (2,3,5,7,11,13,17,19,23,29,31,37,41,
                     43,47,53,59,61,67,71,73,79,83,89,97)
    for y in prime_numbers:
        if p % y == 0:
            return False
    
    def fermat():
        return powmod(a, p - 1, p) == 1

    def miller_rabin():
        d = p - 1
        q = 0
        while not(d & 1):
            q = q+1
            d >>= 1
        m = d

        for i in range(q):
            u = m * (2**i)  
            tmp = powmod(a, u, p)
            if tmp == 1 or tmp == p - 1:
                return True
        return False

    return fermat() and miller_rabin()

需要注意的是，米勒-拉宾素性检验是个证明素数的只是一个必要但不充分条件，有极小的概率误判；不过当我们取多个 a 都能通过测试的话，我们可以将误判的概念忽略不计。

因此用一个函数 prime_test 将 prime_miller_rabin 包装一下，生成多个随机数 a，测试其是否都符合条件。

1
2
3
4
5
6
import random

def prime_test(n, k):
    if k <= 0: return True
    a = random.randint(2, n - 1)
    return prime_miller_rabin(a, n) and prime_test(n, k - 1)

接着用 prime_test 来测试 fake_prime 生成的数是否是质数，如果不是，则重新生成，直到获得一个质数为止。

1
2
3
4
5
def get_prime(bit):
    prime_number = fake_prime(bit)
    if prime_test(prime_number, 5):
        return prime_number
    return get_prime(bit)

实现 RSA 加密

终于可以正式开始实现 RSA 算法了！

我们先来看一下公式：

$$ C = M^{e} \bmod n $$

$ M $ 是需要加密的信息。$ (n, e) $ 是我们的公钥。其中 $ n $ 就是 $ p $ 与 $ q $ 的乘积，即：

$$ n = p \cdot q $$

而 $ e $ 是一个小于 $ r $ 并且与 $ r $ 互质的正整数。

$ r $ 就是所谓的欧拉函数，公式如下：

$$ r = \phi(n) = \phi(p) \cdot \phi(q) = (p - 1) \cdot (q - 1) $$

$ e $ 的大小直接影响了加密的速度，所以一般不会很大，也不会很小，并且里面 (二进制) 的 1 要尽可能的少。这里我们选择 0b1000000000001 (即 4097) 作为默认的 $ e $。

这个 $ r $ 虽然不直接参与 RSA 加密，但其必须 $ e $ 互质；我们通过欧几里得算法求 $ e $ 与 $ r $ 的最大公约数，进而判断它们是否互质；如果不互质，就换一个 $ e $ (e -= 1)：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
def get_r_e(p, q, e = 4097):
    def gcd(a, b):
        if not b: return a
        return gcd(b, a % b)

    r = (p - 1) * (q - 1)

    if gcd(e, r) == 1:
        return (r, e)
    return gen_orla(p, q, e - 1)

套用公式，实现对 $ M $ 的加密：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
if __name__ == "__main__":

    p = get_prime(1024)
    q = get_prime(1024)
    n = p * q

    r, e = get_r_e(p, q)

    M = int("7355608")

    C = powmod(M, e, n)

实现 RSA 解密

除了加密，我们还需要实现解密。

同样的，先来看一看公式：

$$ M = C^{d} \bmod n $$

这里的 $ C $ 是上一步得到的密文。$ (n, d) $ 是我们的私钥。

$ n $ 已经解释过了，是 $ p $ 与 $ q $ 的乘积 ($ n = p \cdot q $)。

这个 $ d $ 是 $ e $ 关于 $ r $ 的逆模元，即：

$$ e \cdot d \equiv 1 \pmod r $$

所以：

$$ d = e^{-1} \bmod r $$

拓展欧几里得算法

这里我们使用拓展欧几里得算法：

令 $ s_0 = e $，$ s_1 = r $，$ x_0 = 1 $，$ x_1 = 1 $，$ y_0 = 0 $，$ y_1 = 1 $，

当 $ s_{i} > 0 $ 时：

$$ x_{i} = \lfloor \dfrac{s_{i-2}}{s_{i-1}} \rfloor \cdot x_{i-1} + x_{i-2} \newline y_{i} = \lfloor \dfrac{s_{i-2}}{s_{i-1}} \rfloor \cdot y_{i-1} + y_{i-2} \newline s_{i} = s_{i-2} \bmod s_{i-1} $$

当 $ s_{n-1} = 0 $ 时：

$$ (e, r) = s_{n} \newline x = (-1)^{n} \cdot x_{n} \newline y = (-1)^{n+1} \cdot y_{n} $$

写成代码就是：

1
2
3
4
5
6
7
def mod_inverse(e, r):
    x0, x1, x2 = r, 0, 1
    while r != 0:
        x1, x2 = x2 - ((e // r) * x1), x1
        e, r = r, e % r
    if x2 < 0: x2 += x0
    return x2

然后套用公式，就可以对 $ C $ 进行解密了：

1
2
3
4
5
if __name__ == "__main__":

    # p, q, C, r, n 来自于上一步
    d = mod_inverse(e, r)
    D = powmod(C, d, n)

到现在为止，虽然 RSA 算法已经可以用了，但仍然很慢；特别是解密，当质数的位数非常大时就非常容易卡死。

因此，我们需要使用中国剩余算法对解密进行优化。

中国剩余算法

中国剩余算法来自于《孙子算经》的第 26 题：

今有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二，问物几何？
答曰：二十三。
术曰：三三数之剩二，置一百四十；五五数之剩三，置六十三，七七数之剩二，置三十，并之。得二百三十三，以二百一十减之，即得。凡三三数之剩一，则置七十；五> 五数之剩一，则置二十一；七七数之剩一，则置十五；一百六以上以一百五减之即得。

其核心思想就是辗转相除法。

对于给定的数 $ k $，令 $ n_{i} (1 \leq i \leq k) $ 与其互质；令 $ a_{i} (1 \leq i \leq k) $ 满足 $ 0 \leq a_{i} \leq n_{i} $。此时，有且只有一个数 $ x $ 同时满足：

$$ 0 \leq x < N (N = n_{1} \cdot n_{2} \cdot \mathellipsis \cdot n_{k}) $$

和

$$ x \equiv a_{1} \bmod{n_{1}} \newline x \equiv a_{2} \bmod{n_{2}} \newline \vdots \newline x \equiv a_{k} \bmod{n_{k}} $$

其在 RSA 解密中的应用如下：

计算 $ d_{p} $ 和 $ d_{q} $：
$$ d_{p} \equiv d \bmod{p - 1} \newline d_{q} \equiv d \bmod{q - 1} $$
计算 $ m_{p} $ 和 $ m_{q} $：
$$ m_{p} \equiv c^{d_{p}} \bmod{p} \newline m_{q} \equiv c^{d_{q}} \bmod{p} $$
使用 $ m $ 将其改写为同余的形式：
$$ m \equiv m_{p} \bmod{p} \newline m \equiv m_{q} \bmod{q} $$
使用中国剩余定理计算 $ m $：
$$ 1 = y_{p}p + y_{q}q $$
因为 $ n = p \cdot q $，所以：
$$ m \equiv (m_{p} y_{q} q + m_{q} y_{p} p) \bmod{n} $$

写成代码就是：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
def chinese_remainder(c, d, p, q):
    dp = powmod(d, 1, p - 1)
    dq = powmod(d, 1, q - 1)

    m1 = powmod(c, dp, p)
    m2 = powmod(c, dq, q)

    qinv = mod_inverse(q, p)
    h = (qinv * (m1 - m2)) % p
    m = m2 + h * q
    return m

然后使用 chinese_remainder 代替 powmod 进行解密：

1
2
3
4
5
if __name__ == "__main__":

    # p, q, C, r, n 来自于上一步
    d = mod_inverse(e, r)
    D = chinese_remainder(C, d, p, q)

用 Java 实现

传送门 -> RSA 加解密的 Java 实现。

References

附录

完整的代码如下：

  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
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#!/usr/bin/env python

import random

def powmod(p, q, m):
    """
    p^q mod n
    """
    r = 1
    while q:
        if q & 1:
            r = (r * p) % m
        q >>= 1
        p = (p * p) % m
    return r

def fake_prime(w):
    list = []
    list.append('1')
    for i in range(w - 2):
        list.append(random.choice(['0', '1']))
    list.append('1')
    return int(''.join(list), 2)

def prime_miller_rabin(a, p):
    prime_numbers = (2,3,5,7,11,13,17,19,23,29,31,37,41,
                     43,47,53,59,61,67,71,73,79,83,89,97)
    for y in prime_numbers:
        if p % y == 0:
            return False
    
    def fermat():
        """
        a^(p-1) = 1 (mod p)
        """
        return powmod(a, p - 1, p) == 1

    def miller_rabin():
        """
        a^2 = 1 (mod p)
        a = 1 or a = p - 1
        """
        d = p - 1
        q = 0
        while not (d & 1):
            q = q+1
            d >>= 1
        m = d

        for i in range(q):
            u = m * (2**i)  
            tmp = powmod(a, u, p)
            if tmp == 1 or tmp == p - 1:
                return True
        return False

    return fermat() and miller_rabin()


def prime_test(n, k):
    if k <= 0: return True
    a = random.randint(2, n - 1)
    return prime_miller_rabin(a, n) and prime_test(n, k - 1)

def get_prime(bit):
    prime_number = fake_prime(bit)
    if prime_test(prime_number, 5):
        return prime_number
    return get_prime(bit)

def get_r_e(p, q, e = 4097):
    def gcd(a, b):
        """
        Euclidean algorithm
        """
        if not b: return a
        return gcd(b, a % b)

    def euler(p, q):
        """
        Euler's totient function
        """
        return (p - 1) * (q - 1)

    r = euler(p, q)

    if gcd(e, r) == 1:
        return (r, e)
    return get_r_e(p, q, e - 1)

def mod_inverse(e, r):
    """
    Modular inverse operation\n
    y1 = e^-1 mod r
    """

    x0, x1, x2 = r, 0, 1
    while r != 0:
        x1, x2 = x2 - ((e // r) * x1), x1
        e, r = r, e % r
    if x2 < 0: x2 += x0
    return x2

def chinese_remainder(c, d, p, q):
    """
    Chinese remainder theorem
    """
    dp = powmod(d, 1, p - 1)
    dq = powmod(d, 1, q - 1)

    m1 = powmod(c, dp, p)
    m2 = powmod(c, dq, q)

    qinv = mod_inverse(q, p)
    h = (qinv * (m1 - m2)) % p
    m = m2 + h * q
    return m

if __name__ == "__main__":

    p = get_prime(1024)
    q = get_prime(1024)
    n = p * q

    r, e = get_r_e(p, q)
    d = mod_inverse(e, r)

    print('Private key:\n')
    print('p: %d\n' % p)
    print('q: %d\n' % q)
    print('d: %d\n\n' % d)

    print('Public key:\n')
    print('n: %d\n' % n)
    print('e: %d\n\n' % e)

    M = int(input("Input message: "))

    C = powmod(M, e, n) # encryption
    print('\nAfter encryption is completed, the ciphertext obtained:\n%d\n' % C)

    D = chinese_remainder(C, d, p, q) # decryption
    print('Decryption is complete, and the plaintext obtained is:\n%d\n' % D)

我们发现 $ 1^{2} \bmod p $ 和 $(-1)^{2} \bmod p $ 总是得到 $ 1 $，我们称 $ −1 $ 和 $ 1 $ 是 $ 1 \bmod p $ 的“平凡平方根” ↩︎