快速幂

快速幂

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

让我们先来思考一个问题：7的10次方，怎样算比较快？

方法1：最朴素的想法，77=49，497=343，… 一步一步算，共进行了9次乘法。
这样算无疑太慢了，尤其对计算机的CPU而言，每次运算只乘上一个个位数，无疑太屈才了。这时我们想到，也许可以拆分问题。
方法2：先算7的5次方，即77777，再算它的平方，共进行了5次乘法。
但这并不是最优解，因为对于“7的5次方”，我们仍然可以拆分问题。
方法3：先算77得49，则7的5次方为4949*7，再算它的平方，共进行了4次乘法。
模仿这样的过程，我们得到一个在 [公式] 时间内计算出幂的算法，也就是快速幂。

递归快速幂

刚刚我们用到的，无非是一个二分的思路。我们很自然地可以得到一个递归方程：

计算a的n次方，如果n是偶数（不为0），那么就先计算a的n/2次方，然后平方；如果n是奇数，那么就先计算a的n-1次方，再乘上a；递归出口是a的0次方为1。
递归快速幂的思路非常自然，代码也很简单（直接把递归方程翻译成代码即可）：
//递归快速幂

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

注意，这个temp变量是必要的，因为如果不把a^n/2记录下来，直接写成qpow(a, n /2)*qpow(a, n /2)，那会计算两次a^n/2，整个算法就退化为了 O(n) 。
在实际问题中，题目常常会要求对一个大素数取模，这是因为计算结果可能会非常巨大，但是在这里考察高精度又没有必要。这时我们的快速幂也应当进行取模，此时应当注意，原则是步步取模，如果MOD较大，还应当开long long。
//递归快速幂（对大素数取模）

#define MOD 1000000007
typedef long long ll;
ll qpow(ll a, ll n)
{
    if (n == 0)
        return 1;
    else if (n % 2 == 1)
        return qpow(a, n - 1) * a % MOD;
    else
    {
        ll temp = qpow(a, n / 2) % MOD;
        return temp * temp % MOD;
    }
}

大家知道，递归虽然简洁，但会产生额外的空间开销。我们可以把递归改写为循环，来避免对栈空间的大量占用，也就是非递归快速幂。

非递归快速幂

我们先看代码，再来仔细推敲这个过程：
//非递归快速幂

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^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;
}

注意，较复杂类型的快速幂的时间复杂度不再是简单的 O(logn) ，它与底数的乘法的时间复杂度有关。

例如，矩阵快速幂的一个经典应用是求斐波那契数列：

（洛谷P1962） 斐波那契数列

#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 + ans.a2) % MOD);
    return 0;
}