RSA算法

RSA是一种非对称的算法，也就是加密所用的私钥和解密的公钥不同。给定明文 $x$，公钥为 $m$ 和 $e$，密文为 $y$，那么加密过程用函数 $\text{encrypt}(x)$表示：

对于公钥 $m$，用于解密的私钥为 $d$，那么解密过程用函数 $\text{decrypt}(x)$表示：

加密和解密过程使用同样的运算，上边两个式子使用了同余关系，而在实际应用中应该是模运算，也就是求一个整数模 $m$ 最小非负剩余。

数学原理

模数 $m$ 和指数 $e$ 是为公钥；而指数 $d$ 是私钥，由公钥计算出来。一般来说指数 $e$ 是固定的，也可以是从多个备选值中选择（比如？）；$m$ 是两个长度相当的质数的乘积 $m=pq$。$m$ 和 $e$ 存在以下关系

其中，$\varphi(m)$ 是欧拉函数，是对 $[1,m]$ 之间和 $m$ 互质的整数的计数；如果 $p$ 和 $q$ 互质，那么 $\varphi(pq)=\varphi(p)\varphi(q)=(p-1)(q-1)$。关于欧拉函数和欧拉定理更多细节，可以参考另一篇文章。

私钥 $d$

私钥 $d$， 和 $e$ 以及 $m$ 存在以下关系

这是贝祖定理的结果：如果两个整数 $a$ 和 $b$ 互质，那么存在两个整数 $s$ 和 $t$ 使得 $as+bt=1$；转换一下就是 $as’=1+bt’$。在这里，$a=e$ 就是公钥之一的指数，而 $s’=d$ 是私钥；为了计算出 $s$ 的值，需要使用扩展欧几里得算法。

加解密证明

对密文 $y$ 进行解密也就是计算 $y$ 的 $d$ 次幂模 $m$ 剩余。代入贝祖定理的结果 $ed=1+k\varphi(m)$

如果 $\gcd(x, m)=1$，这个时候可以应用欧拉定理的结果 $x^{\varphi(m)}\equiv1\pmod m$

加解密的正确性成立。

如果 $\gcd(x, m)\ne 1$，那么 $\gcd(x,m)$ 的取值只有三种情况：$p$，$q$或者$m=pq$。假如 $\gcd(x, m)=pq=m$，这表示密文 $x$ 是 $m$ 的倍数，$x\equiv0\pmod m$。那么

假如 $\gcd(x, m)=p$ 或者 $q$，虽然欧拉定理不适用，但是仍然可以得出一样的结果

因此 $y^d\equiv x^{k\varphi(m)+1}\equiv x\pmod m$，加解密的正确性成立。

快速幂取余运算

RSA算法的加解密过程都需要对一个整数计算高次幂模 $m$ 剩余。$d$ 可能是一个二进制长度超过2000位的整数，使用循环进行乘法计算是不现实的。有一种名为快速幂的高效算法可以计算高次幂。快速幂(repeated squaring,或者叫successive squaring)的主要思路利用幂指数的二进制表示形式，将 $d$ 次乘法转化为 $\log_{2}{d}$ 个指数为 $2^i$ 的幂的乘法，能够将时间复杂度从 $O(n)$ 下降至 $O(\log{n})$。

该算法首先检查幂指数的二进制表示形式。举个例子，计算 $5^{11} \mod 7$

初始值为 $1$ 。从 LSB开始执行循环，对于第 $i$ 个bit，这个 bit 如果为1，就表示需要对上一轮的结果乘以 $5^{2^i}$；而由于 $5^{2^i}=5^{2^{i-1}\cdot2}=5^{2^{i-1}}\cdot 5^{2^{i-1}}$，也就是说 $5^{2^i}$ 可以根据上一次循环 $5^{2^{i-1}}$ 计算出来。因为每轮循环需要额外1次乘法，有 $\log_{2}{n}$ 次循环，所以需要 $2\log_{2}{n}$ 次乘法。

python实现代码如下：

def func(base, exp):
    if exp < 0:
        ret = 1/base
        exp = -1 * exp
    else:
        ret = base

    while exp:
        if exp & 1:
            ret = (ret * base)
        exp >>= 1
        base = (base * base)
    return ret

然而这个实现并不能通过 leetcode，原因是随着 exp 增大、乘积base也在增加，需要想办法优化这个变量。同样是从指数的二进制形式入手，从MSB到第 $i$ 个bit的幂，是MSB到$i-1$个bit的幂的平方和底数的乘积，例如

而MSB到$i-1$的幂又可以如法炮制。按照这个思路可以写出优雅的递归代码：

def bar(base, exp):
#def bar(base, exp, m):
    if exp == 0:
        return 1 if base else 0
    elif exp < 0:
        base = 1 / base
        exp = -1 * exp
    ret = bar(base, exp >> 1)
    #ret = bar(base, exp >> 1, m)
    if exp & 1:
        return ret * ret * base # % m # 取余需要加上 % m
    else:
        return ret * ret # % m # 取余需要加上 % m

这个代码计算的是高次幂。如果要计算幂的模$m$剩余就要按照注释去修改。非递归的循环版本如下：

def bar(base, exp):
#def bar(base, exp, m):
    ret = 1
    if exp < 0:
        base = 1.0 / base
        exp = -1 * exp
    for bit in bin(exp)[2:]:
        ret *= ret # % m
        if int(bit):
            ret *= base # % m
        #ret %= m
    return ret

安全性

