:: Pascal Page :: Массивы и матрицы

Массивы и матрицы

Произведение сцепленных матриц / Определитель матрицы
Макс. последовательности в массиве: одинаковые / возрастающие элементы
Как задать матрицу, чтобы быстро поменять местами ее строки / столбцы
Вычисление заданного многочлена от матрицы	Вычисление квадратного корня из матрицы

	Как заполнить матрицу "по спирали"

Произведение сцепленных матриц

Вычисление произведения "сцепленных" матриц,
Например:

(2, 3) Х (3, 3) = (2, 3)
(3, 2) Х (2, 5) = (3, 5)

const
	nRowFirst  = 3;	{ Количество строк первой матрицы }
	nCommon = 2;	{ Количество столбцов первой и строк второй матрицы }
	nColSecond = 5;	{ Количество столбцов второй матрицы }

type
	TType = Integer;
	TFirst = { Тип данных для хранения первой матрицы }
		array[1 .. nRowFirst, 1 .. nCommon] of TType;
	TSecond = { Тип данных для хранения второй матрицы }
		array[1 .. nCommon, 1 .. nColSecond] of TType;
	TResult = { Тип данных для хранения результата }
		array[1 .. nRowFirst, 1 .. nColSecond] of TType;

{ Собственно процедура перемножения матриц: (res) = (a) X (B) }
procedure matrixMult(var res: TResult;
	a: TFirst; b: TSecond);
var i, j, k: integer;
begin
	for i := 1 to nRowFirst do
		for j := 1 to nColSecond do begin
			res[i, j] := 0;
			for k := 1 to nCommon do
				res[i, j] := res[i, j] + a[i, k] * b[k, j];
		end
end;

var
	iRow, iCol: Integer;
	m1: TFirst; m2: TSecond;
	mr: TResult;
begin
	WriteLn('Ввод первой матрицы: (', nRowFirst, 'x', nCommon, ')');
	for iRow := 1 to nRowFirst do
		for iCol := 1 to nCommon do begin
			write('m1[', iRow, ',', iCol,'] : '); ReadLn(m1[iRow, iCol])
		end;

	writeln('Ввод второй матрицы: (', nCommon, 'x', nColSecond, ')');
	for iRow := 1 to nCommon do
		for iCol := 1 to nColSecond do begin
			write('m2[', iRow, ',', iCol,'] : '); ReadLn(m2[iRow, iCol])
		end;

	matrixMult(mr, m1, m2); { Вызов процедуры перемножения }

	writeln('Результирующая матрица: (', nRowFirst, 'x', nColSecond, ')');
	for iRow := 1 to nRowFirst do begin
		for iCol := 1 to nColSecond do
			write(mr[iRow, iCol]:4);
		writeln
	end;
end.

Определитель матрицы (рекурсивное определение)

Вычисление детерминанта (определителя) матрицы с помощью рекурсии:

const
	max_n = 4;
type
	matrix = array[1 .. max_n, 1 .. max_n] of real;

{ Матрица, для которой будет вычисляться определитель }
const 
	a: matrix = (
		(2, 9, 9, 4),
		(2, -3, 12, 8),
		(4, 8, 3, -5),
		(1, 2, 6, 4)
	);

function minusOne(n: integer): integer;
begin
	minusOne := (1 - 2*Byte(Odd(n)));
end;

function get_addr(i, j: integer; const n: integer): integer;
begin
	get_addr := pred(i) * n + j
end;

{ Рекурсивное определение определителя }
function det(var p; const n: integer): real;
type
	matrix = array[1 .. max_n * max_n] of real;
var
	my_p: matrix absolute p;
	pp: ^matrix;
	s: real;
	i, j, curr: integer;
begin
	s := 0.0;
	if n = 2 then begin
		det := my_p[1]*my_p[4] - my_p[2]*my_p[3]; exit
	end;
	
	for i := 1 to n do begin
		GetMem(pp, Sqr(Pred(n)) * sizeof(real));
		curr := 1;
		for j := 1 to n do
			if j <> i then begin
				move(my_p[get_addr(j, 2, n)],
					pp^[get_addr(curr, 1, Pred(n))], pred(n) * sizeof(real));
				inc(curr);
			end;
		s := s + minusOne(Succ(i)) * my_p[get_addr(i, 1, n)] * det(pp^, Pred(n));
		FreeMem(pp, Sqr(Pred(n)) * sizeof(real))
	end;
	det := s
end;

begin
	writeln( det(a, 4):0:0 );
end.

Как найти в массиве максимальную последовательность одинаковых символов?

Эта функция возвращает индекс первого элемента найденной последовательности. Отсчет ведется с нулевого элемента... Длина цепочки возвращается через параметр Len). Если в массиве несколько последовательностей одинаковых элементов, возвращается начальный индекс и длина первой из них.

