矩阵乘法（一）：基本运算

      矩阵，是线性代数中的基本概念之一。一个m×n的矩阵就是m×n个数排成m行n列的一个数阵。在计算机中，一个矩阵实际上就是一个二维数组。因此，可以将矩阵定义为一个结构体：

      struct matrix
      {
            int  mat[110][110]; // 存储矩阵中各元素
            int  row,col; // 矩阵的大小，row行，col列
      };

     矩阵相乘是矩阵的一种基本运算。

     设a为m×n矩阵，b为n×k矩阵，则它们的乘积ab（有时记做a·b）是一个m×k矩阵。

     其乘积矩阵a·b的第i行第j列的元素为第一个矩阵a第i行上的n个数与第二个矩阵b第j列上的n个数对应相乘后所得的n个乘积之和。即：

 

      需要注意的是：只有当矩阵a的列数与矩阵b的行数相等时，矩阵a×b才有意义。因此，矩阵相乘不满足交换律。设a是3×4矩阵，b是4×5矩阵，a与b相乘后，a·b是3×5矩阵；但b·a根本就无法运算。

      矩阵乘法满足结合律。

【例1】矩阵的乘法。
      输入矩阵a和矩阵b的数据，输出新的矩阵c=a*b。

      例如，样例输入

4 3
1 2 3
4 5 6
7 8 9
10 11 12
3 5
7 8 9 10 11
4 5 6 7 8
1 2 3 4 5
样例输出
18 24 30 36 42
54 69 84 99 114
90 114 138 162 186
126 159 192 225 258

      （1）编程思路。

       按照矩阵乘法的定义，用一个三重循环完成运算。

      （2）源程序。

#include <stdio.h>
#include <string.h>
struct matrix
{
      int mat[110][110]; // 存储矩阵中各元素
      int row,col; // 矩阵的大小，row行，col列
};
matrix matmul(matrix a ,matrix b)      // 矩阵a*b
{
      matrix c;
      c.row=a.row;
      c.col=b.col;
      memset(c.mat,0,sizeof(c.mat));
      int i,j,k;
      for (i = 0; i<=a.row ; i++)
         for (j=0 ;j<b.col; j++)
              for (k = 0 ;k<a.col;k++)
                   c.mat[i][j] += a.mat[i][k] * b.mat[k][j];
      return c;
}
int main()
{
      int i,j,x,y;
      matrix a,b,c;
      scanf("%d%d",&x,&y);
      a.row=x;
      a.col=y;
      for (i=0;i<x;i++)
          for (j=0;j<y;j++)
              scanf("%d" ,&a.mat[i][j]);
      scanf("%d%d",&x,&y);
      b.row=x;
      b.col=y;
      for (i=0;i<x;i++)
           for (j=0;j<y;j++)
               scanf("%d" ,&b.mat[i][j]);
      c=matmul(a,b);
      for (i = 0 ;i <c.row;i++)
      {
           for (j=0;j<c.col;j++)
               printf("%5d" ,c.mat[i][j]);
           printf("\n");
      }
      return 0;
}

      在实际应用中，我们经常会用到矩阵的幂运算。

      n个矩阵a相乘称为a的n次方，或称a^n为矩阵a的n次幂。

      求矩阵a的n次方通常采用快速幂运算。下面我们来探讨快速幂运算的思路。

      由于矩阵乘法具有结合律，因此 a^4 = a * a * a * a = (a*a) * (a*a) = a^2 * a^2。由此可以得到这样的结论：当n为偶数时，a^n = a^(n/2) * a^(n/2)；当n为奇数时，a^n = a^(n/2) * a^(n/2) * a （其中n/2取整）。这样，我们可以采用一种类似于二分的思想快速求得矩阵的幂。

       例如，a^9 = a*a*a*a*a*a*a*a*a     （一个一个乘，要乘9次）

                        = a*(a*a)*(a*a)*(a*a)*(a*a)

                        = a*(a^2)^4  

                        = a*((a^2)^2)^2   （a平方后，再平方，再平方，再乘上剩下的一个a，要乘4次）

       设c=a^k，c初始化为一个单位矩阵，即c矩阵中除了对角线的元素为1外，其余全部元素为0。

       c.mat[i][i]=1 ，          c,mat[i][j]=0  (i!=j)。

       任何一个矩阵乘以单位矩阵就是它本身。即可以把单位矩阵等价为整数1。因此，矩阵快速幂的算法描述为：

       while (k!=0)

       {

               if  (k%2==1)  c=c*a;     // c=c*a，表示矩阵c与a相乘，结果送c

               a=a*a;

               b=b/2;

       }

       为加深理解，以c=a^9模拟手算一下。

       k=9,  k!=0                       c=c*a （运算结果 c=a）   a=a*a  （运算结果 a=a^2）

       k=k/2  k=4!=0  4%2==0                                              a=a*a   （运算结果a=a^4）

       k=k/2 k=2!=0  2%2==0                                               a =a*a   （运算结果a=a^8） 

       k=2/2 k=1!=0  1%2==1   c=c*a （运算结果 c=a*a^8=a^9）   a=a*a  （运算结果 a=a^16）

       k=1/2  k=0  算法结束。

       可以看出，上述手算过程正好和9的二进制数表示1001相契合。

 

【例2】矩阵快速幂。

      给定n*n的矩阵a，求a^k。

      输入格式：
      第1行， n,k

      第2至n+1行，每行n个数，第i+1行第j个数表示矩阵第i行第j列的元素

      输出格式：
      共n行，每行n个数，第i行第j个数表示矩阵第i行第j列的元素，每个元素模10^5+7

      （1）编程思路。

      因为矩阵的幂运算参与运算的矩阵一定是n*n方阵。因此，在下面的程序中我们将结构体定义简化，去掉表示矩阵行列的变量row和col。

      另外，矩阵乘法运算后，所得结果通常会很大，所以一般采用64位整数表示。同时一般会在计算过程中不断取模，避免高精度运算。

      （2）源程序。

#include <stdio.h>
#include <string.h>
#define modnum 100007
struct matrix
{
      __int64 mat[101][101]; // 存储矩阵中各元素
};
matrix matmul(matrix a ,matrix b,int n)
{
      matrix c;
      memset(c.mat,0,sizeof(c.mat));
      int i,j,k;
      for (i = 1; i<=n ; i++)
          for (j=1 ;j<=n ; j++)
             for (k = 1 ;k<=n ;k++)
             {
                   c.mat[i][j]=(c.mat[i][j]+a.mat[i][k] * b.mat[k][j]) % modnum;
             }
      return c;
}
matrix quickmatpow(matrix a ,int n,int b)   // n阶矩阵a快速b次幂
{
      matrix c;
      memset(c.mat ,0 ,sizeof(c.mat));
      int i;
      for (i = 1 ;i <= n ;i++)
            c.mat[i][i] = 1;
      while (b!=0)
      {
           if (b & 1)
                c = matmul(c ,a ,n);      // c=c*a;
           a = matmul(a ,a ,n);          // a=a*a
           b /= 2;
      }
      return c;
}
int main()
{
      int i,j,n,k;
      matrix a;
      scanf("%d%d",&n,&k);
      for (i=1;i<=n;i++)
           for (j=1;j<=n;j++)
               scanf("%i64d",&a.mat[i][j]);
      a=quickmatpow(a,n,k);
      for (i = 1 ;i <=n;i++)
      {
            for (j=1;j<=n;j++)
                 printf("%i64d ",a.mat[i][j]);
            printf("\n");
      }
      return 0;
}