[转载]快速幂与矩阵快速幂

来源：http://huanyouchen.github.io/2018/05/23/Quick-Matrix-Pow/

矩阵快速幂

矩阵的快速幂是用来高效地计算矩阵的高次方的。将朴素的o（n）的时间复杂度，降到log（n）。
本文先学习快速幂和矩阵乘法的基础知识，然后将两者结合实现矩阵快速幂方法。然后举一个例子：使用矩阵快速幂求斐波那契数列。

---

快速幂

一般计算底数x的n次幂${\mathrm{xnxn 的方法： xn=x\times x\times x\dots x\times xxn=x\times x\times x\dots x\times x ,需要做n次乘法运算，代码实现如下：}}^{}$

1 2 3 4 5 6 7

def powx(x, n): res = 1 for i in range(n): res = res * x return res   print(powx(x,n))

时间复杂度为O(N), 当n非常大时候运算效率很低。怎么才能提高运算效率来快速计算底数x的n次幂呢？可以通过快速幂方法。

先用一个具体的例子解释其原理：比如，计算x的11次方${\mathrm{x11x11，而11可以用二进制表示：11=1011=1\times 23+0\times 22+1\times 21+1\times 2011=1011=1\times 23+0\times 22+1\times 21+1\times 20}}^{}$

那么可以把${\mathrm{x11x11转换为：x11=x23\times x21\times x20=x8\times x2\times x1x11=x23\times x21\times x20=x8\times x2\times x1}}^{}$

原先需要做11次乘法运算，转换后只需要3次乘法运算。

通过上面的具体例子来推导一般情况计算${\mathrm{xnxn 的方法，先把n转换为2进制，从低位到高位根据二进制中的0或者1来进行乘法运算，比如上面的x11x11例子，从低位到高位的运算过程：}}^{}$

• 11 = 1011
• 1：res = ${\mathrm{x20x20}}^{}$
• 1: res = ${\mathrm{x20\times x21x20\times x21}}^{}$
• 0: 跳过，不运算
• 1：res = ${\mathrm{x20\times x21\times x23x20\times x21\times x23}}^{}$

得到结果res = ${\mathrm{x1\times x2\times x8=x11x1\times x2\times x8=x11}}^{}$

判断n的二进制低位是0或者1的方法, 也就是判断n是偶数或者奇数的方法，可以通过位运算and来实现：

• n and 1 返回1， 则n是奇数，即n的二进制低位是1
• n and 1 返回0， 则n是偶数，即n的二进制低位是0

从低位到高位的运算过程也可以通过位运算>>实现， n = n >> 1, 把n的二进制位右移过程也就是高位到低位的过程。比如11 = 1011的右移过程:

• n = n >> 1 = 1011 >> 1 = 101
• n = n >> 1 = 101 >> 1 = 10
• n = n >> 1 = 10 >> 1 = 1
• n = n >> 1 = 1 >> 1 = 0

根据上面分析，快速幂的代码实现：

1 2 3 4 5 6 7 8 9 10

def powx(x, n): res = 1 while n: if n&1: res = res * x x = x * x n = n >> 1 return res   print(powx(x,n))

上面快速幂代码中第5,6行涉及到两步乘法运算，当x很大时，比如${\mathrm{98765432119876543211,\; 这样直接乘效率也比较低，可以通过快速乘法进一步优化：}}^{}$

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

#快速乘法 def fast_multi(a,b): multi_res = 0 while b: if b&1: multi_res = multi_res + a a = a+a b >>= 1 return multi_res   #快速幂 def powx(x, n): res = 1 while n: if n&1: res = fast_multi(res, x) x = fast_multi(x, x) n = n >> 1 return res

---

矩阵乘法

矩阵：一个$\mathrm{m\times nm\times n的矩阵就是m\times nm\times n个数排成m行n列的一个数阵。}$

矩阵乘法：设A为$\mathrm{m\times pm\times p的矩阵，B为p\times np\times n的矩阵，那么称m\times nm\times n的矩阵C为矩阵A与B的乘积，记作C=AB。矩阵A与矩阵B相乘需要满足条件：A的列数等于B的行数。}$

