快速幂算法的基础和扩展

忙贬 · 7 天前

快速幂

快速幂（Fast Exponentiation）算法解决这样一个问题：求解自然数的指数运算。计算 $a^b$ 时，按照指数定义的朴素的方法是通过连续相乘：

\[a^b = \underbrace{a \times a \times \cdots \times a}_{b\text{次}}\]
这种方法需要进行 $b-1$ 次乘法，当 $b$ 很大时（如 $10^9$），时间复杂度 $O(b)$ 是完全不可接受的。
快速幂通过巧妙的二进制分解技术，将幂运算的时间复杂度从 $O(b)$ 优化到 $O(\log b)$。
考虑计算 $a^{13}$，将指数 13 用二进制表示：

\[13 = 1101_2 = 2^3 + 2^2 + 0 + 2^0 = 8 + 4 + 0 + 1\]
因此：

\[a^{13} = a^{8 + 4 + 0 + 1} = a^8 \times a^4 \times a^0 \times a^1\]
而 $ a^8 = (a^4) ^2 = ((a^2) ^2) ^2$ ，分解后的幂次很容易计算
算法流程：

初始化结果 1
从最低位开始检查指数的二进制位
如果当前位为 1，将当前的底数（$a^{2^x}$）乘入结果
底数不断平方（不断计算 $a^0,a^1,a^2...$），指数右移一位
重复直到指数的最高位 1 也被遍历

快速幂算法也可以从递归的角度来理解，这种理解方式更加直观。

\[a^b = \begin{cases}1 & \text{if } b = 0 \\(a^{b/2})^2 & \text{if } b \text{ is even} \\a \times (a^{(b-1)/2})^2 & \text{if } b \text{ is odd}\end{cases}\]

long long quick_pow(long long base, long long exp) {
long long res = 1;
while (exp) {
if (exp & 1) {
res *= base;
}
base *= base;
exp >>= 1;
}
return res;
}

复制代码

带模数版本

更多的时候，我们要求解的是 $a^b \bmod m$。这也可以用快速幂思想解决。快速幂模数版本的正确性基于模运算的分配律：

\[(a \times b)\bmod m = [(a \bmod m) \times (b \bmod m)] \bmod m\]
因此，我们可以直接对 $ a $ 取模，同时在算法每一步中，我们都对中间结果取模，这保证了最终结果的正确性，同时防止数值溢出。
不过我们不能直接对指数取模。指数 $b$ 必须保持原值，因为：

\[a^b \bmod m \neq a^{b \bmod m} \bmod m\]
当然，既然复杂度是对数的，所以 $b$ 大一些一般也无所谓。

long long quick_pow(long long base, long long exp, long long mod) {
long long res = 1;
base %= mod; // 先取模，防止初始base过大
while (exp) {
if (exp & 1) {
res = (res * base) % mod;
}
base = (base * base) % mod;
exp >>= 1;
}
return res;
}

复制代码

不过，在某些特定情况下，我们可以使用欧拉定理来化简指数：
欧拉定理：如果 $a$ 和 $m$ 互质（即 $\gcd(a, m) = 1$），那么：

\[a^{\phi(m)} \equiv 1 \pmod{m}\]
其中 $\phi(m)$ 是欧拉函数，表示小于 $m$ 且与 $m$ 互质的正整数的个数。
所以当 $a$ 和 $m$ 互质时，我们可以将指数对 $\phi(m)$ 取模：

\[a^b \mod m = a^{b \mod \phi(m)} \mod m\]
这在某些数学和密码学应用中很有用，但不是快速幂算法的必要部分。代码略。
快速幂方法的时间复杂度是 $O(\log b)$，循环次数等于指数的二进制位数，效率极高。
演示：计算 $3^{13} \mod 100$

指数 13 = 1101(二进制)
初始化: res = 1, base = 3
第1轮 (最低位为1): res = 1×3 = 3, base = 3² = 9
第2轮 (位为0): res = 3, base = 9² = 81
第3轮 (位为1): res = 3×81 = 243 ≡ 43, base = 81² = 6561 ≡ 61
第4轮 (位为1): res = 43×61 = 2623 ≡ 23
结果: 3¹³ ≡ 23 (mod 100)

复制代码

快速乘（防治溢出）

同样的思想也可以应用到乘法本身中。两个 32 位整数相乘，范围将达到 64 位；两个 64 位整数相乘，范围将达到 128 位。同样大小的数无法装入正确的结果。
快速乘（又称"龟速乘"）模仿快速幂的思想，将乘法运算转换为加法运算。核心思路是将 $a \times b$ 看作是 $b$ 个 $a$ 相加，然后利用二进制分解来优化，这样就可以在中间结果下取模。
对于 $a \times b$，将 $b$ 二进制分解：

\[b = \sum_{i=0}^{k} b_i \cdot 2^i \quad \text{其中 } b_i \in \{0,1\}\]
那么：

\[a \times b = a \times \sum_{i=0}^{k} b_i \cdot 2^i = \sum_{i=0}^{k} b_i \cdot (a \cdot 2^i)\]

