用Rust语言开发RSA算法模拟

原文： https://mp.weixin.qq.com/s/0BkRrfoP8Q5WJ_1lGXgB1A

 

作者：徐杨，《 Substrate快速入门与开发实战》开发课第三期助教。NBLTrust  Co-Founder&CTO，上海交通大学计算机科学学士和硕士学位，擅长技术管理和程序开发。

 

每周日晚 8 点，都会进行《Substrate 快速入门与开发实战》开发课的内容知识拓展——助教技术分享会。昨晚，由徐杨助教给我们带来「助教技术分享会」第 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!