矩阵乘法举例：

令$\mathrm{A=[1234]A=[1234] ， B=[5678]B=[5678]，\; 则：}$

$$\mathrm{C=AB=[1234]\times [5678]=[1\times 5+2\times 71\times 6+2\times 83\times 5+4\times 73\times 6+4\times 8]=[19224350]C=AB=[1234]\times [5678]=[1\times 5+2\times 71\times 6+2\times 83\times 5+4\times 73\times 6+4\times 8]=[19224350]}$$

代码实现：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

def MatrixMultiply(matrix_a, matrix_b): n_row = len(matrix_a) n_col = len(matrix_b[0]) # C的行数等于A的行数，C的列数等于B的列数 matrix_c = [[0 for j in range(n_col)] for i in range(n_row)] for i in range(0, n_row): for j in range(0, n_col): for k in range(0, n_row): matrix_c[i][j] = matrix_c[i][j] + matrix_a[i][k] * matrix_b[k][j] return matrix_c     a = [[1,2],[3,4]] b = [[5,6],[7,8]] c = MatrixMultiply(a, b) print(c) # [[19, 22], [43, 50]]

且矩阵乘法满足结合律：

$${\mathrm{A11=A8\times A2\times A1=[1234]8\times [1234]2\times [1234]1A11=A8\times A2\times A1=[1234]8\times [1234]2\times [1234]1}}^{}$$

---

矩阵快速幂

现在知道了两个矩阵乘法运算，那么如果求矩阵A的n次方${\mathrm{AnAn，可以用An=A\times A\times A\times \dots \times AAn=A\times A\times A\times \dots \times A，不过学习了计算底数xx的n次幂xnxn 的快速幂方法，把底数xx换成矩阵A，计算AnAn，使用上面快速幂的方法同样可以实现矩阵快速幂。}}^{}$

结合快速幂和矩阵乘法，实现代码：

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

def MatrixMultiply(matrix_a, matrix_b): n_row = len(matrix_a) n_col = len(matrix_b[0]) # C的行数等于A的行数，C的列数等于B的列数 matrix_c = [[0 for j in range(n_col)] for i in range(n_row)] for i in range(0, n_row): for j in range(0, n_col): for k in range(0, n_row): matrix_c[i][j] = matrix_c[i][j] + matrix_a[i][k] * matrix_b[k][j] return matrix_c     def getUnitMatrix(l): # 构造单位矩阵，l为矩阵a的行数 unit_matrix = [[0 for j in range(l)] for i in range(l)] for k in range(l): unit_matrix[k][k] = 1 return unit_matrix     def QuickMatrixPow(a, n): # res_matrix初始化为单位矩阵 res_matrix = getUnitMatrix(len(a)) while n: if n&1: res_matrix = MatrixMultiply(res_matrix, a) a = MatrixMultiply(a,a) n = n >> 1 return res_matrix     a = [[1,2],[3,4]] n = 11 print(QuickMatrixPow(a,n)) # [[25699957, 37455814], [56183721, 81883678]]

---

扩展：斐波那契数列矩阵算法

关于斐波那契数列的定义：

斐波那契数列（意大利语：Successione di Fibonacci），又译为费波拿契数列、费波那西数列、斐波那契数列、黄金分割数列。

在数学上，费波那契数列是以递归的方法来定义：

• $\mathrm{F(0)=0F(0)=0}$
• $\mathrm{F(1)=1F(1)=1}$
• $\mathrm{F(n)=F(n-1)+F(n-2)(n>=2)F(n)=F(n-1)+F(n-2)(n>=2)}$

用文字来说，就是费波那契数列由0和1开始，之后的费波那契系数就是由之前的两数相加而得出。首几个费波那契系数是：

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233……（OEIS中的数列A000045）

特别指出：0不是第一项，而是第零项。

维基百科:斐波那契数列

将$\mathrm{F(n)=F(n-1)+F(n-2)F(n)=F(n-1)+F(n-2)用矩阵表示：[F(n)=F(n-1)+F(n-2)F(n-1)]=[1110]\times [F(n-1)F(n-2)][F(n)=F(n-1)+F(n-2)F(n-1)]=[1110]\times [F(n-1)F(n-2)]}$

