学步园

【题意】

给定一个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.