幂运算(Exponentiation)又称指数运算，写作\(c= m^e \)，即\(c\)等于\(e\)个\(m\)相乘。当\(m\)或者\(e\)的数值不大的时候，可以通过重复的乘法运算来获得\(c\)。当\(m\)和\(e\)都是一个大数的时候，比如说512位长。尽管有了计算机，显然对\(m\)重复\(e\)次乘法运算的效率不高。那么该怎样快速进行大数幂运算呢？
平方乘
为了找到快速计算\(c= m^e\)的方法，首先检查一下\(m\)和\(e\)的特点。在计算机存储或者计算数字的时候，它们都是以二进制的形式存在。例如\(m^{5}\)在计算机看来应该是 $$ m^{5} = m^{101}。$$ 先简单描述一个新的算法来求 \(c = m^5 \)：令 \(c = 1\)，然后从101的最高位开始，先求其平方\(c = c^2 \)，如果该位的值为1，则令\(c = c \cdot m\)，否则进入下一位重复这般操作。现在\(c = m^5 \)的运算过程如下 $$ c = m^{101} = ((1^2 \cdot m)^2)^2 \cdot m。$$ 通过观察二进制指数，可以发现并不是\(e\)里面每一个比特都参与乘法运算，只有为1的比特才会引起一次乘法运算。对比原来的5次乘法，新的算法只需要2次平方(1的平方不算)和2次乘法。 当\(e=25\)时，优化效果更加明显 $$ m^{25} = m^{11001} = ((((1 \cdot m)^2 \cdot m)^2)^2)^2 \cdot m，$$ 原本需要25次乘法运算，现在通过4次平方和3次乘法即可得到结果。这个算法被称作Square&Multiply Alogrithm，译作平方乘算法。在RSA加密系统中，密文\(c = m^e \mod N\)。通常，此时\(m, e 和N\)都是大数，高效运算对于RSA的应用有着重要意义。让我们以RSA为例，用伪代码更精确地描述一下原算法。
Input: Base \(m\), modulus \(N\), exponent \( e = \sum_{i=0}^l e_i2^i, \ e_l = 1\)Output: \( c = m^e \mod N\)
- \( c = m\)
- FOR \(i = l-1\) DOWN TO \(0\) DO:
- 2.1 \( c = c^2 \mod N\)
- 2.2 IF \(e_i = 1\) THEN:
- 2.2.1 \(c = c \cdot m \mod N\)
- RETURN \(c\)
现在来研究一下平方乘算法需要的计算步骤。先看下表
| 操作 | 次数 |
|---|---|
| SQR | \(l\) |
| MUL | HW(e)-1 |
注意一下，里面的\(l\)不是指数的长度，而是其长度减1。HW代表的是Hamming Weight，HW(\(e\))的值为在\(e\)里面存在的1的数量。在算法的第一步只是赋值操作，并没有平方也没有乘法操作，所以最终SQR的次数为\(l\)，MUL次数为HW(\(e\))-1。
还以\( c = m^{25} = m^{11001}\)为例，仔细看一下算法的运行过程。
| 迭代 | 指数 | 结果 | 操作 |
|---|---|---|---|
| 0 | 1 | \( c = m \) | |
| 1a | 1 | \(c = c^2 \) | SQR |
| 1b | 11 | \( c = c \cdot m \) | MUL |
| 2a | 11 | \( c = c^2 \) | SQR |
| 2b | 110 | - | - |
| 3a | 1100 | \( c = c^2 \) | SQR |
| 3b | 1100 | - | - |
| 4a | 11001 | \( c = c^2 \) | SQR |
| 4b | 11001 | \( c = c \cdot m \) | MUL |
在表中步骤a代表是一次平方操作，步骤b代表当指数的该位值为1时候进行的一次乘法操作。
k比特幂运算
在上述的算法里面，一次迭代只向右推进一位比特，速度还不够快。现在假设\(k\)个比特是一个整体，也就是将原来的二进制指数转换成\(b=2^k\)进制数，即\(e=(E_tE_{t-1} \cdots E_1E_0)_b\)。这样就可以迭代一次推进\(k\)个比特。然后不是乘以\(m\)，而是\(m^{E_i}\), \(E_i \in \{0,1,\cdots,2^k-1\}\)。通过事先计算出所有的\(m^{E_i}\)，并存在表格LUT里面，当进行乘法操作的时候，直接根据当时的指数\(E_i\)，查表得出\(m^{E_i}\)然后相乘，可以进一步加快速度。k比特算法的伪代码如下
Input: Base \(m\), modulus \(N\), Exponent \(e = (E_tE_{t-1} \cdots E_1E_0)_b\), where \(b = 2^k\) for some \(k > 1\), \(E_i \in \{0, 1, \cdots, 2^k -1\}\)Output: \(c = m^e \mod N\)
- Precomputation
- 1.1 \(m_0 = 1\)
- 1.2 FOR \(i\) FROM \(1\) TO \((2^k -1)\) DO:
- 1.2.1 \(m_i = m_{i-1} \cdot m \mod N\)
- \(c = m^{E_t}\)
- FOR \(i\) FROM \(t-1\) DOWN TO \(0\) DO:
- 3.1 \(c = c^{2^k} \mod N\)
- 3.2 \(c = c \cdot m^{E_i} \mod N\)
- RETURN \(c\)
同样，也用表格来展示迭代的过程。
| 迭代 | 结果 | 操作 |
|---|---|---|
| 0 | \(m^{E_t} \) | TLU |
| 1a | \( (m^{E_t})^{2^k} = m^{(E_t0)_b} \) | k·SQR |
| 1b | \(m^{(E_t0)_b} \cdot m^{E_{t-1}} = m^{(E_tE_{t-1})_b} \) | TLU&MUL |
| 2a | \((m^{(E_tE_{t-1})_b})^{2^k} = m^{(E_tE_{t-1}0)_b} \) | k·SQR |
| 2b | \(m^{(E_tE_{t-1}0)_b} \cdot m^{E_{t-2}} = m^{(E_tE_{t-1}E_{t-2})_b}\) | TLU&MUL |
| \(\vdots\) | \(\vdots\) | \(\vdots\) |
| ta | \((m^{(E_t \cdots E_0)_b})^{2^k} = m^{(E_t \cdots E_1 0)_b} \) | k·SQR |
| tb | \(m^{(E_t \cdots E_1 0)_b} \cdot m^{E_0} = m^{(E_t \cdots E_1 E_0)_b}\) | TLU&MUL |
更小的表格
在\(b\)进制的情况下面，再仔细看一下指数里面的\(E_i\)。如前所述，任何数字在计算机里面都是以二进制存储。可以将其表示为 $$ E_i = (\overset{k-h_i}{\overline{e_{k-1} e_{k-2} \cdots e_{h_i}}} \ \underset{h_i}{\underline{0 \cdots 0}})_2$$ 其中\(0 \leq h_i \leq k-1\)，它代表在\(E_i\)里面值为1的最低位，\(e_{h_i} = 1\)。定义 \( u_i = e_{k-1} e_{k-2} \cdots e_{h_i} \)，于是把\(E_i\)分成了两部分，全为0的低位，长度是\(h_i\)，值为\(u_i\)的高位，长度是\(k-h_i\)。当\(E_i = 0\)时，则令\(u_i = h_i = 0 \)。于是， $$ E_i = u_i \cdot 2^{h_i}。$$因为\(u_i < 2^k\)，所以\(m^{u_i}\)存在于表LUT中。同时\(u_i\)总是奇数，因此并不需要查找表中的偶数项。在预先计算的步骤，只需要存储奇数项索引和值即可，因此关于表格的运算量减少了一半。
也许有同学会在这里感到奇怪，难道就不需要偶数项指数了吗？回想一下最开始1比特的情况：0位不会引起乘法运算，只有1位才会，但是每一位比特都会引起一次平方运算。现在再来回过头，看一下\(k\)比特情况：我们把指数\(E_i\)分成了奇数的高位，与值为0的低位。高位比特的值为\(u_i\)，它会引起一次乘以\(m^{u_i}\)的操作，同时也有\(k-h_i\)次平方操作。低位比特只会引起\(h_i\)次平方操作。
现在修改算法，它的伪代码如下
Input: Base \(m\), modulus \(N\), Exponent \(e = (E_tE_{t-1} \cdots E_1E_0)_b\), where \(b = 2^k\) for some \(k > 1, E_i = u_i \cdot 2^{h_i}\)Output: \(c = m^e \mod N\)
- Precomputation
- 1.1 \(m_0 = 1, m_1 = m, m_2 = m^2 \mod N\)
- 1.2 FOR \(i\) FROM \(1\) to \((2^{k-1} - 1)\) DO:
- 1.2.1 \(m_{2i+1} = m_{2i-1} \cdot m_2\)
- \(c = 1\)
- FOR \(i\) FROM \(t\) DOWN TO \(0\) DO:
- 3.1 \(c = c^{2^{k-h_i}} \mod N\)
- 3.2 \(c = c \cdot m_{u_i} \mod N\)
- 3.3 \(c = c^{2^{h_i}} \mod N\)
表格可以更加直观地展示算法运行过程
| 迭代 | 结果 | 操作 |
|---|---|---|
| 3.1 | \((c^{(E_t \cdots E_{i+1})_b})^{2^{k - h_i}} = c ^{(E_t \cdots E_{i+1})_b(\overset{k-h_i}{\overline{0 \cdots 0}})_2} \) | \((k − h_i ) \cdot \)SQR |
| 3.2 | \((c ^{(E_t \cdots E_{i+1})_b(0 \cdots 0)_2}) \cdot m^{u_i} = c^{(E_t \cdots E_{i+1})_b(u_i)_2}\) | TLU&MUL |
| 3.3 | \((c^{(E_t \cdots E_{i+1})_b(u_i)_2})^{2^{h_i}} =c^{(E_t \cdots E_{i+1})_b(\underset{ = E_i}{\underline{u_i \overset{h_i}{\overline{0 \cdots 0}}} )_2}}\) | \( h_i \cdot \)SQR |
滑动窗口
所谓滑动窗口指的是，每次向右移动的距离并不固定，窗口的大小根据窗口内的比特信息而变化。滑动窗口技术是最高效的大数幂运算。现在来看一下究竟它是怎样工作的。
首先给窗口定义一个最大长度用来限制其范围: \(w \leq k\)。然后规定窗口的滑动规则：窗口向右滑动过程中，遇见0位直接跳过，直到遇到1，这是窗口开始的位置，然后继续向右寻找窗口的最后一位。在不超过最大长度的情况下，窗口最后一位比特也是1。在算法运行过程中，窗首和窗尾都是在不停地滑动，直到其遇到一个首尾均为1的最长序列。如下例所示，\(k = 4\) $$ \cdots \underset{w_i}{\boxed{1101}} \ \underset{w_{i+1}}{\boxed{101}} \ 000 \ \underset{w_{i+2}}{\boxed{1001}} \cdots$$
有了滑动窗口的辅助，算法对指数的操作粒度变得更细，但又不至于逐位比特计算，平衡了时间和空间的效率。
同样地，最后用伪代码描述一下大数幂运算最高效的算法
Input: Base \(m\), modulus \(N\), Exponent \(e = (e_le_{l-1} \cdots e_1e_0)_2\), with \(e_l = 1\) and \(k \geq 1\)Output: \(c = m^e \mod N\)
- Precomputation
- 1.1 \(m_0 = 1, m_1 = m, m_2 = m^2 \mod N\)
- 1.2 FOR \(i\) FROM \(1\) to \((2^{k-1} - 1)\) DO:
- 1.2.1 \(m_{2i+1} = m_{2i-1} \cdot m_2 \mod N \)
- \(c = 1, i = l\)
- WHILE \(i \geq 0\) DO:
- 3.1 IF \(e_i = 0\) THEN DO:
- 3.1.1 \(c = c^2 \mod N \)
- 3.1.2 \(i = i-1 \)
- 3.2 ELSE find the longest bitstring \(e_ie_{i-1} \cdots e_j\) such that \(i-j+1 \leq k\) and \(e_j = 1\) and DO:
- 3.2.1 \(c = c^{2^{i-j+1}} \cdot m_{(e_ie_{i-1} \cdots e_j)_2}\mod N\)
- 3.2.2 \(i = j -1\)
- 3.1 IF \(e_i = 0\) THEN DO:
- RETURN \(c\)
结束语
笔者在做密码学实践课作业的时候，老师要求实现一个大数幂运算算法。在C语言有一个厍叫做GMPLib，即gmp.h头文件。里面的函数可以操作任意精度的数值，而且速度非常得快。其中一个名叫mpz_powm()的函数，就可以实现快速幂运算。笔者的工作其实就是自己实现了这样一个函数的算法，代码在这里可以参阅。这篇介绍详细介绍大数幂运算的小文，撰写难度令我几度停手，写 \(\LaTeX \) 确实废指 🤘
特别希望此文能够帮助到一些后来者。