基础算法—快速幂详解

幂运算是非常常见的一种运算，求取$a^n$,最容易想到的方法便是通过循环逐个累乘，其复杂度为$O(n)$，这在很多时候是不够快的，所以我们需要一种算法来优化幂运算的过程。

一、快速幂——反复平方法

该怎样去加速幂运算的过程呢？既然我们觉得将幂运算分为n步进行太慢，那我们就要想办法减少步骤，把其中的某一部分合成一步来进行。

比如，如果$n$能被2整除，那我们可以先计算一半，得到$a^{n/2}$的值，再把这个值平方得出结果。这样做虽然有优化，但优化的程度很小，仍是线性的复杂度。

再比如，如果我们能找到$2^k = n$，那我们就能把原来的运算优化成$((a^2)^2)^2...$，只需要$k$次运算就可以完成，效率大大提升。可惜的是，这种条件显然太苛刻了，适用范围很小。不过这给了我们一种思路，虽然我们很难找到$2^k = n$，但我们能够找到$2^{k_1} + 2^{k_2} + 2^{k_3} +......+ 2^{k_m} = n$。这样，我们可以通过递推，在很短的时间内求出各个项的值。

我们都学习过进制与进制的转换，知道一个$b$进制数的值可以表示为各个数位的值与权值之积的总和。比如，2进制数$1001$，它的值可以表示为10进制的$1 imes2^3 + 0 imes2^2 + 0 imes2^1 + 1 imes2^0$，即$9$。这完美地符合了上面的要求。可以通过2进制来把$n$转化成$2^{k_m}$的序列之和，而2进制中第$i$位（从右边开始计数，值为$1$或是$0$）则标记了对应的$2^{i - 1}$是否存在于序列之中。譬如，$13$为二进制的$1101$，他可以表示为$2^3 + 2^2 + 2^0$，其中由于第二位为$0$，$2^1$项被舍去。

如此一来，我们只需要计算$a、a^2、a^4、a^8......a^{2^{k_m}}$的值（这个序列中的项不一定都存在，由$n$的二进制决定）并把它们乘起来即可完成整个幂运算。借助位运算的操作，可以很方便地实现这一算法，其复杂度为$O(log n)$。

typedef long long ll;
ll mod;
ll qpow(ll a, ll n)//计算a^n % mod
{
    ll re = 1;
    while(n)
    {
        if(n & 1)//判断n的最后一位是否为1
            re = (re * a) % mod;
        n >>= 1;//舍去n的最后一位
        a = (a * a) % mod;//将a平方
    }
    return re % mod;
}

取模运算一般情况下是需要的，当然也可以省去。

二、矩阵快速幂

快速幂只是通过二进制拆分$n$来加速幂运算的手段，当然并不只适用于求取数字的幂次，对于矩阵的$n$次方，也可以用同样的手段求取。除了乘法的规则与上面的快速幂不同之外，其他方面并没有太大的差别。

不过这有什么意义呢？利用矩阵的幂次，我们可以快速地完成递推。

比如，在POJ3070 Fibonacci中，就需要我们快速地求取斐波那契数列的第$n$项（取模），对于$n=10^{18}$，一步步推过去显然太慢了，那么我们可以考虑构造矩阵来帮我们完成递推。

首先复习一下矩阵的乘法:

$$left[egin{array}{cccc}{a_{11}} & {a_{12}} & {dots} & {a_{1 n}} \ {a_{21}} & {a_{22}} & {dots} & {a_{2 n}} \ {dots} & {dots} & {dots} & {dots} \ {a_{n 1}} & {a_{n 2}} & {dots} & {a_{n n}}end{array} ight] left[egin{array}{cccc}{b_{11}} & {b_{12}} & {dots} & {b_{1 n}} \ {b_{21}} & {b_{22}} & {dots} & {b_{2 n}} \ {dots} & {dots} & {dots} & {dots} \ {b_{n 1}} & {b_{n 2}} & {dots} & {b_{n n}}end{array} ight]=left[egin{array}{cccc}{c_{11}} & {c_{12}} & {dots} & {c_{1 n}} \ {c_{21}} & {c_{22}} & {dots} & {c_{2 n}} \ {dots} & {dots} & {dots} & {dots} \ {c_{n 1}} & {c_{n 2}} & {dots} & {c_{n n}}end{array} ight]$$

其中$C_{i j}=sum_{k=1}^{n} a_{i k} * b_{k j}$

对于$i ge 3$

有$Fib_i = Fib_{i-1} + Fib_{i-2}$

