昨晚，由徐杨助教给我们带来「助教技术分享会」第 4 讲，题为「浅谈RSA 算法」。

 

前言

本文分为上下两个篇幅，第一部分是基础数据知识，包括了：模运算、最大公约数、辗转相除法、扩展欧几里得算法、模运算的逆元、欧拉函数、欧拉定理等。

第二部分是RSA算法的介绍与代码部分，包括了：算法准备、数据加解密、RSA算法证明、蒙哥马利算法、RSA算法代码与示例等。

希望读者通过阅读此文，能够将RSA算法的来龙去脉搞清楚，同时可以基于RUST语言开发一个简单的RSA算法模拟。

 

基础数学知识点

 

模运算

已知正整数a，b，我们求出a除以b的商与余数：

b * q + r = a

其中q是商，r是余数，也可以记做：

a ≡ r (mod b)，这里的mod即称为模运算，以上表达式读作“a模b等于r”或者“a模b余r”。

这里举几个举例：

9 ≡ 1 (mod 4)

100 ≡ 2 (mod 7)

所有mod b值相同的整数构成的集合，称为一个剩余类。

例如：集合{3n+1}即为mod 3的一个剩余类。对于整数n，其剩余类个数为n个，分别对应mod n余0、1、2...n-1。任何整数a在且仅在mod n的一个剩余类中。

所谓“剩余系”，就是指对于某一个特定的正整数n，一个整数集中的数模n所得的余数域。例如：集合{2, 5, 9}模3的剩余系为：{2, 2, 0}。

从模n的每个剩余类中各取一个数，得到一个由n个数组成的集合，叫做模n的一个完全剩余系。例如，集合{3, 4, 5}为模3的一个完全剩余系。

在与模n互素的全体剩余类中，从每一个类中各任取一个数作为代表组成的集合，叫做模m的一个简化剩余系，或称作缩系。

例如，模5的一个简化剩余系是{1，2，3，4}，模10的一个简化剩余系是{1，3，7，9}，模18的一个简化剩余系是{1，5，7，11，13，17}。

和模运算相关的定理有蛮多的，例如：威尔逊定理，费马小定理，欧拉定理，感兴趣的同学可以移步仔细阅读体会。本文稍后会介绍欧拉定理。

 

最大公约数

最大公约数，是指能够被两个正整数m，n整除的最大的正整数，记做：gcd(m, n)

这里举几个举例：

gcd(2, 6) = 2

gcd(6, 39) = 3

 

辗转相除法

辗转相除法又称欧几里德算法，是指用于计算两个正整数m，n的最大公约数。应用领域有数学和计算机两个方面。

计算公式：

gcd(m,n) = gcd(n, m mod n)

例如，我们要计算gcd(62, 26)

步骤	公式	m mod n
1	gcd(62, 26)	10
2	gcd(26, 10)	6
3	gcd(10, 6)	4
4	gcd(6, 4)	2
5	gcd(4, 2)	0

这样，gcd(62, 26) = 2

具体证明请参考这里。

这里是辗转相除法的RUST实现代码。

pub fn gcd(x: u64, y: u64) -> u64 {

let remainder = x % y;

if remainder == 0 {

return y;

} else {

return gcd(y, remainder);

}

}

#[cfg(test)]

mod tests {

use super::*;

#[test]

fn gcd_works() {

assert_eq!(gcd(2, 4), 2);

assert_eq!(gcd(6, 27), 3);

assert_eq!(gcd(4, 2), 2);

assert_eq!(gcd(27, 6), 3);

}

}

 

扩展欧几里得算法

扩展欧几里得算法（Extended Euclidean Algorithm，下面简称EEA）是欧几里得算法（又叫辗转相除法）的扩展。

该算法基于这个一样假设，对于整数a，b，除了计算其gcd值外，我们是否可以找到整数x、y（其中一个很可能是负数），使得：ax + by = gcd(a, b)。

在欧几里得算法计算过程中，我们仅仅使用到了每次迭代的余数，EEA则使用到了商。

下面我们使用一个例子，来看看具体的一个步骤：

假设：a = 240，b = 46