typedef long long ll;
// 快速乘：返回 (a * b) % mod，防止中间过程溢出
ll quick_mul(ll a, ll b, ll mod) {
ll res = 0;
a %= mod;
while (b) {
if (b & 1) {
res = (res + a) % mod;
}
a = (a + a) % mod; // a = 2a，相当于左移一位
b >>= 1;
}
return res;
}
// 使用快速乘的快速幂
ll quick_pow_safe(ll base, ll exp, ll mod) {
ll res = 1;
base %= mod;
while (exp) {
if (exp & 1) {
res = quick_mul(res, base, mod); // 关键替换！
}
base = quick_mul(base, base, mod); // 关键替换！
exp >>= 1;
}
return res;
}

复制代码

方法时间复杂度空间复杂度防溢出能力直接乘法$O(1)$$O(1)$无快速乘$O(\log n)$$O(1)$有快速乘通过 $O(\log n)$ 次加法替代 $O(1)$ 次乘法，实际上更慢了，所以也叫做“龟速乘”，这属于用时间换取了数值安全性。
浮点幂

如果是底数浮点，指数自然数，那么直接应用快速幂没有任何问题。但若指数是浮点数，这个问题会麻烦的多：浮点数指数无法直接进行二进制位操作，且误差会随着运算的拆分不断累积。
相比之下浮点幂的主要思想是利用自然对数变换法来计算浮点幂：

\[a^b = e^{b \cdot \ln(a)}\]
其中，自然对数和指数是常见且重要的函数，有快速且精确的办法来实现。
常见的库（如C++ 、Intel MKL、GNU Scientific Library）采用此类核心思路。

// 伪代码示意
if (a == 0.0) {
if (b > 0) return 0.0;
if (b == 0) return 1.0; // 或 NaN，依标准而定
return INFINITY; // 或报错
}
if (a == 1.0) return 1.0;
if (b == 0.0) return 1.0;
if (b == 1.0) return a;
result = exp(b * log(a)); # 对于一般情况

复制代码

矩阵快速幂

已知矩阵 $A$，由于矩阵乘法满足结合律，指数为自然数时，仍可以利用快速幂思想求解 $A^n$。这最其深刻、最实用的扩展之一。它将快速幂的核心理念从标量运算成功迁移到了线性代数领域。

\[A^n = \begin{cases} I & \text{if } n = 0 \\(A^{n/2})^2 & \text{if } n \text{ is even} \\A \times (A^{(n-1)/2})^2 & \text{if } n \text{ is odd}\end{cases}\]

typedef vector<vector<long long>> Matrix;
Matrix matrixMultiply(const Matrix& A, const Matrix& B, long long mod) {
int n = A.size();
Matrix C(n, vector<long long>(n, 0));
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
for (int k = 0; k < n; k++) {
C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % mod;
}
}
}
return C;
}
Matrix matrixPow(Matrix base, long long exp, long long mod) {
int n = base.size();
// 初始化单位矩阵
Matrix res(n, vector<long long>(n, 0));
for (int i = 0; i < n; i++) {
res[i][i] = 1;
}
while (exp > 0) {
if (exp & 1) {
res = matrixMultiply(res, base, mod);
}
base = matrixMultiply(base, base, mod);
exp >>= 1;
}
return res;
}

复制代码

经典案例：斐波那契数列的矩阵解法

斐波那契数列的递推关系：

\[F(n) = F(n-1) + F(n-2)\]
可以表示为矩阵形式：

\[\begin{bmatrix} F(n) \\ F(n-1) \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \times \begin{bmatrix} F(n-1) \\ F(n-2) \end{bmatrix}\]
递推得到：

\[\begin{bmatrix} F(n) \\ F(n-1) \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-1} \times \begin{bmatrix} F(1) \\ F(0) \end{bmatrix}\]
由于我们可以快速计算矩阵的幂，我们就绕过了斐波那契数列的定义，使用对数次矩阵乘法的时间直接计算出了某一项。

long long fibonacci_matrix(long long n, long long mod) {
if (n == 0) return 0;
if (n == 1) return 1;
Matrix base = {{1, 1}, {1, 0}};
Matrix result = matrixPow(base, n - 1, mod);
return result[0][0]; // F(n)
}

复制代码

更一般的，对于 k 阶线性递推：

\[a_n = c_1a_{n-1} + c_2a_{n-2} + \cdots + c_ka_{n-k}\]
构造转移矩阵：

\[M = \begin{bmatrix}c_1 & c_2 & \cdots & c_{k-1} & c_k \\1 & 0 & \cdots & 0 & 0 \\0 & 1 & \cdots & 0 & 0 \\\vdots & \vdots & \ddots & \vdots & \vdots \\0 & 0 & \cdots & 1 & 0\end{bmatrix}\]
则：

\[\begin{bmatrix} a_n \\ a_{n-1} \\ \vdots \\ a_{n-k+1} \end{bmatrix} = M^{n-k+1} \times \begin{bmatrix} a_{k-1} \\ a_{k-2} \\ \vdots \\ a_0 \end{bmatrix}\]

来源：程序园用户自行投稿发布，如果侵权，请联系站长删除
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！

账号		自动登录	找回密码
密码			立即注册

快速幂算法的基础和扩展

相关帖子

浏览过的版块

签约作者

快速幂算法的基础和扩展

相关帖子

相关推荐

浏览过的版块

签约作者