type TType = integer;
const
	n = 10;
	a: array[1 .. n] of TType = (1, 2, 2, 3, 4, 5, 5, 5, 26, 6);

function get_chain(var arr: array of TType;
		const n: integer; var Len: integer): integer;
var
	i, start, count: integer;
	T: TType;
begin
	Len := 0; i := 0;
	while i <= pred(n) do begin
		count := 1; T := arr[i];
		while arr[i] = T do begin
			inc(count); inc(i)
		end;
		
		if Len < count then begin
			start := i - count;
			Len := count
		end;
		inc(i)
	end;
	get_chain := start
end;

{ Пример использования: }
var _start, _max: Integer;
begin
	_start := get_chain(a, n, _max);
	writeln('start = ', _start, ' len = ', _max)
end.

Как найти в массиве максимальную возрастающую последовательность символов?

Эта функция возвращает индекс первого элемента найденной последовательности. Отсчет ведется с нулевого элемента... Длина цепочки возвращается через параметр max). Если в массиве несколько возрастающих последовательностей элементов, возвращается начальный индекс и длина первой из них.

type TType = integer;
const
	n = 15;
	a: array[1 .. n] of TType = (1, 2, 5, 6, 3, 7, 8, 9, 8, 11, 23, 33, 45{0}, 56, 90);

function find_ascend(var arr: array of ttype;
		const n: integer; var max: integer): integer;
var curr, start: integer;

	function check_max(i: integer): integer;
	begin
		if max < curr then begin
			max := curr;
			check_max := i - curr + byte(i = pred(n))
		end
		else check_max := start
	end;

var i: integer;
begin
	max := 0; curr := 1;
	start := 1;

	for i := 1 to pred(n) do begin
		if arr[i - 1] < arr[i] then inc(curr)
		else begin
			start := check_max(i);
			curr := 1
		end;
	end;
	start := check_max(i);
	find_ascend := start
end;

{ Пример использования: }
var _start, _max: integer;
begin
	_start := find_ascend(a, n, _max);
	writeln('start = ', _start, ' len = ', _max)
end.

Как задать матрицу, чтобы быстро поменять местами ее строки ?

Для того, чтобы поменять местами строки в матрице максимально быстро, матрицу нужно задавать как массив строк:

const
	size_row = 10;	{ число строк }
	size_col = 10;	{ число столбцов }
type
	tvector = array[1 .. size_col] of integer;
	tmatrix = array[1 .. size_row] of tvector;

При таком определении для того, чтобы поменять местами строки матрицы, достаточно сделать следующее:

procedure swap_rows(var mx: tmatrix; const i, j: integer);
var T: tvector;
begin
	T := mx[i]; mx[i] := mx[j]; mx[j] := T
end;

procedure print(var mx: tmatrix);
var i, j: integer;
begin
	for i := 1 to size_row do begin
		for j := 1 to size_col do
			write(mx[i][j]:4);
		writeln
	end;
end;

var
  mx: tmatrix;
  i, j: integer;

begin
	{ Заполнение матрицы }
	for i := 1 to size_row do
		for j := 1 to size_col do
			mx[i][j] := random(100);

	{ Матрица до обмена }
	writeln('before:'); print(mx);

	{ Обмен строк }
	swap_rows(mx, 1, 5);

	{ Матрица после обмена }
	writeln('after:'); print(mx);
end.

Как задать матрицу, чтобы быстро поменять местами ее столбцы ?

Для того, чтобы поменять местами столбцы в матрице максимально быстро, матрицу нужно задавать как массив столбцов:

const
	size_row = 10;	{ число строк }
	size_col = 10;	{ число столбцов }
type
	trow = array[1 .. size_row] of integer;
	tmatrix = array[1 .. size_col] of trow;

При таком определении для того, чтобы поменять местами столбцы матрицы, достаточно сделать следующее:

procedure swap_cols(var mx: tmatrix; const i, j: integer);
var T: trow;
begin
	T := mx[i]; mx[i] := mx[j]; mx[j] := T
end;

procedure print(var mx: tmatrix);
var i, j: integer;
begin
	for i := 1 to size_row do begin
		for j := 1 to size_col do
			write(mx[j][i]:4); { <<< Обратите внимание на индексы !!! }
		writeln
	end;
end;

var
	mx: tmatrix;
	i, j: integer;
begin
	{ Заполнение матрицы }
	for i := 1 to size_row do
		for j := 1 to size_col do
			mx[i][j] := random(100);

	{ Матрица до обмена }
	writeln('before:'); print(mx);

	{ Обмен столбцов }
	swap_cols(mx, 1, 5);

	{ Матрица после обмена }
	writeln('after:'); print(mx);
end.