序号	商 Qi	余数 Ri	Xi	Yi	Sumi
1		240	1	0
2		46	0	1
3	240 / 46 = 5	240 % 46 = 10	1 - 5 * 0 = 1	0 - 5 * 1 = -5	1 * 240 + (-5) * 46 = 10
4	46 / 10 = 4	46 % 10 = 6	0 - 4 * 1 = -4	1- 4 * (-5) = 21	(-4) * 240 + 21 * 46 = 6
5	10 / 6 = 1	10 % 6 = 4	1 - 1 * (-4) = 5	(-5) - 1 * 21 = -26	5 * 240 + (-26) * 46 = 4
6	6 / 4 = 1	6 % 4 = 2	(-4) - 1 * 5 = -9	21 - 1 * (-26) = 47	(-9) * 240 + 47 * 46 = 2
7	4 / 2 = 2	4 % 2 = 0	5 - 2 * (-9) = 23	(-26) - 2 * 47 = -120	23 * 240 + (-120) * 46 = 0

在上面的步骤中，如果只看左侧的三列，这就是辗转相除法的计算步骤。第七步算出余数R7是零，那么这两个数的gcd则为R6 = 2。

这里，在右侧又多了三列，其中Xi与Yi都使用到了商Qi, 其计算公式为：

Xi= Xi-2- Qi * Xi-1

Yi= Yi-2- Qi * Yi-1

Sumi= Xi* a + Yi* b

算法迭代的第一步，我们看第三行，X3 = 1, Y3 = -5, 代入公式算出：Sum3= 10 = R3

根据Xi与Yi的计算公式，不难得出每一步中，Sumi = Ri。具体推演过程中也验证了这个想法。

这样，根据这个迭代的过程，到了第7步，我们算出了gcd是2，那么反推到第6步算出的X6与Y6，即是我们需要的解，使得：(-9) * 240 + 47 * 46 = 2。

EEA的具体证明请参考这里。

好了，说了这么多，这个EEA究竟用来干啥用？？？

 

模运算的逆元

逆元是抽象代数中的一个概念，指在群G中任意一个元素a，都在G中有唯一的逆元a'，具有性质：

a · a' = a' · a = e ( · 为该群中定义的运算)。其中，e为该群的单位元。

将这个概念类推到模运算中，模运算单位元是 1。

那么对于整数n来说，逆元的定义为，对于任意a，存在a'，使得：

a * a' ≡ 1 (mod n)

定义逆元函数：

inv(a, n) = a'

我们举几个例子：

inv(5, 11) = 9

inv(7, 13) = 2

inv(10, 31) = 28

inv(12, 29) = 17

定理：假设p为质数，对于任何正整数a，a < p，那么存在唯一的逆元a'，a' < p，使得：

a * a' ≡ 1 (mod p)

使用反证法可以证明该定理，同学们自己开一下脑洞。

假设已知一个大质数p，例如质数10888869450418352160768000001，我们又知道a，怎么求出a的逆元呢？我们可以穷举，从2开始，一直到p-2，计算a * a' mod p的值，直到等于1为止。

那么，是否还有更好的方法呢？

有，方法之一就是使用EEA算法。

因为p是质数，对于任意a < p，gcd(a, p) = 1

又根据EEA，我们知道可以求得：x和y，使得：a * x + p * y = gcd(a, p) = 1

所以：a * x ≡ 1 (mod p)

所以a的逆元就是x。这里需要注意，有可能此方法算出的x是负数，或者大小超过了p，所以最后需要对x做一个优化。详见代码。

这里是EEA和INV的RUST实现代码。

fn ext_euclid(x: i64, y: i64) -> (i64, i64, i64) {

let (mut old_r, mut r) = (x, y);

let (mut old_s, mut s) = (1, 0);

let (mut old_t, mut t) = (0, 1);

if y == 0 {

return (1, 0, x);

} else {

while r != 0 {

let q = old_r / r;

let t_r = r;

r = old_r - q * r;

old_r = t_r;

let t_s = s;

s = old_s - q * s;

old_s = t_s;

let t_t = t;

t = old_t - q * t;

old_t = t_t;

}

return (old_s, old_t, old_r);

}

}