公钥 $e$ 和 $m$ 对所有人公开，所有人都可以获取到 $e$ 和 $m$ 的值；私钥 $d$ 保持隐匿，其值由欧几里得算法根据 $e$ 和 $\phi(m)$ 计算出来，而 $\varphi(m)$ 和 $m$ 的整数分解方式有关。RSA算法的关键在于，如果其他人获得了 $m$ 但不知道如何分解 $m$，就无法在短时间里计算出欧拉函数的值、也就是私钥。

一个简单的、用来分解整数的方式就是测试从 $2$ 到 $\sqrt{m}$ 的所有质数因子的暴力算法。如果 $m$ 是一个很小的整数，比如说32位的整数，即使是一台1C2G的服务器也可以在几毫秒内用暴力程序计算出所有因子。但是对于超过2000位的整数来说，使用暴力算法在一个可接受的时间内计算出结果几乎是不可能的。对于一个比特长度为 $n$ 的整数变量，它的可能取值范围是$[0, 2^n-1]$，暴力运算的最差情况是 $O(\sqrt{2^n})=O(2^{\frac{n}{2}})$，属于指数级别的复杂度。

在特定的场景下，Pollard $p-1$ 算法可以较快地分解给定的整数。对于一个合数 $N$ 和它的一个质数因子 $p$，如果 $p-1$ 的因子包含一些比较小的质数 $q$，那么 Pollard 算法可以在较短时间内计算出 $p$。

先提一个概念：

$k$-smooth

如果一个整数的最大质数因子不超过 $k$，那么可以说这个整数是 $k$-smooth 的。比如，

• $15=3\times5$ 是 $5$-smooth的，不是 $3$-smooth 的，但同时也是 $7$-smooth的，这是因为 $3\lt5\lt7$；
• $180=2^2\cdot3^2\cdot5$，$180$ 也是 $5$-smooth的。

Pollard $p-1$ 分解算法

对于一个底数 $a$，由费马小定理可以知道 $a^{p-1}\equiv1\pmod p$，换句话说 $p$ 整除 $a^{p-1}-1$。加上 $p$ 整除 $N$，也就说 $p$ 是 $a^{p-1}-1$ 和 $N$ 的公因数，这个时候运用辗转相除法就可以计算出 $N$ 的一个因子。不过到这里为止因为不清楚 $p$，也就无法获得 $a^{p-1}-1$。

但假如 $B$ 是 $p-1$ 的倍数，费马小定理也是成立的。Pollard 算法的主要思路是，首先假定 $N$ 的一个质数因子 $p$ 而且 $p-1$ 是 $k$-smooth 的，然后寻找一个足够大的$p-1$的倍数 $B$ ，然后尝试通过辗转相除法计算出质数因子 $p$。

假设存在质数因子 $p$ 且 $p-1$ 是 $k$-smooth

对于所有质数 $q<k$，令 $E_q=q^{t_q} \le N\le q^{t_q+1}$ ，构造出

关键在于 $p-1$ 必然可以整除 $B$。这是因为，基于 $p-1$ 是 $k$-smooth 的假设，根据算术基本定理

考虑 $p-1$ 的质数因子 $q$，$q^{s_q}\mid {p-1}$ 且

因此 $q^{s_q}\lt q^{t_q+1}$。这两个数都是 $q$ 的幂，所以 $s_q\lt t_q+1$，进一步地，$q^{s_q}\le q^{t_q}$。符号 $p$ 的任意性说明，$p-1$的任何一个质数因子的幂都要小于 $B$ 的对应的质数因子的幂，换句话说 $p-1$ 整除 $B$。接着就可以用辗转相除法计算 $\gcd(2^{B}-1, N)$。

这个时候 $2^{B}-1$ 仍然是一个非常大的数字，可以用前面提到的快速幂模$N$取余算法。举个例子，如果需要分解 $N=221=13*17$，假设 $N$ 有一个质数因子是 $2$-smooth 的，那么：

计算 $\gcd(2^B-1, 221)$ ：

bar(2, 2**7, 221)
# 35
gcd(bar(2, 2**7, 221)-1, 221)
# 17

因此 $17$ 是 $221$ 的一个因子，是一个质数，同时 $16=17-1$ 也是 $2$-smooth 的。

假设所有质数因子 $p -> p-1$ 都不满足 $k$-smooth

其实算法并不能提前发现是否所有质数因子 $p$ 的 $p-1$ 都不满足 $k$-smooth。实际上，指定 $k$ 之后并不总能发现质数因子。在上面 $N=221$ 的例子中，如果一开始假设存在质数因子 $p$ 且 $p-1$ 是 $3$-smooth，$B=E_2E_3=2^73^4$，

pairs = [
    (2, 7),
    (3, 4),
]

value = 2 
for p, e in pairs:
    value = bar(value, p**e, 221)
print(f'value={value}, gcd={gcd(value-1, 221)}')

就会出现$\gcd(2^{B-1}, N)$ 等于 $1$ 或 $N$ 的情况。假如运气不好，指定的 $k$ 无法让算法找到正确的因子，那么可以调整、增大 $k$ 的值，利用上一个 $k’$ 对应的 $2^{B_k’}$ 为底数、$E_k$ 为指数计算出剩余 $2^{B_k’E_k}-1\mod(N)$，然后再计算公因子。上面的python代码就是实现了这个逻辑。

但即便是反复增大 $k$ 的值，这个算法也不见得可以在极端情况下高效地算出结果。因为 $p>2$ 是质数，所以 $p-1$ 是一个偶数：

如果 $q$ 的值非常大、甚至接近 $p$ 的长度（比如1000bits以上），算法得提前找到大于这个长度的质数；しかし、质数测试也不是一件容易的事。

• Previous
  HTTPS 反向代理
• Next
  数论学习笔记 - Group, ring and field