Решение практических задач на Паскале. Выпуск 6

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

Преамбула

На нашем форуме (http://www.yourpascal.com/viewforum.php?f=6) был задан вопрос о возведении матрицы в степень в неподходящем для этого месте "нужно напечатать период дроби (на паскале)". Отвечу в рассылке, в которой давно уже ничего не писал. Увы! Времени очень мало

В этом выпуске я представлю способ решения задачи с использованием рекурсии, так как давно искал повод для этого :)). Вот и пришло время ...

Сейчас будет предложен вариант решения, в котором ранг матрицы фиксирован. Он задается константой max в модуле uMatrics. Если Вам захочется изменить ранг матрицы, то нужно изменить значение константы и заново откомпилировать программу. Это не очень удобно и в Delphi для решения такой проблемы существуют динамические массивы (которых в Паскале нет :(( ).

Но если нельзя, но хочется, то добрый старый Паскаль позволяет сделать все ... Правда, для этого матрицу следует расположить в динамической памяти компьютера и работать с ней с помощью указателя. Чтобы было легче понять более сложный вариант, я и в этот раз буду размещать матрицу в динамической памяти. Это вовсе не обязательно. Можно обойтись "традиционными" двумерными массивами, объявленными в разделе VAR программы. Но используя указатели на массивы, размещенные в динамической памяти, позволяет использовать функции, а не процедуры. А это тоже плюс

Если описанный здесь способ покажется сложным, то напишите мне - я опишу и тот, более простой, способ решения с "обычными" массивами-матрицами.

Краткое описание

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

Достоинство метода рекурсий в том, что код получается маленький (и, вернее, он гусарский).  Недостаток в том, что неэффективно используется память и при возведении в большую степень может не хватить стека. Тогда стек нужно увеличить за счет кучи. Как это делается, можно посмотреть в стандартной справке, если набрать $M и вызвать контекстную помощь, нажав сочетание Ctrl+F1.

Замечу, что испытания моей программы показали: переполнение стека наблюдается при возведении матрицы в степень больше 650. А для меньших степеней гораздо вероятнее получить ошибку переполнения. Она  возникает, когда в переменную типа Real делается попытка записать значение больше допустимого.

Сначала приведу исходный код программы, а потом дам пояснения. Хотите - читайте ...

Исходный код программы

Все основное находится в модуле uMatrics. Он показан в листинге 1.

Для тестирования используется программа Solve006.pas, текст которой приведен в листинге 2.

{Листинг 1}
unit uMatrix;
interface

CONST
 max = 2;
TYPE
  PMatrix = ^TMatrix;
  TMatrix = array [1..max, 1..max] of Real;

function NewMatrix: PMatrix;
function PowMatrix(Source: PMatrix; pow: Integer): PMatrix;
procedure ShowMatrix(matrix: PMatrix; Dest: String; Size, Digs: Byte; mes: String);
procedure AutoFillMatrix(var Matrix: PMatrix);
function MulMatrix(source1, source2: PMatrix): PMatrix;
function UnitaryMatrix: PMatrix;

implementation

{Возведение матрицы в степень pow}
function PowMatrix(Source: PMatrix; pow: Integer): PMatrix;
var p: PMatrix;
begin
 p:=NewMatrix;
 case pow of
   -MaxInt..-1: Exit;
    0: p:=UnitaryMatrix;
    else
       p:=MulMatrix(Source, PowMatrix(Source, Pow-1));
 end;
 PowMatrix:=p
end;

{Локальная подпрограмма для форматированного вывода числа типа Real
 При Size <= 0 вывод не форматируется
 При Dig < 0 используется только Size}
function RealToStrFmt(Value: Real; Size, Digs: Byte): String;
var s: String;
begin
 if Size <= 0 then Str(Value, s) else
    if Digs < 0 then Str(Value:Size, s) else Str(Value:Size:Digs, s);
 RealToStrFmt:=s
end;

procedure ShowMatrix(matrix: PMatrix; Dest: String;
     Size, Digs: Byte; mes: String);
var i, j: Integer;
    f: Text;
begin
  if Matrix = nil then Exit;
  if Dest <> '' then
  begin
     Assign(f, Dest);
     {$I-} Reset(f); {$I+}
     if IOResult <> 0 then begin
       {$I-} Rewrite(f); {$I+}
       if IOResult <> 0 then
       begin
         WriteLn('Невозможно создать файл ', DEST);
         WriteLn('Вывод будет происходить на экран');
         Dest:='';
       end;
     end
     else
     begin
       Close(f);
       Append(f);
       WriteLn(f);
       WriteLn('Результаты будут дописаны в конец файла ', Dest)
     end;
  end;
  if Dest = '' then WriteLn(mes) else WriteLn(f, mes);
  for i:=1 to max do begin
   for j:=1 to max do
     if Dest = '' then
       Write(RealToStrFmt(Matrix^[i,j],Size,Digs),' ')
     else
       Write(f, RealToStrFmt(Matrix^[i,j],Size,Digs),' ');
  if Dest = '' then
     WriteLn
  else
     WriteLn(f);
  end;
  if Dest <> '' then
  begin
    Flush(f);
    Close(f);
  end;
end;

procedure AutoFillMatrix(var Matrix: PMatrix);
var i, j: Integer; k: Integer;
begin
  if Matrix = nil then Matrix := NewMatrix;
  k:=1;
  for i:=1 to max do
   for j:=1 to max do
   begin Matrix^[i,j]:= -k; inc(k) end;
end;

{Перемножение матриц}
function MulMatrix(source1, source2: PMatrix): PMatrix;
var i, j, k: Integer;
     Dest: PMatrix;
begin
   New(Dest);
   for i:=1 to max do
     for j:=1 to max do begin
       Dest^[i, j]:=0;
       for k:=1 to max do
          Dest^[i, j] := Dest^[i, j] + Source1^[i,k]*Source2^[k,j]
     end;
   MulMatrix := Dest;
{   ShowMatrix(Dest,'');}
end;

function UnitaryMatrix: PMatrix;
var i, j: Integer;
    p: PMatrix;
begin
  New(p);
  for i:=1 to max do
    for j:=1 to max do
      if i=j then p^[i,j]:=1 else p^[i,j]:=0;
  UnitaryMatrix:=p
end;

function NewMatrix: PMatrix;
var p: PMatrix;
begin  New(p); NewMatrix:=p end;

{"Служебная" часть модуля для обеспечения
  правильной работы с динамической памятью}
function HeapFunc(Size: Word): Integer; far;
begin HeapFunc:=1 end;

VAR
   varForHeapRelease: Pointer;
   SaveExit: Pointer;

procedure FreeAllMemory; far;
begin
   Release(varForHeapRelease);
   ExitProc:=SaveExit;
end;

BEGIN
  Mark(varForHeapRelease);
  HeapError:=@HeapFunc;
  SaveExit:=ExitProc;
  ExitProc:=@FreeAllMemory;
END.

{Листинг 2}
uses uMatrix;
VAR
   m: PMatrix;
BEGIN
  WriteLn;
  AutoFillMatrix(m);
  ShowMatrix(m, 'test.txt',12, 2, 'Исходная матрица');
  m:=PowMatrix(m, 2);
  ShowMatrix(m,'test.txt', 12, 2, 'Матрица после возведения в степень 2');
END.

Пояснения

Как отмечено вначале, здесь использован самый дубовый прямой метод: матрица каждый раз умножается сама на себя столько раз, в какую степень хотим возвести ее. По определению: A⁰ = E; A¹ = А; A^N = A^N-1·A.  Здесь E - единичная матрица.

Расчет производится по формуле:    C_ij = ∑ _ik* B_kj;  где C, A и B - матрицы. ∑ - означает суммирование по индексу k от 1 до max. При этом число строк в матрице A должно быть равно числу столбцов в матрице B. Следовательно, возводить в некоторую степень можно только квадратные матрицы, число столбцов в которых равно числу строк.

Краткое описание модуля uMatrix

В модуле определены шесть подпрограмм для непосредственного решения задачи и две вспомогательные для обеспечения правильной работы с динамической памятью.

В разделе объявлений модуля "объявляются" новые типы переменных:

1. TMatrics = array [1..max, 1..max] of Real; - двумерный массив, то есть, матрица ранга max,
2. PMatrix = ^TMatrix; - указатель на двумерный массив. Благодаря такой переменной программа будет знать о структуре массива-матрицы, размешенной в динамической памяти. Это позволит обращаться к его элементам по их индексам.

В модуле созданы следующие подпрограммы:

1. function NewMatrix: PMatrix; - размещает в динамической памяти ЭВМ массив и возвращает указатель на него. Очень простая функция. Заключается просто в вызове "стандартной" процедуры New. Легко обойтись без нее. Но с этой функцией код программы читается легче.
2. function PowMatrix(Source: PMatrix; pow: Integer): PMatrix; - самая главная. Она организует возведение матрицы, указатель на который передан ее в качестве первого параметра Source, в степень pow. Рекурсия происходит в строках

   else

   p:=MulMatrix(Source, PowMatrix(Source, Pow-1));
   По-моему, подпрограмма простая

3. procedure ShowMatrix(matrix: PMatrix; Dest: String; Size, Digs: Byte; mes: String); - выводит содержимое матрицы matrix на экран, если параметр Dest пустой, или в файл. Если файл с заданным именем существует, то информация дописывается в конец этого файла. Параметры Size, Digs предназначены для форматирования вывода. Сообщение mes можно добавить перед выводом матрицы.
4. procedure AutoFillMatrix(var Matrix: PMatrix); - подпрограмма автоматического заполнения матрицы. Она же проверяет, была ли уже создана матрица? Если нет, в этом случае указатель Matrix равен nil, то матрица инициализируется с помощью функции NewMatrix. Если захотите заполнить матрицу вводя значения элементов с клавиатуры, то лучше написать свою процедуру.
5. function MulMatrix(source1, source2: PMatrix): PMatrix; - выполняет умножение матрицы source1 на матрицу source2.
6. function UnitaryMatrix: PMatrix; - генерирует единичную матрицу. Все элементы такой матрицы равны нулю, кроме тех, что расположены на главной диагонали. Они равны единице.

При работе с динамической памятью могут возникнуть две неприятности:

• система может не выделить блок нужного Вам размера. Если Вы запускаете программу под Windows, то это скорее всего произойдет, когда попытаетесь разместить матрицу, занимающую в памяти больше 64К. Паскаль - это ДОСовская программа. Она не умеет работать с большими блоками. (Хотя справедливости ради замечу, что Borland Pascal может. Я описывал это в одном из выпусков рассылки по Турбо Паскалю)
• в предложенном простом варианте в динамической памяти могут остаться лишние массивы. Можно было бы и не делать их, точнее, своевременно удалять ненужные. Но тогда каждая из подпрограмм модуля была бы чуть сложнее. Я решил сделать так:
  • в начале работы программы я "запоминаю" вызовом самый "верхний" адрес в куче вызовом Mark(varForHeapRelease) в исполняемом блоке модуля;
  • когда программа завершает работу, нужно вызвать Release(varForHeapRelease). И я не хочу объяснять это каждому встречному и поперечному. Я использую тот факт, что у Паскаля есть процедура, которая вызывается независимо от Вашего желания. Я подменил ее своею, которая имеет такой код:
    procedure FreeAllMemory; far;
    begin
       Release(varForHeapRelease);
       ExitProc:=SaveExit;
    end;
    Она сначала освобождает память от запомненного уровня до самого верха, а потом "вспоминает" стандрантую процедуру выхода (ExitProc:=SaveExit). В переменной SaveExit для этого перед установкой своей процедуры (ExitProc:=@FreeAllMemory;) я запомнил адрес стандартной (SaveExit:=ExitProc;) в исполняемом блоке модуля.
  • КОНЕЧНО, это способ нельзя признать удачным и абсолютно надежным, но до изобретения объектов (или, по крайней мере, до их использования) приходится мириться с недостатками

Описание тестирующей программы

Программа простая. Сначала определена переменная типа PMatrix, указатель на двумерный массив. Затем вызывается процедура автоматического заполнения матрицы. Если хотите, можете написать другую процедуру, которая значение элементов будет считывать с клавиатуры. Затем все просто:

• сначала записываем содержимое исходной матрицы в файл;
• возводим матрицу во вторую степень;
• в тот же файл записываем матрицу после возведения в степень.

В результате должен получиться файл с содержимым, показанным ниже

Тестовый пример

Это результаты работы программы: содержимое файла test.txt

Исходная матрица
       -1.00        -2.00 
       -3.00        -4.00 

Матрица после возведения в степень 2
        7.00        10.00 
       15.00        22.00

Заключительное замечание (реклама объектно-ориентированного программирования)

Недостатком предложенного метода является то, что пользователь может вызвать некоторые из подпрограмм с неверными параметрами. В худшем случае ошибок времени выполнения может и не произойти. Тогда человек будет думать, что расчет выполнен правильно. Избавиться от этого недостатка можно только с помощью объектов. Там можно и удобный интерфейс предоставить и ошибки легко отслеживать.  Добавьте еще отмеченные проблемы с использованием памяти

Исходный текст программы и модуля можно скачать в виде архива по адресу http://www.borlpasc.narod.ru/Boris/Solve006.zip

Мы приглашаем Вас и Ваших друзей к сотрудничеству. Напишите, какая проблема Вас лично интересует - и мы постараемся помочь Вам. Поделитесь со всеми, если Вам удастся найти красивое решение. Присылайте свои программы, и если они хороши, то опубликуем их с обязательным указанием Вашего авторства.

По всем вопросам можно писать либо в Гостевую книгу нашего сайта на www.turbopascal.tk, либо

мне, автору этого выпуска, © Сурину Борису: surin_bp@mail.ru, ICQ: 320096696.

Постараюсь ответить на все вопросы и учесть все разумные предложения

Рассылка поддерживается сайтом www.turbopascal.tk. При перепечатке ссылка на сайт обязательна

Обращаем внимание: наш форум размещается на www.yourpascal.com !!!

Внимание: сессия и экзамены еще не начались - самое время подписаться на нашу рассылку:

Рассылки@Mail.ru Шпаргалки