上周做了一道简单 C 语言练习题，做完了才反应过来是快速幂...😶

记得上一次看到快速幂相关的内容已经是好几年前的事情了？当时看的一知半解，也没太在意，现在又做到了，竟然莫名觉得很清晰，是我变强了吗？哈哈~
不管怎么说，这次得把这个简单实用的算法搞清楚了！

Intro

快速幂算法，也叫二进制求幂（Binary Exponentiation），算法的名字就是算法的用途。算法的思路也比较简单，就是通过减少乘法操作的次数来达成快速求幂。在常规求幂的过程中，$a^n=\underbrace{a \times a \cdots \times a}_n$，而快速幂则是将每次得到的幂作为下一次的乘数和被乘数进行相乘，这样无形之中减少了乘法的次数。

Sample

以计算 7 的 10 次方为例，说明一下常规思路和快速幂思路。
**常规思路：**那就是 $7 \times 7 = 49$，$49 \times 7 = 343$，... 一步一步算，一共要进行 9 次乘法。

**快速幂：**先算 $7 \times 7 = 49$，然后算 $49 \times 49 = 343$，再算 $7^4$ 的平方，最后在乘以一个 $7$ 就可以得出结果了，一共是 4 次乘法。

从这里也可以看出，常规思路的时间复杂度是 $O(n)$，快速幂的时间复杂度是 $O(log n)$。

Implementation

Recursion

快速幂也是基于二分的，可以考虑从递归的角度来实现，这里直接借用上周题目中所给的递归函数：
$$
fun(a, b) =
\begin{cases}
{1} & \text{if b = 0;} \
{fun(a \times a, b / 2)} & \text{else if b mod 2 = 0;} \
{fun(a \times a, b / 2) \times a} & \text{else.} \
\end{cases}$$

这是当时的实现：

ver1

#define MOD 1000007
int fun(int a, int b) {
  if(b == 0) return 1;
  else if(b % 2 == 0) return (fun(a * a, b / 2) % MOD);
  else return ((fun(a * a, b / 2) % MOD) * (a % MOD));
}

因为在计算过程中，结果可能会非常大，为了保证最后结果的正确性需要进行取余操作，另外还建议改用long long来避免溢出，所以可以得到下面的版本：

ver2

#define MOD 1000007
long long fun(long long a, long long b) {
  if(b == 0) return 1;
  else if(b % 2 == 0) return (fun(a * a, b / 2) % MOD);
  else return ((fun(a * a, b / 2) % MOD) * (a % MOD));
}

Non-Recursion

想要将上述递归算法改为非递归，就需要从二进制的角度来思考问题。
比如，要计算 $7^{10}$，就有 $7^{(1010)} = 7^{(1000)} \times 7^{(0010)}$，所以非递归算法的关键在于利用位运算计算乘数，最后可以得到下面的版本：

ver3

long long fun(long long a, long long n) {
    long long ans = 1;
    while(n) {
        if(n & 1)
            ans = ((ans % MOD) * (a % MOD)) % MOD;
        a *= a;
        n >>= 1;
    }
    return ans;
}

这里，借用他人文章中的一张图来说明整个计算过程。
$7^{10}$

回头来看，递归快速幂的思路其实也是二进制，也即对指数不断进行取余和除 2 操作，这个过程其实就是除二取余法。

Extension

扩展部分主要引用他人文章中的内容，包括一些自己的理解和应用。
首先是这类思路的模板：

// 泛型的非递归快速幂
template <typename T>
T qpow(T a, ll n) {
    T ans = 1; // 赋值为乘法单位元，可能要根据构造函数修改
    while (n) {
        if (n & 1)
            ans = ans * a; // 这里就最好别用自乘了，不然重载完*还要重载*=，有点麻烦。
        n >>= 1;
        a = a * a;
    }
    return ans;
}

还有一个矩阵快速幂的例子：洛谷 P1962 斐波那契数列以及对应的代码：

P1962

#include <cstdio>
#define MOD 1000000007
typedef long long ll;

struct matrix
{
    ll a1, a2, b1, b2;
    // constructor with initialize list
    matrix(ll a1, ll a2, ll b1, ll b2) : a1(a1), a2(a2), b1(b1), b2(b2) {}
    // operator * overloading
    matrix operator*(const matrix &y)
    {
        matrix ans((a1 * y.a1 + a2 * y.b1) % MOD,
                   (a1 * y.a2 + a2 * y.b2) % MOD,
                   (b1 * y.a1 + b2 * y.b1) % MOD,
                   (b1 * y.a2 + b2 * y.b2) % MOD);
        return ans;
    }
};

matrix qpow(matrix a, ll n)
{
    matrix ans(1, 0, 0, 1); //单位矩阵
    while (n)
    {
        if (n & 1)
            ans = ans * a;
        a = a * a;
        n >>= 1;
    }
    return ans;
}

int main()
{
    ll x;
    matrix M(0, 1, 1, 1);
    scanf("%lld", &x);
    matrix ans = qpow(M, x - 1);
    printf("%lld\n", (ans.a1 + ans.a2) % MOD);
    return 0;
}

这个小题用到了三个重要的知识点：

矩阵快速幂
构造函数初始化列表
运算符重载

真要用纯用 C 语言的东西，估计代码量是两倍了...🤔

参考内容：
1. 算法学习笔记(4)：快速幂

bachzart's Blog

浅析快速幂算法