Как вычислить заданный многочлен от матрицы A

Задан порядок квадратной матрицы A, ее элементы и коэффициенты многочлена (a, b, c, ...). Необходимо вычислить заданный многочлен от матрицы A:

P(A) = a*A^2 + b*A + c (пример многочлена 2-го порядка)

{ Порядок матрицы }
const size = 4;
type
	TMatrix = array[1 .. size, 1 .. size] of real;

(* Умножение матриц *)
procedure matrixMult(Var m: TMatrix; a, b: TMatrix);
var i, j, k: Integer;
begin
	for i := 1 to size do
		for j := 1 to size do begin
			m[i, j] := 0;
			for k := 1 to size do
				m[i, j] := m[i, j] + a[i, k] * b[k, j]
		end;
end;

(* Возведение матрицы в степень *)
procedure matrixPower(var m: TMatrix; a: TMatrix; pow: integer);
var
	i, j: Integer;
	T: TMatrix;
begin
	if pow = 0 then begin
		for i := 1 to size do
			for j := 1 to size do m[i, j] := Byte(i = j);
		exit
	end;

	move(a, T, sizeof(T));
	for i := 1 To Pred(pow) do
		matrixMult(T, T, a);
	move(T, m, sizeof(T))
end;

(* Сложение матриц *)
procedure matrixAdd(var m: TMatrix; a, b: TMatrix);
var
	i, j: Integer;
begin
	for i := 1 to size do
		for j := 1 to size do
			m[i, j] := a[i, j] + b[i, j]
end;

(* Умножение матрицы на число *)
procedure matrixScale(var m: TMatrix; a: TMatrix; f: real);
var
	i, j: Integer;
begin
	for i := 1 to size do
		for j := 1 to size do
			m[i, j] := f * a[i, j]
end;

(* Печать матрицы *)
procedure matrixPrint(a: TMatrix);
var
	i, j: Integer;
begin
	for i := 1 to size do begin
		for j := 1 to size do
			write(a[i, j]:9:2);
		writeln
	end
end;

const
	n = 3; { Порядок многочлена }
	p: array[1 .. n] of real = { Коэффициенты многочлена }
		(1.0, -2.0, 3.0);
var
	a: TMatrix;
	T, Res: TMatrix;
	i, j: Integer;
begin
	for i := 1 to size do
		for j := 1 to size do
			a[i, j] := Random(20);
	matrixPrint(a);

	for i := 1 to n do begin
		matrixPower(T, a, n - i);
		matrixScale(T, T, p[i]);
		matrixAdd(Res, Res, T)
	end;
	matrixPrint(Res)
end.

Пример вычисления этого же задания с использованием Free Pascal Compiler-а и перегрузки функций приведен здесь:
Перегрузка операций - пример использования

Вычисление квадратного корня из матрицы

Пусть нам задана квадратная матрица, и необходимо вычислить ее корень (то есть, матрицу, при возведении в квадрат дающую исходную).

Для того, чтобы это сделать, воспользуемся Denman–Beavers iteration.

Если Y₀ - исходная матрица, а Z₀ - единичная матрица того же порядка, то итерация заключается в следующем:

Y_k+1 = (1/2)*(Y_k + Z_k^-1)
Z_k+1 = (1/2)*(Z_k + Y_k^-1)
при этом матрица Y_k сходится к корню из Y₀, а матрица Z_k - к обратной ей матрице, т.е. к Y₀^-1/2

Прграммная реализация выглядит так (компилятор FPC 2.2.0 - для того, чтобы иметь возможность перегружать операции работы с матрицами - сложение, умножение - но в принципе эта программа довольно просто портируется и на Турбо Паскаль):

{$mode objfpc}

const
  size = 3;
  eps = 0.001;

type
  TMatrix = array[1 .. size, 1 .. size] of double;

const
  E: TMatrix = (
    (1, 0, 0),
    (0, 1, 0),
    (0, 0, 1)
  );


operator * (const a, b: TMatrix) m: TMatrix;
var i, j, k: integer;
begin
  for i := 1 to size do
    for j := 1 to size do begin
      m[i, j] := 0;
      for k := 1 to size do
        m[i, j] := m[i, j] + a[i, k] * b[k, j]
    end;
end;

operator * (const a: TMatrix; const f: double) m: TMatrix;
var i, j: integer;
begin
  for i := 1 to size do
    for j := 1 to size do
      m[i, j] := f * a[i, j]
end;

operator + (const a, b: TMatrix) m: TMatrix;
var i, j: integer;
begin
  for i := 1 to size do
    for j := 1 to size do
      m[i, j] := a[i, j] + b[i, j]
end;

