快速幂运算

最近看到的一些题里有利用快速矩阵幂来求快速解递归函数（比如Fibonacci数列），在很久之前看和RSA算法里也涉及到了快速求幂的算法，所以今天就想一块来介绍一下。
在介绍开始之前，先对本文的算法做一个声明，下文所提到的幂运算的指数n，满足条件n >= 0，这一点非常重要。
首先幂运算，大家其实都很好理解，其实就是：
无非就是把一个数自乘n次，所以我们通常计算数字的幂的算法就是：
（代码1）

int Pow(int A, int n)
{
    if(n == 0) return 1;
    int rslt(1);
    for(int i(0); i < n; ++i)
    {
        rslt *= A;
    }
    return rslt;
}

这个算法很简单，其复杂度是O(n)的，一般来说，这样的复杂度并不会使我们困惑，但是一般应用幂运算的地方，指数都会非常非常的大，比如1 000 000 000这个级别的，这时候我们会遇到两个问题，第一个就是我们不能再用int来存储整数，必须用高精度整数类型来进行存储，另一个就是在指数是如此变态的数量级之下，我们的计算量会骤然上升，结果也会异常惊人的大。
快速幂算法就是为了解决这个问题而发明的，其实它本身很简单，就是利用二分法，在具体介绍之前，我们需要提一个小小的前提：
那就是能使用快速幂算法的数据类型必须满足结合律，这一点非常重要，用数学描述就是(A * B) * C = A * (B * C)。
在这个前提下，所谓的快速幂运算就是将一个长串分解成一些重复的字串，通过这种方式来避免重复的运算，以提高效率，举个例子，比如：

对于一般的解法：
A^8 = A * A * A * A * A * A * A * A
总共需要7次乘法运算；

将其平均分解：
A^8 = (A * A * A * A) * (A * A * A * A) = (A * A * A * A) ^ 2
这样我们就只需要4次乘法运算了；

我们还可以将其再分解：
A^6 = [(A * A) * (A * A)] ^ 2 = [(A * A) ^ 2] ^ 2
这样就将乘法运算的次数减少为了3次。

通过这样的分解求解，到最后原本需要7次的运算减少为了3次，事实上，这种二分解法可以将原本n次的运算减少为logn / log2，这样的效果是惊人的，在1 000 000 000这样数量级的指数运算下，该方法可以将运算次数减少到30次（OAO!!!），不可思议吧？这就是对数的威力。在这个思想的指导下，我们可以将算法的复杂度修改为不可思议的O(logn)。
进而我们可以得出整数的快速幂算法：
（代码2）

int qPow(int A, int n)
{
    if(n == 0) return 1;
    int rslt(1);
 
    while(n)
    {
        if(n & 1) // 若幂为奇数
        {
            rslt *= A;
        }
        A *= A;
        n >>= 1; // 右位移等价于除以2
    }
    return rslt;
}

在上面这段代码中，我们还考虑了奇数次幂出现的情况，关于这一点很容易理解，我就不加赘述了。
回到我们的文章开头，其实上面所说的只能算是个引子，接下来我想要介绍的是关于矩阵的快速求幂的方法。其实聪明的读者应该已经明白了，矩阵和整数的快速幂运算算法在代数上应该是等价的，矩阵也具备快速幂运算所必需的条件：结合律。那么我们的任务就是构造一个矩阵类，并实现该类的乘法运算，那基本就算是完事了，当然，根据矩阵乘法的性质，如果想要有乘幂的话，则矩阵必须为方阵。
下面给出矩阵类的代码：
（代码3）

class Matrix
{
public:
    int N; // 矩阵维数
    int** m; // 存储矩阵的二维数组
 
    Matrix(int n = 2)
    {
        m = new int*[n];
        for(int i(0); i < n; ++i)
        {
            m[i] = new int[n];
        }
        N = n;
        clear();
    }
 
    // 将矩阵清空为零矩阵
    void clear()
    {
        for(int i(0); i < N; ++i)
        {
            memset(m[i], 0, sizeof(int) * N);
        }
    }
 
    // 将矩阵设定为单位矩阵
    void unit()
    {
        clear();
        for(int i(0); i < N; ++i)
        {
            m[i][i] = 1;
        }
    }
 
    // 矩阵的赋值
    Matrix operator= (Matrix &othr)
    {
        Matrix(othr.N);
        for(int i(0); i < othr.N; ++i)
        {
            for(int j(0); j < othr.N; ++j)
            {
                m[i][j] = othr.m[i][j];
            }
        }
        return *this;
    }
 
    // 矩阵的乘法
    //!假设所有因子均为同阶方阵
    Matrix operator* (Matrix &othr)
    {
        Matrix rslt(othr.N);
        for(int i(0); i < othr.N; ++i)
        {
            for(int j(0); j < othr.N; ++j)
            {
                for(int k(0); k < othr.N; ++k)
                {
                    rslt.m[i][j] += m[i][k] * othr.m[k][j];
                }
            }
        }
        return rslt;
    }
};

代码中提到的单位矩阵就是主对角线全为1，其它元素全为0的矩阵，这个矩阵的性质类似于整数中的1，按代数的讲法，它属于矩阵的幺元。
有了矩阵类，我们就可以开始编写矩阵快速幂的算法了，直接看代码吧：
（代码4）

Matrix qMPow(Matrix &A, int n)
{
    Matrix rslt(A.N);
    rslt.unit();
    if(n == 0) return rslt;
    while(n)
    {
        if(n & 1) // 若幂为奇数
        {
            rslt = rslt * A;
        }
        A = A * A;
        n >>= 1; // 右位移等价于除以2
    }
    return rslt;
}

现在请你回过头来看代码2，是不是几乎完全一样？就这样，我们用同样的思想完成了矩阵的快速幂算法，至于矩阵的快速幂运算有何用途，那就请留意下下一篇文章吧~
这次就先到这了。

导航

公告