fn inv(a: i64, p: i64) -> i64 {

let (r, _, _) = ext_euclid(a, p);

return ((r % p) + p) % p

}

#[cfg(test)]

mod tests {

use super::*;

#[test]

fn ext_euclid_works() {

assert_eq!(inv(5, 11), 9);

assert_eq!(inv(7, 13), 2);

assert_eq!(inv(10, 31), 28);

assert_eq!(inv(12, 29), 17);

}

}

 

欧拉函数

对正整数n，欧拉函数是小于或等于n的正整数中与n互质的数的数目（因此φ(1) = 1）。此函数以其首名研究者欧拉命名(Euler's totient function)，它又称为Euler's totient function、φ函数、欧拉商数等。 例如φ(8) = 4，因为1,3,5,7均和8互质。

欧拉函数如下所示：

 

其中，x有n个质因数，pi是第i个质因数。例如我们需要计算φ(8)，8的质因数只有2，那么：

φ(8) = 8 * (1 - 1/2) = 8 * 1/2 = 4，正确

对于任何质数p，φ(p) = p - 1

对于任何质数p和q，φ(p * q) = (p - 1) * (q - 1)

 

欧拉定理

欧拉定理（Euler Theorem，也称费马-欧拉定理或欧拉函数定理）是一个关于同余的性质。欧拉定理得名于瑞士数学家莱昂哈德·欧拉，该定理被认为是数学世界中最美妙的定理之一。欧拉定理实际上是费马小定理的推广。

欧拉定理表明，若n,a为正整数，且n,a互质，则：

 

 

其中，φ(n)是欧拉函数。

这里做一个简洁的证明：

将1-n中，与n互质的数，定义为集合A：{X1, X2 ... , Xφ(n)}，可知总共有φ(n)个数。同时它们都与n互质。所以集合A是mod n的一个缩系。

构造集合B：{Y1 = a * X1, Y2 = a * X2 ... , Yφ(n)= a * Xφ(n)}，因为a、Xi都与n互质，可知Yi也都与n互质。

这样，集合B是否也是一个缩系？假设不是一个缩系，那么必存在i、j，使得：Yi≡ Yj(mod n)

可得：a * Xi- a * Xj≡ 0 (mod n) => a * (Xi- Xj) ≡ 0 (mod n)

由于a与n互质，可知：Xi - Xj ≡ 0 (mod n)，显然，这与Xi的定义相矛盾，Xi、Xj是小于n的不同整数。

所以，集合B也是一个缩系。

这样，我们有两个mod n的缩系集合A、B，根据缩系的定义，将两个集合的各自元素求其积，可知这两个积mod n是相等的：

(a * X1) * (a * X2) * (a * Xφ(n)) ≡ X1* X2 * Xφ(n)mod n

=> aφ(n)* (X1X2...Xφ(n)) ≡ X1X2...Xφ(n)mod n

=> (aφ(n)- 1) * (X1X2...Xφ(n)) ≡ 0 mod n

由于Xi与n互质，可得：aφ(n)- 1 ≡ 0 mod n => aφ(n)≡ 1 mod n，定理得证。

 

RSA算法

 

综述

RSA算法是由Ron Rivest、Ad Shamir、Leonard Adleman，于1977年提出的，RSA就是他们三位名字的首字母拼在一起组成的。

RSA是一种非对称加密算法，其核心思路是，alice公开自己的公钥Kpub，同时仅仅只有alice知道自己的私钥Kpri；bob使用Kpub对数据m进行加密，得到秘文c，并将c传输给alice。

此时，alice使用私钥Kpri对密文c进行解密，得到明文数据m。由于只有alice知道私钥Kpri，就算其他人得到了密文c，也无法解密出其明文来。对大整数进行因数分解的难度，是RSA算法安全性的保证。

 

算法准备

1. 选取两个质数p、q，p != q，计算其乘积：N = p * q

2. 计算：φ(N) = (p - 1) * (q - 1)

3. 让：r = φ(N)

4. 任选一整数e，1 < e < r，且：gcd(e, r) = 1

5. 根据算法EEA，求得e的逆元d，使得：ed ≡ 1 (mod r)

6. 最终，(N, e）即是公钥，（N, d）即是私钥

7. alice公开公钥（N，e），同时将私钥（N，d）自己私密保管

 

数据加密

bob对于任何消息m，将其转化成整数n，计算：c ≡ ne (mod N)，这里c即是密文。

 

数据解密

alice拿到密文，计算：n ≡ cd (mod N)，这里n即是明文。

 

算法证明

因为：ed ≡ 1 (mod r)

所以：ed = rm + 1，m为正整数

所以：cd≡ ned≡ nrm + 1≡ n * nrm(mod N)

根据欧拉定理，可知：nr ≡ 1 (mod N)，代入上式：

n * 1m≡ n (mod N)

所以结论为：cd≡ n (mod N)

 

蒙哥马利算法

在RSA算法的实现过程中，首先需要随机生成两个大质数，目前推荐至少是2048位的质数，这样可以保证安全。我们在程序模拟过程中，不会使用这么大的质数，仅仅做一个演示。

RSA需要的几个外部算法，gcd与EEA，在之前的部分都已经介绍过。可以注意到，其加密、解密的过程，牵涉到模运算与指数计算，因为真实环境下，e或者d这两个值，可能会非常的大（例如超过21000），所以这里还需要有一些优化，不能直接计算ne、cd。

所以，这里还需要介绍一个算法 - 蒙哥马利算法，蒙哥马利(Montgomery)幂模运算是快速计算ab mod p的一种算法，是RSA加密算法的核心之一。需要注意的是，该算法仅仅是作为提升效率而使用，与RSA算法的安全性本身没有关系。

我们可以看一个简单的例子，例如需要计算A13 mod P，我们可以这样：

步骤	公式	值	操作	op_code
1	A1 = 1 * 1 mod p	1	二次方	1
2	A2 = A1 * A mod p	A1	乘A	-
3	A3 = A2 * A2 mod p	A2	二次方	1
4	A4 = A3 * A mod p	A3	乘A	-
5	A5 = A4 * A4 mod p	A6	二次方	0
6	A6 = A5 * A5 mod p	A12	二次方	1
7	A7 = A6 * A mod p	A13	乘A	-

我们将13按照2进制显示为：0b1101，大家可以观察到，这里的四个数字，1101，和上表中的op_code的1101相对应，其意思是：

• 首先设置结果的初始值为1

• 将幂值的二进制表达式从左向右做循环

  • 遇到1，则做两个操作：二次方操作 + 乘A操作

  • 遇到0，则仅仅做：二次方操作

二次方操作，其实就是自乘，所以其核心思路是，将模幂运算降为了模乘运算，算法难度从O(2n)降为了O(n)。

 

RSA算法代码与示例

至此所有的准备工作都好了，我们可以愉快的编码了。

 

初始化

RSA算法初始化的部分如下。

extern crate cryptography_algo;

use cryptography_algo::gcd;

use cryptography_algo::ext_euclid;

use rand::Rng;

use std::env;

fn main() {

let p;

let q;

if env::args().len() == 3 {

p = i64::from_str_radix(&env::args().nth(1).unwrap(), 10).unwrap();

q = i64::from_str_radix(&env::args().nth(2).unwrap(), 10).unwrap();

if !prime_tools::is_u64_prime(p as u64) {

panic!("p is not a prime number!");

}

if !prime_tools::is_u64_prime(q as u64) {

panic!("q is not a prime number!");

}

} else {

p = 2134324421;

q = 1990843139;

}

let N: i64 = p * q;

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

println!("p: {} ,q: {} N: {}, r: {}", p, q, N, r);

let mut rng = rand::thread_rng();

let mut e;

loop {

e = rng.gen::<u16>() as i64;

if gcd::gcd(e, r) == 1 {

break;

}

}

let d = ext_euclid::inv(e, r);

println!("e: {}, d: {}", e, d);

println!("public key is(N, e): ({}, {})", N, e);

println!("private key is(N, d): ({}, {})", N, d);

}

核心部分很简单，根据p、q，计算N、r；随机生成e，需要gcd(e, r) = 1，同时求得e的逆元d。

 

加密解密

下面的是蒙哥马利算法的RUST代码。

pub fn power_mod(base: i64, mut power: i64, N: i64) -> i64 {

let mut bits = Vec::new();

while power != 0 {

match power & 1 {

1 => bits.push(true),

0 => bits.push(false),

_ => {}

}

power = power >> 1;

}

let mut result: i64 = 1;

while let Some(bit) = bits.pop() {

result = mod_multiply(result, result, N);

if bit {

result = mod_multiply(result, base, N);

}

}

result

}

pub fn mod_multiply(a: i64, b:i64, N: i64) -> i64 {

((a as i128 * b as i128) % (N as i128)) as i64

}

#[cfg(test)]

mod tests {

use super::*;

#[test]

fn mod_multiply_works() {

assert_eq!(mod_multiply(10, 20, 17), 13);

assert_eq!(mod_multiply(6, 27, 5), 2);

assert_eq!(mod_multiply(4, 2, 19), 8);

assert_eq!(mod_multiply(27, 6, 13), 6);

}

#[test]

fn power_mod_works() {

assert_eq!(power_mod(3, 4, 5), 1);

assert_eq!(power_mod(6, 10, 5), 1);

assert_eq!(power_mod(5, 13, 19), 17);

assert_eq!(power_mod(27, 6, 51), 9);

assert_eq!(power_mod(45, 13, 19), 7);

assert_eq!(power_mod(1234567, 100, 199), 29);

assert_eq!(power_mod(66887799, 1000000, 1001), 22);

assert_eq!(power_mod(1357924680, 999999, 666889), 355775);

assert_eq!(power_mod(113355778866, 99999, 981287), 797582);

assert_eq!(power_mod(113355778866, 999999, 5050404053), 31958690);

}

}

加解密算法入口程序如下。

extern crate cryptography_algo;

use cryptography_algo::power_mod;

use std::env;

fn main() {

if env::args().len() < 3 {

println!("program should be start with: rsa_enc_dec base power modula");

return;

}

let base = i64::from_str_radix(&env::args().nth(1).unwrap(), 10).unwrap();

let power = i64::from_str_radix(&env::args().nth(2).unwrap(), 10).unwrap();

let modula = i64::from_str_radix(&env::args().nth(3).unwrap(), 10).unwrap();

println!("base: {}, power: {}, modula: {}", base, power, modula);

println!("result: {}", power_mod::power_mod(base, power, modula));

}

 

运行示例

可以使用openssl命令得到质数，命令如下：

openssl prime -generate -bits 32 -safe

算法初始化：

target/release/rsa_init 2134324421 1990843139

p: 2134324421 ,q: 1990843139 N: 4249105129947997519, r: 4249105125822829960

e: 1343, d: 1920481616808978247

public key is(N, e): (4249105129947997519, 1343)

private key is(N, d): (4249105129947997519, 1920481616808978247)

初始化需要两个质数参数p、q（不输入则使用默认值），这里注意，由于程序限制p、q的值需要小于：sqrt(263)，p、q最大为3037000493。

初始化的结果是，得到了公钥（N，e），与私钥（N，d）。这时，alice将公钥给到bob，私钥自己私密保管。

数据加密：

target/release/rsa_enc_dec 50412164937805327 1343 4249105129947997519

base: 50412164937805327, power: 1343, modula: 4249105129947997519

result: 1099769683952491905

bob拿到alice的公钥后，对自己的数据，进行加密计算。

该程序有三个参数，第一个是数据明文[50412164937805327]，后面两个参数分别为公钥中的e与N。

加密过程使用蒙哥马利算法进行模幂计算，最终得到了加密后的密文。

这时，bob将密文[1099769683952491905]给到alice。

数据解密：

target/release/rsa_enc_dec 1099769683952491905 1920481616808978247 4249105129947997519

base: 1099769683952491905, power: 1920481616808978247, modula: 4249105129947997519

result: 50412164937805327

alice拿到bob的密文数据后，结合自己的私钥数据，进行数据解密。

该程序有三个参数，第一个是密文数据[1099769683952491905]，后面两个参数分别为私钥中的d与N。

解密过程使用蒙哥马利算法进行模幂计算，最终得到了解密后的明文[50412164937805327]。

可以看到，最终解密出的明文数据，与最初bob生成的数据是一致的。

Bingo!