可以构造出矩阵递推式$$left[egin{array}{ll}{1} & {1} \ {1} & {0}end{array} ight]left[egin{array}{c}{F i b_{i}} \ {F i b_{i-1}}end{array} ight]=left[egin{array}{c}{F i b_{i+1}} \ {F i b_{i}}end{array} ight]$$

那么$$left[egin{array}{ll}{1} & {1} \ {1} & {0}end{array} ight]^{n-2}left[egin{array}{l}{1} \ {1}end{array} ight]=left[egin{array}{c}{F i b_{n}} \ {F i b_{n-1}}end{array} ight]$$

我们就可以利用矩阵快速幂以$O(log n)$求取斐波那契数列的第$n$项了。

const ll mod = 10000;
const int maxv = 2;

struct Matrix {
    ll a[maxv][maxv]; //矩阵

    Matrix operator*(const Matrix &b) const& { 
        //矩阵乘法，复杂度O(maxv^3)，也可看作常数，但maxv较大（大于5）时会使运算时间提高好几个数量级
        Matrix ans;
        for (int i = 0; i < maxv; ++i) {
            for (int j = 0; j < maxv; ++j) {
                ans.a[i][j] = 0;
                for (int k = 0; k < maxv; ++k) {
                    ans.a[i][j] += a[i][k] * b.a[k][j] % mod;
                    ans.a[i][j] %= mod;
                }
            }
        }
        return ans;
    }
    static Matrix qpow(Matrix x, ll n) {//矩阵快速幂，将乘法复杂度看作常数则复杂度为O(log n)
        Matrix ans;
        for (int i = 0; i < maxv; ++i) {
            for (int j = 0; j < maxv; ++j) {
                if (i == j)
                    ans.a[i][j] = 1;
                else
                    ans.a[i][j] = 0;
            }
        }//初始化为单位矩阵，参考普通数字的快速幂这里初始化为1
        while (n) {//其余细节基本相同
            if (n & 1)
                ans = ans * x;
            x = x * x;
            n >>= 1;
        }
        return ans;
    }
    Matrix(ll temp[maxv][maxv]) {//构造方法
        for (int i = 0; i < maxv; ++i) {
            for (int j = 0; j < maxv; ++j) {
                a[i][j] = temp[i][j];
            }
        }
    }
    Matrix() { }
};

ll fib(ll n) {//求取斐波那契数列第n项（本题取模）
    if (n == 0)
        return 0;
    if (n <= 2)
        return 1;
    
    ll temp[maxv][maxv] = {
        1, 1,
        1, 0
    };
    Matrix m(temp);
    m = Matrix::qpow(m, n - 2);
    return (m.a[0][1] + m.a[0][0]) % mod;
}

 当然了，构造矩阵的方法是多种多样的，例如上面的题目中就给出了另一种构造出斐波那契数列的方法。

矩阵快速幂能解决的递推也远不止斐波那契数列，下面根据我个人的习惯列举几个比较常见、比较“套路”的情况：

（1）对于$a_n = xcdot a_{n-1} + ycdot a_{n-2} + c^n$，可以得到

$$
left[egin{array}{ccc}{x} & {y} & {c} \ {1} & {0} & {0} \ {0} & {0} & {c}end{array} ight]left[egin{array}{c}{a_{i}} \ {a_{i-1}} \ {c^{i}}end{array} ight]=left[egin{array}{c}{a_{i+1}} \ {a_{i}} \ {c^{i+1}}end{array} ight]
$$

（2）对于$a_n = a_{n-1} + a_{n-2} + n^2$，可以得到

$$
left[egin{array}{ccccc}{1} & {1} & {1} & {2} & {1} \ {1} & {0} & {0} & {0} & {0} \ {0} & {0} & {1} & {2} & {1} \ {0} & {0} & {0} & {1} & {1} \ {0} & {0} & {0} & {0} & {1}end{array} ight]left[egin{array}{c}{a_{i}} \ {a_{i-1}} \ {i^{2}} \ {i} \ {1}end{array} ight]=left[egin{array}{c}{a_{i+1}} \ {a_{i}} \ {(i+1)^{2}} \ {i+1} \ {1}end{array} ight]
$$

（3）对于$a_n=xcdot a_{n-1} + ycdot a_{n-2}$，求$a_n$的前缀和$s_n$

$$
left[egin{array}{lll}{1} & {x} & {y} \ {0} & {x} & {y} \ {0} & {1} & {0}end{array} ight]left[egin{array}{c}{s_{i}} \ {a_{i}} \ {a_{i-1}}end{array} ight]=left[egin{array}{c}{s_{i+1}} \ {a_{i+1}} \ {a_{i}}end{array} ight]
$$