【题意】
给定一个n*n的矩阵A,求S=A+A^2+A^4+..+A^k(每个元素mod m)
【输入】
第一行n、m、k
接下来n行每行n个数字描述矩阵
【输出】
输出S矩阵
A的几次方快速幂可求
但求的是A^1加到A^k,k十分大
这个还是采取二分思想
先算出来A^1+A^2+..A^(k/2)
再总体乘以(I+A^(k/2)),若k为奇数再加上A^k即可
算矩阵乘法的时候先枚举k再枚举j,判断一下[i,k]是不是0,是0就跳过,可以快上许多,没有这个优化的时候程序是tle的……
program poj3233;
type
square=array [1..30,1..30] of longint;
var
n,m,k,i,j:longint;
ans,root,xxx:square;
function matrixadd (a,b:square):square;inline;
var
i,j:longint;
begin
for i:=1 to n do
for j:=1 to n do
if a[i,j]+b[i,j]>=m then
matrixadd[i,j]:=a[i,j]+b[i,j]-m
else
matrixadd[i,j]:=a[i,j]+b[i,j];
end;
function matrixmul (a,b:square):square;inline;
var
i,j,k:longint;
ans:square;
begin
fillchar(ans,sizeof(ans),0);
for i:=1 to n do
for k:=1 to n do
if a[i,k]>0 then
for j:=1 to n do
ans[i,j]:=(ans[i,j]+a[i,k]*b[k,j]) mod m;
exit(ans);
end;
function power (now:longint):square;
var
i,j:longint;
ans:square;
begin
ans:=xxx;
for i:=30 downto 1 do
begin
ans:=matrixmul(ans,ans);
if now or (1 shl (i-1))=now then ans:=matrixmul(ans,root);
end;
exit(ans);
end;
procedure quick (now:longint);
var
i,j:longint;
begin
if now=1 then
begin
ans:=root;
exit;
end;
quick(now div 2);
ans:=matrixadd(ans,matrixmul(ans,power(now div 2)));
if now and 1 = 1 then ans:=matrixadd(ans,power(now));
end;
begin
read(n,k,m);
for i:=1 to n do
for j:=1 to n do
begin
read(root[i,j]);
root[i,j]:=root[i,j] mod m;
end;
for i:=1 to n do
xxx[i,i]:=1;
quick(k);
for i:=1 to n do
begin
for j:=1 to n do
write(ans[i,j],' ');
writeln;
end;
end.