operator ** (const a: TMatrix; const pow: integer) m: TMatrix;
var i, j: Integer;
begin
  if pow = 0 then begin
    for i := 1 to size do
      for j := 1 to size do
        m[i, j] := Byte(i = j);
    exit
  end;

  m := a;
  for i := 1 to pred(pow) do
    m := m * a;
end;

procedure matrixPrint(a: TMatrix);
var i, j: integer;
begin
  for i := 1 to size do begin
    for j := 1 to size do
      write(a[i, j]:9:2);
    writeln
  end
end;

function matrixInvert(a: TMatrix): TMatrix;

  function sign(r: double): shortint;
  begin
    if r >= 0 then sign := 1 else sign := -1;
  end;

  procedure AddStrings(var a, b: TMatrix;
            i1, i2: integer; r: double);
  var j: integer;
  begin
    for j := 1 to size do begin
      a[i1,j] := a[i1,j] + r*a[i2,j];
      b[i1,j] := b[i1,j] + r*b[i2,j];
    end;
  end;

  procedure MultString(var a, b: TMatrix;
            i1: integer; r: double);
  var j: integer;
  begin
    for j := 1 to size do begin
      a[i1,j] := a[i1,j]*r;
      b[i1,j] := b[i1,j]*r;
    end;
  end;


var
  i, j: integer;
  b: TMatrix;
begin
  b := E;

  for i :=1 to size do begin
    for j := i + 1 to size do
      AddStrings(a, b, i, j,sign(a[i,i])*sign(a[j,i]));

    if abs(a[i,i]) > Eps then begin
      multString(a, b, i, 1 / a[i,i]);
      for j := i+1 to size do
        AddStrings(a, b, j, i, -a[j,i]);
    end
    else begin
      writeln('error inverting matrix.'); halt(100);
    end
  end;

  if (a[size, size]>eps) then begin
    for i := size downto 1 do
      for j := 1 to i-1 do begin
        AddStrings(a, b, j, i, -a[j,i]);
      end;
  end
  else begin
    writeln('error inverting matrix.'); halt(100);
  end;

 result := b;
end;

const
  a: TMatrix = (
    (3, 2, 1),
    (1, 0, 2),
    (4, 1, 3)
  );

var
  Ypred, Y, Zpred, Z: TMatrix;
  i: integer;

begin
  Ypred := a; Zpred := E;

  for i := 1 to 4 do begin
    Y := (Ypred + matrixInvert(Zpred)) * 0.5;
    Z := (Zpred + matrixInvert(Ypred)) * 0.5;
    Ypred := Y; Zpred := Z;

    writeln('iteration #', i);
    matrixPrint(Ypred ** 2); { <--- для проверки печатаем возведенную в квадрат Ypred }
  end;

  writeln('result:');
  matrixPrint(Ypred);

end.

После запуска программы получаем вот такой результат:

iteration #1
     5.50     2.75     3.00
     3.25     1.25     2.75
     8.25     3.25     5.50
iteration #2
     3.37     2.26     1.15
     1.15    -0.01     2.26
     4.55     1.15     3.37
iteration #3
     3.01     2.01     1.00
     1.01    -0.00     2.01
     4.02     1.01     3.01
iteration #4
     3.00     2.00     1.00
     1.00     0.00     2.00
     4.00     1.00     3.00
result:
     1.74     0.97     0.02
    -0.05     0.32     0.97
     1.15    -0.05     1.74

, то есть, за 4 итерации был найден квадратный корень из матрицы...

Как заполнить матрицу "по спирали"

Начиная с левого верхнего угла, по часовой стрелке... Матрица может быть как квадратной, так и прямоугольной.

const
  m = 5; { Количество строк    }
  n = 8; { Количество столбцов }
type
  TMx = array[1 .. m, 1 .. n] of integer;

procedure Print(mx : TMx);
var i, j: byte;
begin
  writeln;
  for i := 1 to m do begin
    writeln;
    for j := 1 to n do write(mx[i,j]:3);
 end;
 writeln;
end;

procedure Spiral(var mx : TMx);
var
  i, j, c: byte;
  A: integer;
begin
   i := 1; j := 1; c := 0; A := 1;
   repeat
       while (A <= m*n) and (j <= n - c) do begin mx[i,j] := A; inc(j); inc(A); end;
       inc(i); dec(j);
       while (A <= m*n) and (i <= m - c) do begin mx[i,j] := A; inc(i); inc(A); end;
       dec(j); dec(i);
       while (A <= m*n) and (j >= 1 + c) do begin mx[i,j] := A; dec(j); inc(A); end;
       inc(c); inc(j); dec(i);
       while (A <= m*n) and (i >= 1 + c) do begin mx[i,j] := A; dec(i); inc(A); end;
       inc(j); inc(i);
   until A > m * n;
end;

var
  mx : TMx;
begin
  Spiral(mx);
  print(mx);
end.