而$\mathrm{F(n-1)、F(n-2)、F(n-3)F(n-1)、F(n-2)、F(n-3)也可以用同样的矩阵表示方式，于是有：[F(n)F(n-1)]=[1110]\times [F(n-1)F(n-2)]=[1110]2\times [F(n-2)F(n-3)]=\dots =[1110]n-1\times [F(1)F(0)][F(n)F(n-1)]=[1110]\times [F(n-1)F(n-2)]=[1110]2\times [F(n-2)F(n-3)]=\dots =[1110]n-1\times [F(1)F(0)]}$

得到： $[F(n)F(n-1)]=[1110]n-1\times [F(1)F(0)][F(n)F(n-1)]=[1110]n-1\times [F(1)F(0)]$

我们知道$[F(1)=1F(0)=0][F(1)=1F(0)=0],且[1110]n-1[1110]n-1可以用上面介绍的矩阵快速幂计算出来，两者相乘得到矩阵的第1行第1列就是F(n)F(n)了。$

实现代码：

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

def MatrixMultiply(matrix_a, matrix_b): n_row = len(matrix_a) n_col = len(matrix_b[0]) # C的行数等于A的行数，C的列数等于B的列数 matrix_c = [[0 for j in range(n_col)] for i in range(n_row)] for i in range(0, n_row): for j in range(0, n_col): for k in range(0, n_row): matrix_c[i][j] = matrix_c[i][j] + matrix_a[i][k] * matrix_b[k][j] return matrix_c     def get_unit_matrix(l): unit_matrix = [[0 for j in range(l)] for i in range(l)] for k in range(l): unit_matrix[k][k] = 1 return unit_matrix     def QuickMatrixPow(a, n): res_matrix = get_unit_matrix(len(a)) while n: if n&1: res_matrix = MatrixMultiply(res_matrix, a) a = MatrixMultiply(a,a) n = n >> 1 return res_matrix     def get_Fib_n(n): if n == 0: return 0 elif n == 1: return 1 else: a = [[1,1],[1,0]] base = [[1],[0]] Fib_n = MatrixMultiply(QuickMatrixPow(a,n-1),base) return Fib_n[0][0]   if __name__ == '__main__': # 求斐波那契数列第F(n)个数 n = 8 print(get_Fib_n(n)) # 21

结束。

---------------------------------------

这里附上我个人的c++实现代码：

#include<bits/stdc++.h>
using namespace std;
typedef int** matrix;
/*此处简化，只计算方阵 */
/*矩阵乘法 */
matrix MatrixMultiply(matrix a,matrix b,int n){
    matrix c=new int*[n];
    for(int i=0;i<n;i++){
        c[i]=new int[n];
        memset(c[i],0,n*sizeof(int));
    }
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            for(int k=0;k<n;k++){
                c[i][j]=c[i][j]+(a[i][k])*(b[k][j]);
            }
        }
    }
    return c;
}
/*构造单位矩阵 */
matrix getUnitMatrix(int n){
        matrix c=new int*[n];
    for(int i=0;i<n;i++){
        c[i]=new int[n];
        memset(c[i],0,n*sizeof(int));
    }
    for(int k=0;k<n;k++){
        c[k][k]=1;
    }
    return c;
}
matrix QuickMatrixPow(matrix a,int n,int times){
    matrix c=getUnitMatrix(n);
    while(times){
        if(times&1)
        c=MatrixMultiply(c,a,n);
        a=MatrixMultiply(a,a,n);
        times>>=1;
    }
    return c;
}
int main(){
    matrix a=new int*[2];
    int times=11;
    a[0]=new int[2];
    a[1]=new int[2];
    a[0][0]=1;
    a[0][1]=2;
    a[1][0]=3;
    a[1][1]=4;
    matrix ans=QuickMatrixPow(a,2,times);
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            cout<<ans[i][j]<<' ';
        }
        cout<<endl;
    }
    getchar();
}