定义

快速幂（Exponentiation by squaring，平方求幂）是一种简单而有效的小算法，它可以以 $O(\log n)$ 的时间复杂度计算乘方。快速幂不仅本身非常常见，而且后续很多算法也都会用到快速幂。

背景问题

例如计算一个 $7^{10}$ ，有以下几种方法。

最朴素的想法，7*7=49，49*7=343，… 一步一步算，共进行了9次乘法。
- 这样算无疑太慢了，尤其对计算机的CPU而言，每次运算只乘上一个个位数，无疑太屈才了。这时我们想到，也许可以拆分问题。
先算7的5次方，即7*7*7*7*7，再算它的平方，共进行了5次乘法。
- 但这并不是最优解，因为对于“7的5次方”，我们仍然可以拆分问题。
先算7*7得49，则7的5次方为49*49*7，再算它的平方，共进行了4次乘法。

模仿这样的过程，我们得到一个在 $O(\log n)$ 时间内计算出幂的算法，也就是快速幂。

具体实现

递归快速幂

用上面的背景问题可以得到下面这个递归公式：

$a^n = \begin{cases} a^{n-1}\times a, & n\ is\ odd \\ a^{\frac{n}{2}}\times a^{\frac{n}{2}}, & n\ is\ even \ but \ not \ 0 \\ 1,&n=0 \end{cases}$

依据这个递归公式，可以写出递归版本的快速幂。

int qpow(int a,int n){
    if(n == 0){
        return 1;
    }else if(n%2 == 1){
        return qpow(a,n-1)*a; //n是奇数
    }else{
        int temp=qpow(a,n/2); //n是偶数
        return temp*temp;
    }
}

代码中有一个temp变量，暂存 $a^{\frac{n}{2}}$ 的值，这样就能避免计算两次。

非递归快速幂

虽然递归版本很容易写，但是会产生额外的空间开销。因此可以把递归版本的改写为循环版本的。
将 $7^{10}$ 转化为 $y^2 \times 7^8$ 。

int qpow(int a, int n){
    int ans = 1;
    while(n){
        if(n&1){       //如果n的当前末位为1
            ans *= a;  //ans乘上当前的a
        }
        a *= a;        //a自乘
        n >>= 1;       //n往右移一位
    }
    return ans;
}

如果加上取模过程：
有 $(a*b)\%c = ((a\%c)*(b\%c))\%c$

int qpow(int a, int n){
    int ans = 1;
    while(n){
        if(n&1){       //如果n的当前末位为1
            ans *= a;  //ans乘上当前的a
            ans %= mod;
        }
        a *= a%mod;        //a自乘
        a %= mod;
        n >>= 1;       //n往右移一位
    }
    return ans;
}

感觉有点类似计算 $2*10$ 时转换为 $2*8+2*2==(2<<3)+(2<<1)$ 。 $(10D==1010B)$

矩阵快速幂

其实，在算 $a^n$ 时，只要a的数据类型支持乘法且满足结合律，快速幂的算法都是有效的。矩阵、高精度整数，都可以照搬这个思路。

//泛型的非递归快速幂
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 斐波那契数列
设矩阵

$A=\begin{bmatrix} 0&1\\ 1&1 \end{bmatrix}$

我们有：

$A\begin{bmatrix} F_n\\F_{n+1} \end{bmatrix}=\begin{bmatrix} F_{n+1}\\F_n+F_{n+1} \end{bmatrix}=\begin{bmatrix} F_{n+1}\\F_{n+2} \end{bmatrix}$

于是 :

$\begin{bmatrix} F_n\\F_{n+1} \end{bmatrix}= A\begin{bmatrix} F_{n-1}\\F_{n} \end{bmatrix}=A^2\begin{bmatrix} F_{n-2}\\F_{n-1} \end{bmatrix}=\cdots=A^{n-1}\begin{bmatrix} F_1\\F_{2} \end{bmatrix}=A^{n-1}\begin{bmatrix} 1\\1 \end{bmatrix}$

这样，我们把原来较为复杂的问题转化成了求某个矩阵的幂的问题，这就可以应用快速幂求解了。

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

struct matrix
{
    ll a1, a2, b1, b2;
    matrix(ll a1, ll a2, ll b1, ll b2) : a1(a1), a2(a2), b1(b1), b2(b2) {}
    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 % MOD);
    return 0;
}