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

1
2
3
4
5
6
7
8
9
10
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)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
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)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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,是不是几乎完全一样?就这样,我们用同样的思想完成了矩阵的快速幂算法,至于矩阵的快速幂运算有何用途,那就请留意下下一篇文章吧~
这次就先到这了。