HDU 1575 Tr A 矩阵快速幂

Problem Description
A为一个方阵，则Tr A表示A的迹（就是主对角线上各项的和），现要求Tr(A^k)%9973。

Input
数据的第一行是一个T，表示有T组数据。
每组数据的第一行有n(2 <= n <= 10)和k(2 <= k < 10^9)两个数据。接下来有n行，每行有n个数据，每个数据的范围是[0,9]，表示方阵A的内容。

Output
对应每组数据，输出Tr(A^k)%9973。

Sample Input
2
2 2
1 0
0 1
3 99999999
1 2 3
4 5 6
7 8 9

Sample Output
2
2686

问题分析
挺好理解的矩阵快速幂问题，智障的我一开始直接理解为把矩阵的每个元素k 次方，忘得真快啊。顺便借地整理一下矩阵乘法的思路（示例：a*b=c）。
1. 常用的基本矩阵乘法，即由三重循环进行累加求和。
代码示例：

void multi(int a[m][n],b[n][p],c[m][p])
{
    int temp;
    for(int i=0;i<m;i++)
    {
         for(int j=0;j<p;j++)
         {
             temp=0;
             for(int k=0;k<n;k++)
                temp+=a[i][k]*b[k][j]);
             c[i][j]=temp;
         }
    }
}

2.对矩阵 a,b,c逐行进行访问，以提高访存效率。即对于矩阵 a 的每一行，第 j 列的元素分别和矩阵 b 的第 j 个元素相乘，然后对于每一列的结果进行累加求和，从而得到矩阵 c 第 i 行第 k 列的结果。
代码示例：

void multi(int a[m][n],b[n][p],c[m][p])
{
    int temp[3]={0};
    for(int i=0;i<m;i++)
    {
        for(int j=0;j<p;j++)
            temp[j]=0;
        for(int j=0;j<n;j++)
        {
            for(int k=0;k<p;k++)
            {
                temp[k]+=a[i][j]*b[j][k];
            }
        }
        for(int k=0;k<p;k++)
            c[i][k]=temp[k];
    }
}

代码实现

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
#define mod 9973
typedef __int64 ll;
ll t;
#define n 16
int m;
struct matrix
{
    int a[n][n];
}z0,z1;
matrix multi(matrix x,matrix y)
{
    matrix temp;
    memset(temp.a,0,sizeof(temp.a));
    for(int i=0;i<m;i++)
        for(int j=0;j<m;j++)
            for(int k=0;k<m;k++)
                temp.a[i][j]=(temp.a[i][j]+x.a[i][k]*y.a[k][j])%mod;
    return temp;
}
matrix quick_pow(ll p)
{
    memset(z1.a,0,sizeof(z1.a));
    for(int i=0;i<m;i++)
        z1.a[i][i]=1;
    while(p)
    {
        if(p%2==1)
            z1=multi(z1,z0);
        p/=2;
        z0=multi(z0,z0);
    }
    return z1;
}
int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%I64d",&m,&t);
        for(int i=0; i<m; i++)
            for(int j=0; j<m; j++)
                scanf("%d",&z0.a[i][j]);
        matrix result =quick_pow(t);
        int r=0;
        for(int i=0; i<m; i++)
            r+=result.a[i][i];
        printf("%d\n",r%mod);

    }
    return 0;
}