快速幂算法的基础和扩展

快速幂

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

[a^b = underbrace{a times a times cdots times a}_{btext{次}} ]

这种方法需要进行 (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
2. 从最低位开始检查指数的二进制位
3. 如果当前位为 1，将当前的底数（(a^{2^x})）乘入结果
4. 底数不断平方（不断计算 (a^0,a^1,a^2...)），指数右移一位
5. 重复直到指数的最高位 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++ <cmath>、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}]