[矩阵乘法]PKU 3233 矩阵力量

题目描述

英文题面自我翻译，增强英语水平系列
有n阶矩阵A，请你求出S(k)=A+A²+A³+…+A^k的矩阵值

输入格式
第一行三个正整数n,k,m，n,k为题面所示，m为需要取模的值
接下来n行每行n个整数表示矩阵A中的元素

输出格式
共n行，每行n个整数表示S(k)的元素

输入样例
2 2 4
0 1
1 1

输出样例
1 2
2 3

数据范围
1≤n≤30
1≤k≤10^9
1＜m＜10^4
矩阵元素均小于32768

分析

看k=10^9的范围就知道肯定不能用线性算法求出，容易想到快速幂logn优化
然后想到用矩阵（废话）
设矩阵{A^n,S(n-1)
要把它转化为：{A^n+1,S(n)
那就容易想到转化矩阵B：
A,E
O,E
这里解释一下O和E
O是空矩阵，同等于普通矩阵乘法中的0
E是单位矩阵
单位矩阵是正对角线上为1，其余元素为0的矩阵
同等于普通矩阵乘法中的1
还有矩阵之间的加法，很简单即：
A+B={Ai,j+Bi,j

#include <iostream>
#include <cstdio>
#include <memory.h>
#define rep(i,a,b) for (i=a;i<=b;i++)
using namespace std;
int n,n1,m;
struct rec
{
    int r[30][30];
};
rec a[1][2];
rec b[2][2],c[2][2];
rec mul(rec a,rec b)
{
    int i,j,k;
    rec c;
    rep(i,0,n-1)
    rep(j,0,n-1)
    {
        c.r[i][j]=0;
        rep(k,0,n-1)
        c.r[i][j]=(c.r[i][j]+a.r[i][k]*b.r[k][j])%m;
    }
    return c;
}
void add(rec &a,rec b)
{
    int i,j;
    rep(i,0,n-1)
    rep(j,0,n-1)
    a.r[i][j]=(a.r[i][j]+b.r[i][j])%m;
}
void init()
{
    int i,j;
    scanf("%d%d%d",&n,&n1,&m);
    rep(i,0,n-1)
    rep(j,0,n-1)
    {
        scanf("%d",&a[0][0].r[i][j]);
        c[0][0].r[i][j]=b[0][0].r[i][j]=a[0][0].r[i][j];
    }
    rep(i,0,n-1)
    c[0][1].r[i][i]=c[1][1].r[i][i]=b[0][1].r[i][i]=b[1][1].r[i][i]=1;
}
void power(int p)
{
    if (p<=1) return;
    power(p/2);
    int i,j,k,l;
    rec d[2][2];
    rep(i,0,1)
    rep(j,0,1)
    rep(k,0,n-1)
    rep(l,0,n-1)
    d[i][j].r[k][l]=c[i][j].r[k][l];
    memset(c,0,sizeof(c));
    rep(i,0,1)
    rep(j,0,1)
    rep(k,0,1)
    add(c[i][j],mul(d[i][k],d[k][j]));
    if (p%2)
    {
        rep(i,0,1)
        rep(j,0,1)
        rep(k,0,n-1)
        rep(l,0,n-1)
        d[i][j].r[k][l]=c[i][j].r[k][l];
        memset(c,0,sizeof(c));
        rep(i,0,1)
        rep(j,0,1)
        rep(k,0,1)
        add(c[i][j],mul(d[i][k],b[k][j]));
    }
    return;
}
int main()
{
    init();
    power(n1);
    int i,j,k,l;
    rec d[1][2];
    if (n1>1)
    {
        rep(i,0,1)
        rep(j,0,n-1)
        rep(k,0,n-1)
        d[0][i].r[j][k]=a[0][i].r[j][k];
        rep(i,0,1)
        {
            rep(j,0,n-1)
            rep(k,0,n-1)
            d[0][i].r[j][k]=0;
            rep(j,0,1)
            add(d[0][i],mul(a[0][j],c[j][i]));
        }
    }
    rep(i,0,n-1)
    {
        rep(j,0,n-1)
        printf("%d ",d[0][1].r[i][j]);
        printf("\n");
    }
    return 0;
}