矩陣，是線性代數中的基本概念之一。一個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相契合。
再以C=A^25模擬手算的情況驗證一下。
k=25, k!=0 C=C*A （運算結果 C=A） A=A*A （運算結果 A=A^2）
k=k/2=12!=0 12%2==0 A=A*A （運算結果A=A^4）
k=k/2=6 !=0 6%2==0 A =A*A （運算結果A=A^8）
k=6/2=3!=0 3%2==1 C=C*A （運算結果 C=A*A^8=A^9） A=A*A （運算結果 A=A^16）
k=3/2=1!=0 1%2==1 C=C*A （運算結果 C=A^9*A^16=A^25）
k=1/2=0 算法結束。 正好與25的二進制數 11001相契合。
【例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;
}
【例3】矩陣的跡。
設A是一個n*n矩陣，Tr（A）表示矩陣A的跡（就是主對角線上各項的和）。輸入n(2 <= n <= 10)、k(2 <= k < 10^9)和矩陣A的n*n個元素，求Tr(A^k)%9973的值。
（1）編程思路。
由於k值較大，因此求A^k直接采用矩陣快速冪運算。
（2）源程序。
#include <stdio.h> #include <string.h> #define MOD 9973 struct Matrix { int mat[11][11]; // 存儲矩陣中各元素 }; Matrix matMul(Matrix a ,Matrix b,int n) { Matrix c; memset(c.mat,0,sizeof(c.mat)); int i,j,k; for (k = 1; k<=n ; k++) for (i=1 ;i<=n ; i++) if (a.mat[i][k]!=0) for (j = 1 ;j<=n ;j++) c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % MOD; 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 t,n,k,sum,i,j; Matrix p; scanf("%d" ,&t); while (t--) { scanf("%d%d",&n,&k); for (i=1;i<=n;i++) for (j=1;j<=n;j++) scanf("%d",&p.mat[i][j]); p = quickMatPow(p,n,k); sum=0; for (i=1;i<=n;i++) sum=(sum+p.mat[i][i])%MOD; printf("%d\n" ,sum); } return 0; }
將此源程序提交給 HDU1575 ”Tr A“，可以Accepted。