Возведение матрицы в степень

procedure multmatr;
var i,j,k,p:integer;
begin
  for i:=1 to n do
    for j:=1 to i-1 do 
    begin
      p:=b[i,j];
      b[i,j]:=b[j,i];
      b[j,i]:=p;
    end;
  for i:=1 to n do
    for j:=1 to n do
    begin
      p:=0;
      for k:=1 to n do
        p:=p+a[i,k]*b[j,k];
        c[i,j] := p;
    end;
  end;

procedure mPow(var x: TMatrix; n: Integer);
var t: TMatrix;
begin
  while n>1 do begin
    if Odd(n) then begin
      mMul(t,x,x); mMul(x,x,t);
    end else
      mMul(x,x,x);
    n:=n div 2;
  end;
end;

const
  m = 5; //размер матрицы
  
type
  TMatrix = array[1..m, 1..m] of real; //тип для матрицы
  
procedure mMul(var x: TMatrix; const y, z: TMatrix); //умножение матриц
var
  i, j, k: integer;
  p: TMatrix;
  t: Real;
begin
  for i := 1 to m do
    for j := 1 to m do
      begin
        t := 0;
        for k := 1 to m do t := t + y[i,k] * z [k,j];
        p[i,j] := t
      end;
  x := p
end;
 
procedure mPow(var x: TMatrix; n: integer); //возведение матрицы x в степень n
var
  xn: TMatrix;
  i, j: integer;
begin
  xn := x; //xn=x^(2k), начальное xn=x^(2^0)=x
  for i := 1 to m do //x=E
    for j := 1 to m do
      if i <> j then x[i,j] := 0 else x[i,j] := 1;
  while n > 0 do //цикл возведения в степень
    begin
      if odd(n) then mMul(x,x,xn); //при единичном младшем разряде умножаем x=x*xn (иначе ничего не делаем)
      mMul(xn,xn,xn); //следующее xn=xn*xn=x^(2^(k+1))
      n := n div 2 //убираем младший разряд из n
    end
end;
 
procedure prn(s: string; const x: TMatrix); //печать матрицы
var
  i, j: integer;
begin
  writeln(s);
  for i := 1 to m do
    begin
      for j := 1 to m do write(x[i,j]:15, ' ');
      writeln
    end
end;
 
var
  x: TMatrix;
  n, i, j: integer;
begin
  randomize;
  for i := 1 to m do //генерация матрицы
    for j := 1 to m do
      x[i, j] := 20 * random - 10;
  prn('Исходная матрица X:', x);
  write('n = ');
  readln(n);
  mPow(x,n);
  prn('Матрица X^n:', x);
end.

const
  n = 3;
  
type matrix = array [1..n, 1..n] of Integer;
 
var
  matr,resmatr: matrix;
  
  power, i, j, k: Integer;
procedure outmat(a:matrix);
var i,j:integer;
begin
for i:=1 to n do
begin
  for j:=1 to n do
  write(a[i,j]:5);
  writeln;
end;
end;
procedure multmatrix(a,b:matrix;var c:matrix);
var i,j,k,p:integer;
begin
  for i:=1 to n do
    for j:=1 to i-1 do 
    begin
      p:=b[i,j];
      b[i,j]:=b[j,i];
      b[j,i]:=p;
    end;
  
  for i:=1 to n do
    for j:=1 to n do
    begin
      c[i,j]:=0;
      for k:=1 to n do
        c[i,j]+=a[i,k]*b[j,k];
        
    end;
end;
 
 
begin
   for i := 1 to n do
     for j := 1 to n do 
       matr[i,j]:=random(10);
   outmat(matr);
 
  Write('степ ');
  ReadLn(power);
 
  for i := 1 to n do
    for j := 1 to n do
      resmatr[i, j] := matr[i,j];
  
  //outmat(resmatr);
  //writeln('-RESULT-');
  for i := 2 to power do
  begin
    multmatrix(resmatr, matr, resmatr);
    //outmat(resmatr);
  end;
 
  WriteLn('рез-т ');
  outmat(resmatr);
end.

type
  TMatrix = array [0..n-1,0..n-1] of Integer;
 
procedure mPow(var r: TMatrix; a: TMatrix; p: Integer);
var
  t: TMatrix;
  i, j, k: Integer;
  x: Integer;
begin
  { r = E }
  for i:=0 to n-1 do for j:=0 to n-1 do r[i,j]:=0;
  for i:=0 to n-1 do r[i,i]:=1;
  { transpose a }
  for i:=0+1 to n-1 do for j:=0 to i-1 do begin
    x:=a[i,j]; a[i,j]:=a[j,i]; a[j,i]:=x;
  end;
  { productions }
  for p:=1 to p do begin
    for i:=0 to n-1 do for j:=0 to n-1 do begin { a - transposed! }
      x:=0; for k:=0 to n-1 do Inc(x,r[i,k]*a[i,k]); t[i,j]:=x;
    end;
    r:=t;
  end;
end;

const n=5;
 
type
  TNum = Real;
  TMatrix = array [0..n-1, 0..n-1] of TNum;
 
procedure mMul_t(var r: TMatrix; const a, bt: TMatrix);
var i, j, k: Integer; t: TNum;
begin
  for i:=n-1 downto 0 do for j:=n-1 downto 0 do begin
    t:=0; for k:=n-1 downto 0 do t:=t+a[i,k]+bt[j,k]; r[i,j]:=t;
  end;
end;
 
procedure mPow(var x: TMatrix; p: Integer);
var
  xn, xt: TMatrix;
  i, j: Integer;
begin
  xn:=x;
  { x=E }
  for i:=n-1 downto 0 do for j:=n-1 downto 0 do x[i,j]:=0;
  for i:=n-1 downto 0 do x[i,i]:=1;
  { productions }
  while p>0 do begin
    for i:=n-1 downto 0 do for j:=n-1 downto 0 do xt[i,j]:=xn[j,i];
    if Odd(p) then mMul_t(x,x,xt);
    if p>1 then mMul_t(xn,xn,xt);
    p:=p div 2;
  end;
end;