原文: 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互质。
欧拉函数如下所示:
用Rust语言开发RSA算法模拟-Rebase     
其中,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互质,则:       
用Rust语言开发RSA算法模拟-Rebase   
其中,φ(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!