Algorytmy i Struktury Danych - Potęgowanie macierzy

Zdefiniujmy następujące operacje:

A⁰ = I (macierz jednostkowa)
A¹ = A
A² = A×A
A³ = A×A×A
A^k = A×A×…×A (k-1 mnożeń)

Pierwszy algorytm będzie algorytmem naiwnym, który wykonuje k-1 mnożeń macierzy A.

Algorytm potęgowania macierzy – wersja naiwna

Wejście:

n : stopień macierzy; n ∈ N, n > 0
A : potęgowana macierz; A ∈ R.
k : wykładnik potęgowy; k ∈ N, k ≥ 0.

Wyjście:

A : zawiera k-tą potęgę wejściowej macierzy A.

Elementy pomocnicze:

W : macierz pomocnicza, przechowuje wyniki częściowe mnożenia; W ∈ R.
P : macierz pomocnicza, przechowuje macierz A; P ∈ R.
i : do sterowania pętlą; i ∈ N.

Lista kroków:

K01: Jeśli k > 0, 
     to idź do kroku K04
K02: Ustaw macierz jednostkową w A ; A⁰ = I
K03: Zakończ
K04: Jeśli k = 1, ; A¹ = A
     to zakończ
K05: P ← A ; zapamiętaj oryginalną macierz A
K06: Dla i = 2, 3, …, k:
     wykonuj kroki K07…K08
K07:   W ← P×A ; wykonaj mnożenie
K08:   A ← W ; w A wynik mnożenia (i potęgowania)
K09: Zakończ

Przykładowe programy

Uwaga:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Program wczytuje wykładnik potęgowy k, stopień macierzy n oraz elementy macierzy w postaci n² liczb (elementy wczytywane są wierszami). Następnie oblicza k-tą potęgę macierzy i wyświetla wynik.

Dane przykładowe (przekopiuj do schowka i wklej do okna konsoli):

5
4
1 2 3 4
5 6 7 8
9 8 7 6
5 4 3 2

Pascal

// Potęgowanie macierzy
// Data: 9.02.2011
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

program mpow;

// definiujemy typy dla macierzy
// dynamicznych

type
  NReal = array of double; // typ wierszy
  MReal = array of NReal;  // typ macierzy

// procedura mnożenia macierzy
// C = A x B
//----------------------------
procedure mnoz(n : integer; A, B, C : MReal);
var
  i, j, k : integer;
  s     : double;
begin

 for i := 0 to n-1 do
    for j := 0 to n-1 do
    begin
      s := 0;
      for k := 0 to n-1 do
        s := s+A[i][k]*B[k][j];
      C[i][j] := s;
    end;
end;

// procedura przepisuje macierz A do macierzy B
//---------------------------------------------
procedure przepisz(n : integer; A, B : MReal);
var
  i, j : integer;
begin
  for i := 0 to n-1 do
    for j := 0 to n-1 do
      B[i][j] := A[i][j];
end;

// procedura ustawia w macierzy A macierz jednostkową
//---------------------------------------------------
procedure jednostkowa(n : integer; A : MReal);
var
  i, j : integer;
begin
  for i := 0 to n-1 do
  begin
    for j := 0 to n-1 do
      A[i][j] := 0;
    A[i][i] := 1;
  end;
end;

// procedura wylicza potęgę k-tą macierzy A
//-----------------------------------------
procedure potega(k, n : integer; A : MReal);
var
  P, W : MReal;
  i   : integer;
begin
  if k = 0 then jednostkowa(n, A)
  else if k > 1 then
  begin

    // tworzymy macierze pomocnicze P i W
    SetLength(P, n);
    SetLength(W, n);
    for i := 0 to n-1 do
    begin
      SetLength(P[i], n);
      SetLength(W[i], n);
    end;

    // macierz A zapamiętujemy w P
    przepisz(n, A, P);

    // w pętli wykonujemy kolejne mnożenia, 
    // wynik zawsze w A
    for i := 2 to k do
    begin
      mnoz(n, P, A, W);   // W <- P x A
      przepisz(n, W, A); // A <- P x A
    end;

    // usuwamy macierze P i W
    for i := 0 to n-1 do
    begin
      SetLength(P[i], 0);
      SetLength(W[i], 0);
    end;
    SetLength(P, 0);
    SetLength(W, 0);
  end;
end;

//*** PROGRAM GŁÓWNY ***
//----------------------

var
  A : MReal;
  n, i, j, k : integer;

begin

  // wczytujemy wykładnik k
  // oraz stopień macierzy n
  read(k, n);

  // tworzymy macierz dynamiczną
  // i wczytujemy dane wierszami
  SetLength(A, n);
  for i := 0 to n-1 do
  begin
    SetLength(A[i], n);
    for j := 0 to n-1 do
      read(A[i][j]);
  end;

  // obliczamy k-tą potęgę A
  potega(k, n, A);

  // wyświetlamy wyniki
  writeln;

  for i := 0 to n-1 do
  begin
    for j := 0 to n-1 do
      write(A[i][j]:10:2, ' ');
    writeln;
  end;

  // usuwamy macierz A
  for i := 0 to n-1 do
    SetLength(A[i], 0);
  SetLength(A, 0);

end.

C++

// Potęgowanie macierzy
// Data: 9.02.2011
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

#include <iostream>
#include <iomanip>

using namespace std;

// procedura mnożenia macierzy
// C = A x B
//----------------------------
void mnoz(int n, 
          double ** A, 
          double ** B, 
          double ** C)
{
  int i, j, k;
  double s;

  for(i = 0; i < n; i++)
    for(j = 0; j < n; j++)
    {
      s = 0;
      for(k = 0; k < n; k++)
        s += A[i][k]*B[k][j];
      C[i][j] = s;
    }
}

// procedura przepisuje macierz A do macierzy B
//---------------------------------------------
void przepisz(int n, 
              double ** A, 
              double ** B)
{
  int i, j;

  for(i = 0; i < n; i++)
    for(j = 0; j < n; j++)
      B[i][j] = A[i][j];
}

// procedura ustawia w macierzy A macierz jednostkową
//---------------------------------------------------
void jednostkowa(int n, 
                 double ** A)
{
  int i, j;

  for(i = 0; i < n; i++)
  {
    for(j = 0; j < n; j++)
      A[i][j] = 0;
    A[i][i] = 1;
  }
}

// procedura wylicza potęgę k-tą macierzy A
//-----------------------------------------
void potega(int k, 
            int n, 
            double ** A)
{
  double ** P, ** W;
  int i;

  if(!k) jednostkowa(n, A);
  else if(k > 1)
  {
    // tworzymy macierze pomocnicze P i W
    P = new double * [n];
    W = new double * [n];
    for(i = 0; i < n; i++)
    {
      P[i] = new double[n];
      W[i] = new double[n];
    }

    // macierz A zapamiętujemy w P
    przepisz(n, A, P);

    // w pętli wykonujemy kolejne mnożenia, 
    // wynik zawsze w A
    for(i = 2; i <= k; i++)
    {
      mnoz(n, P, A, W);   // W <- P x A
      przepisz(n, W, A); // A <- P x A
    }

    // usuwamy macierze P i W
    for(i = 0; i < n; i++)
    {
      delete [] P[i];
      delete [] W[i];
    }
    delete [] P;
    delete [] W;
  }
}

//*** PROGRAM GŁÓWNY ***
//----------------------

int main()
{
  double ** A;
  int n, i, j, k;

  cout << fixed << setprecision(2);

  // wczytujemy wykładnik k oraz stopień macierzy n
  cin >> k >> n;

  // tworzymy macierz dynamiczną
  // i wczytujemy dane wierszami
  A = new double * [n];
  for(i = 0; i < n; i++)
  {
    A[i] = new double[n];
    for(j = 0; j < n; j++)
      cin >> A[i][j];
  }

  // obliczamy k-tą potęgę A
  potega(k, n, A);

  // wyświetlamy wyniki
  cout << endl;

  for(i = 0; i < n; i++)
  {
    for(j = 0; j < n; j++)
      cout << setw(10) << A[i][j] << " ";
    cout << endl;
  }

  // usuwamy macierz A
  for(i = 0; i < n; i++)
    delete [] A[i];
  delete [] A;

  return 0;
}

Basic

' Potęgowanie macierzy
' Data: 9.02.2011
' (C)2020 mgr Jerzy Wałaszek
'---------------------------

' procedura mnożenia macierzy
' C = A x B
'----------------------------
Sub mnoz(n As Integer, _
         A As Double Ptr Ptr, _
         B As Double Ptr Ptr, _
         C As Double Ptr Ptr)

  Dim As Integer i, j, k
  Dim As Double s

  For i = 0 To n-1
    For j = 0 To n-1
      s = 0
      For k = 0 To n-1
        s += A[i][k]*B[k][j]
      Next
      C[i][j] = s
    Next
  Next
End Sub

' procedura przepisuje macierz A
' do macierzy B
'-------------------------------
Sub przepisz(n As Integer, _
             A As Double Ptr Ptr, _
             B As Double Ptr Ptr)

  Dim As Integer i, j

  For i = 0 To n-1
    For j = 0 To n-1
      B[i][j] = A[i][j]
    Next
  Next
End Sub

' procedura ustawia w macierzy A
' macierz jednostkową
'-------------------------------
Sub jednostkowa(n As Integer, _
                A As Double Ptr Ptr)

  Dim As Integer i, j

  For i = 0 To n-1
    For j = 0 To n-1
      A[i][j] = 0
    Next
    A[i][i] = 1
  Next
End Sub

' procedura wylicza potęgę k-tą
' macierzy A
'------------------------------
Sub potega(k As Integer, _
           n As Integer, _
           A As Double Ptr Ptr)

  Dim As Double Ptr Ptr P, W
  Dim As Integer i

  If k = 0 Then
    jednostkowa(n, A)
  ElseIf k > 1 Then

    ' tworzymy macierze pomocnicze P i W
    P = New Double Ptr[n]
    W = New Double Ptr[n]
    For i = 0 To n-1
      P[i] = New Double[n]
      W[i] = New Double[n]
    Next

    ' macierz A zapamiętujemy w P
    przepisz(n, A, P)

    ' w pętli wykonujemy kolejne mnożenia, 
    ' wynik zawsze w A
    For i = 2 To k
      mnoz(n, P, A, W)   ' W <- P x A
      przepisz(n, W, A) ' A <- P x A
    Next

    ' usuwamy macierze P i W
    For i = 0 To n-1
      Delete [] P[i]
      Delete W[i]
    Next
    Delete [] P
    Delete [] W

  End If
End Sub

Dim As Double Ptr Ptr A
Dim As Integer n, i, j, k

Open Cons For Input As #1

' wczytujemy wykładnik k
' oraz stopień macierzy n
Input #1, k, n

' tworzymy macierz dynamiczną
' i wczytujemy dane wierszami

A = New Double Ptr[n]
For i = 0 To n-1
  A[i] = New Double[n]
  For j = 0 To n-1
    Input #1, A[i][j]
  Next
Next

Close #1

' obliczamy k-tą potęgę A
potega(k, n, A)

' wyświetlamy wyniki
For i = 0 To n-1
  For j = 0 To n-1
    Print Using "#######.## ";A[i][j];
  Next
  Print
Next

' usuwamy macierz A
For i = 0 To n-1
  Delete [] A[i]
Next
Delete [] A

End

Python (dodatek)

# Potęgowanie macierzy
# Data: 12.04.2024
# (c)2024 mgr Jerzy Wałaszek
#---------------------------


# procedura mnożenia macierzy
# c = a x b
#----------------------------
def mnoz(a, b, c):
    n = len(a)
    for i in range(n):
        for j in range(n):
            s = 0.0
            for k in range(n):
                s += a[i][k]*b[k][j]
            c[i][j] = s


# procedura przepisuje macierz a
# do macierzy b
#-------------------------------
def przepisz(a, b):
    n = len(a)
    for i in range(n):
        b[i] = a[i].copy()


# procedura ustawia w macierzy a
# macierz jednostkową
#-------------------------------
def jednostkowa(a):
    n = len(a)
    for i in range(n):
        for j in range(n):
            a[i][j] = 0.0
        a[i][i] = 1.0


# procedura wylicza k-tą
# potęgę macierzy a
#-----------------------
def potega(k, a):
    if k == 0: jednostkowa(a)
    elif k > 1:
        # macierze P i W
        n = len(a)
        arr = [0] * n
        p = [arr.copy()] * n
        w = p.copy()
        # macierz a zapamiętujemy w P
        przepisz(a, p)
        # w pętli wykonujemy
        # kolejne mnożenia, 
        # wynik zawsze w a
        for i in range(k-1):
            mnoz(p, a, w)  # W <- P x a
            przepisz(w, a) # a <- P x a


# *** PROGRAM GŁÓWNY ***
# ----------------------

# wczytujemy wykładnik k
# oraz stopień macierzy n
k = int(input())
n = int(input())

# tworzymy macierz dynamiczną
# i wczytujemy dane wierszami
a = []
for i in range(n):
    arr = input().split()
    arr = [float(arr[j]) for j in range(n)]
    a.append(arr)
# obliczamy k-tą potęgę a
potega(k, a)
# wyświetlamy wyniki
print()
for i in range(n):
    for j in range(n):
        print("%10.2f " % a[i][j], end="")
    print()
# usuwamy macierz a
a = None
arr = None

Wynik:

5
4
1 2 3 4
5 6 7 8
9 8 7 6
5 4 3 2

 412928.00  413184.00  413440.00  413696.00 
1052928.00 1053184.00 1053440.00 1053696.00 
1187072.00 1186816.00 1186560.00 1186304.00 
 547072.00  546816.00  546560.00  546304.00

Na początek: podrozdziału strony

Zauważmy następującą własność mnożenia macierzy:

A^n+m = Aⁿ×A^m

Dalej:

A² = A×A
A⁴ = A²×A²A⁸ = A⁴×A⁴…

Jak widzimy, potęgi o wykładnikach będących potęgami liczby 2 możemy łatwo otrzymywać mnożąc macierz przez siebie. Idea nowego algorytmu potęgowania przedstawię na prostym przykładzie. Załóżmy, iż chcemy obliczyć A¹³. Zatem k = 13. Możemy zapisać:

A¹³ = A⁸⁺⁴⁺¹ = A⁸×A⁴×A

Zwróć uwagę, iż mnożymy przez siebie macierze, które są potęgami pierwotnej macierzy A o wykładnikach równych 8, 4 i 1. Na takie liczby rozkłada się liczba 13. Potęgi te bardzo łatwo otrzymać.

Potrzebujemy 3 mnożeń do wyliczenia wszystkich niezbędnych potęg macierzy A. Aby ostatecznie otrzymać A¹³ potrzebne są jeszcze dwa dodatkowe mnożenia, czyli w sumie 5 mnożeń. Poprzedni algorytm doszedłby do tego samego wyniku po wykonaniu aż 12 mnożeń. Przy większych potęgach zysk jest jeszcze większy.

Wykładnika k wcale nie musimy rozbijać na potęgi liczby 2. W pamięci komputera k jest przechowywane w postaci binarnej. Wystarczy testować stan odpowiednich bitów k:

kolejne potęgi liczby 2 =	2⁴	2³	2²	2¹	2⁰
bity k =	0	1	1	0	1
numery pozycji bitów kolejnych =	4	3	2	1	0

Zasada pracy algorytmu jest następująca:

W macierzy wyniku W ustawiamy macierz jednostkową. Dopóki k jest większe od zera, testujemy najmłodszy bit k i jeśli ma on wartość 1, to wykonujemy mnożenie W = W×A. Bity k przesuwamy o 1 w prawo, aby pozbyć się testowanego bitu. Jeśli po tej operacji k ma wartość 0, to potęgowanie jest zakończone i przerywamy pętlę. W przeciwnym razie macierz A podnosimy do kwadratu: A = A×A i pętlę kontynuujemy. Gdy pętla się zakończy, przepisujemy macierz W do A.

Algorytm potęgowania macierzy – wersja ulepszona

Wejście:

n : stopień macierzy; n ∈ N, n > 0.
A : potęgowana macierz; A ∈ R.
k : wykładnik potęgowy; k ∈ N, k ≥ 0.

Wyjście:

A : zawiera k-tą potęgę wejściowej macierzy A.

Elementy pomocnicze:

W : macierz pomocnicza, tworzony jest w niej wynik potęgowania; W ∈ R.
P : macierz pomocnicza; P ∈ R.

Lista kroków:

K01: Ustaw macierz jednostkową w W
K02: Dopóki k > 0, ; w pętli obliczamy k-tą potęgę A
     wykonuj kroki K03…K09
K03:   Jeśli (k and 1) = 0, ; testujemy kolejne bity k
       to idź do kroku K06
K04    P ← W×A ; jeśli bit jest ustawiony, to do W dołączamy A
K05:   W ← P
K06:   k ← k shr 1 ; przesuwamy w prawo bity k
K07:   Jeśli k = 0, 
       to idź do kroku K10
K08:   P ← A×A ; wyznaczamy kolejny kwadrat A
K09    A ← P
K10: A ← W ; wynik do A
K11: Zakończ

Przykładowe programy

Uwaga:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Dane przykładowe (przekopiuj do schowka i wklej do okna konsoli):

5
4
1 2 3 4
5 6 7 8
9 8 7 6
5 4 3 2

Pascal

// Szybkie potęgowanie macierzy
// Data: 19.02.2012
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

program mpow;

// definiujemy typy
// dla macierzy dynamicznych
type
  NReal = array of double; // typ wierszy
  MReal = array of NReal;  // typ macierzy

// procedura mnożenia macierzy
// C = A x B
//----------------------------
procedure mnoz(n : integer;
               A, B, C : MReal);
var
  i, j, k : integer;
  s     : double;
begin

 for i := 0 to n-1 do
    for j := 0 to n-1 do
    begin
      s := 0;
      for k := 0 to n-1 do
        s := s+A[i][k]*B[k][j];
      C[i][j] := s;
    end;
end;

// procedura przepisuje macierz A
// do macierzy B
//-------------------------------
procedure przepisz(n : integer;
                   A, B : MReal);
var
  i, j : integer;
begin
  for i := 0 to n-1 do
    for j := 0 to n-1 do
      B[i][j] := A[i][j];
end;

// procedura ustawia w macierzy A
// macierz jednostkową
//-------------------------------
procedure jednostkowa(n : integer;
                      A : MReal);
var
  i, j : integer;
begin
  for i := 0 to n-1 do
  begin
    for j := 0 to n-1 do
      A[i][j] := 0;
    A[i][i] := 1;
  end;
end;

// procedura wylicza potęgę k-tą
// macierzy A
//------------------------------
procedure potega(k, n : integer;
                 A : MReal);
var
  P, W : MReal;
  i   : integer;
begin
  // tworzymy macierze W i P

  SetLength(W, n);
  SetLength(P, n);

  for i := 0 to n-1 do
  begin
    SetLength(W[i], n);
    SetLength(P[i], n);
  end;

  // w macierzy W ustawiamy
  // macierz jednostkową
  jednostkowa(n, W);

  // w pętli obliczamy potęg
  // macierzy A w W
  while k > 0 do
  begin
    if k and 1 = 1 then // testujemy
    begin               // najmłodszy bit k
      mnoz(n, W, A, P); // jeśli ustawiony, to
      przepisz(n, P, W); // W = W x A
    end;
    k := k shr 1; // przesuwamy bity k w prawo
    if k = 0 then break; // jeśli brak bitów 1, 
                         // przerywamy
    mnoz(n, A, A, P); // A = A x A
    przepisz(n, P, A);
  end;

  // wynik potęgowania wraca do macierzy A
  przepisz(n, W, A);

  // usuwamy macierze W i P
  for i := 0 to n-1 do
  begin
    SetLength(W[i], 0);
    SetLength(P[i], 0);
  end;
  SetLength(W, 0);
  SetLength(P, 0);
end;

//*** PROGRAM GŁÓWNY ***
//----------------------

var
  A : MReal;
  n, i, j, k : integer;

begin

  // wczytujemy wykładnik k
  // oraz stopień macierzy n
  read(k, n);

  // tworzymy macierz dynamiczną
  // i wczytujemy dane wierszami
  SetLength(A, n);
  for i := 0 to n-1 do
  begin
    SetLength(A[i], n);
    for j := 0 to n-1 do
      read(A[i][j]);
  end;

  // obliczamy k-tą potęgę A
  potega(k, n, A);

  // wyświetlamy wyniki
  writeln;
  for i := 0 to n-1 do
  begin
    for j := 0 to n-1 do
      write(A[i][j]:10:2, ' ');
    writeln;
  end;

  // usuwamy macierz A
  for i := 0 to n-1 do
    SetLength(A[i], 0);
  SetLength(A, 0);

end.

C++

// Szybkie potęgowanie macierzy
// Data: 19.02.2012
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

#include <iostream>
#include <iomanip>

using namespace std;

// procedura mnożenia macierzy
// C = A x B
//----------------------------
void mnoz(int n, 
          double ** A, 
          double ** B, 
          double ** C)
{
  int i, j, k;
  double s;

  for(i = 0; i < n; i++)
    for(j = 0; j < n; j++)
    {
      s = 0;
      for(k = 0; k < n; k++)
        s += A[i][k]*B[k][j];
      C[i][j] = s;
    }
}

// procedura przepisuje
// macierz A do macierzy B
//------------------------
void przepisz(int n, 
              double ** A, 
              double ** B)
{
  int i, j;

  for(i = 0; i < n; i++)
    for(j = 0; j < n; j++)
      B[i][j] = A[i][j];
}

// procedura ustawia w macierzy A
// macierz jednostkową
//-------------------------------
void jednostkowa(int n, 
                 double ** A)
{
  int i, j;

  for(i = 0; i < n; i++)
  {
    for(j = 0; j < n; j++)
      A[i][j] = 0;
    A[i][i] = 1;
  }
}

// procedura wylicza potęgę
// k-tą macierzy A
//-------------------------
void potega(int k, 
            int n, 
            double ** A)
{
  double ** P, ** W;
  int i;

  // tworzymy macierze W i P
  W = new double * [n];
  P = new double * [n];
  for(i = 0; i < n; i++)
  {
    W[i] = new double[n];
    P[i] = new double[n];
  }

  // w macierzy W ustawiamy
  // macierz jednostkową
  jednostkowa(n, W);

  // w pętli obliczamy potęgę
  // macierzy A w W
  while(k)
  {
    if(k&1) // testujemy najmłodszy bit k
    {
      mnoz(n, W, A, P); // jeśli ustawiony, 
      przepisz(n, P, W); // to W = W x A
    }
    k >>= 1; // przesuwamy bity k w prawo
    if(!k) break; // jeśli brak bitów 1, 
                  // przerywamy
    mnoz(n, A, A, P); // A = A x A
    przepisz(n, P, A);
  }

  // wynik potęgowania wraca
  // do macierzy A
  przepisz(n, W, A);

  // usuwamy macierze W i P
  for(i = 0; i < n; i++)
  {
    delete [] W[i];
    delete [] P[i];
  }
  delete [] W;
  delete [] P;
}

//*** PROGRAM GŁÓWNY ***
//----------------------
int main()
{
  double ** A;
  int n, i, j, k;

  cout << fixed << setprecision(2);

  // wczytujemy wykładnik k
  // oraz stopień macierzy n
  cin >> k >> n;

  // tworzymy macierz dynamiczną
  // i wczytujemy dane wierszami
  A = new double * [n];
  for(i = 0; i < n; i++)
  {
    A[i] = new double[n];
    for(j = 0; j < n; j++)
      cin >> A[i][j];
  }

  // obliczamy k-tą potęgę A

  potega(k, n, A);

  // wyświetlamy wyniki

  cout << endl;

  for(i = 0; i < n; i++)
  {
    for(j = 0; j < n; j++)
      cout << setw(10)
           << A[i][j] << " ";
    cout << endl;
  }

  // usuwamy macierz A

  for(i = 0; i < n; i++)
    delete [] A[i];
  delete [] A;

  return 0;
}

Basic

' Szybkie potęgowanie macierzy
' Data: 9.02.2011
' (C)2020 mgr Jerzy Wałaszek
'-----------------------------

' procedura mnożenia macierzy
' C = A x B
'----------------------------
Sub mnoz(n As Integer, _
         A As Double Ptr Ptr, _
         B As Double Ptr Ptr, _
         C As Double Ptr Ptr)

  Dim As Integer i, j, k
  Dim As Double s

  For i = 0 To n-1
    For j = 0 To n-1
      s = 0
      For k = 0 To n-1
        s += A[i][k]*B[k][j]
      Next
      C[i][j] = s
    Next
  Next
End Sub

' procedura przepisuje macierz A do macierzy B
'---------------------------------------------
Sub przepisz(n As Integer, _
             A As Double Ptr Ptr, _
             B As Double Ptr Ptr)

  Dim As Integer i, j

  For i = 0 To n-1
    For j = 0 To n-1
      B[i][j] = A[i][j]
    Next
  Next
End Sub

' procedura ustawia w macierzy A
' macierz jednostkową
'-------------------------------
Sub jednostkowa(n As Integer, _
                A As Double Ptr Ptr)

  Dim As Integer i, j

  For i = 0 To n-1
    For j = 0 To n-1
      A[i][j] = 0
    Next
    A[i][i] = 1
  Next
End Sub

' procedura wylicza potęgę k-tą macierzy A
'-----------------------------------------
Sub potega(k As Integer, _
           n As Integer, _
           A As Double Ptr Ptr)

  Dim As Double Ptr Ptr P, W
  Dim i As Integer

  ' tworzymy macierze W i P
  W = New Double Ptr[n]
  P = New Double Ptr[n]
  For i = 0 To n-1
    W[i] = New Double[n]
    P[i] = New Double[n]
  Next

  ' w macierzy W ustawiamy macierz jednostkową
  jednostkowa(n, W)

  ' w pętli obliczamy potęgę macierzy A w W
  While k > 0
    if(k And 1) = 1 Then ' najmłodszy bit
      mnoz(n, W, A, P)   ' jeśli ustawiony, to
      przepisz(n, P, W) ' W = W x A
    End If
    k Shr= 1          ' przesuwamy bity k w prawo
    If k = 0 Then Exit While ' jeśli brak bitów 1, 
                      ' przerywamy
    mnoz(n, A, A, P)     ' A = A x A
    przepisz(n, P, A)
  Wend

  ' wynik potęgowania wraca do macierzy A
  przepisz(n, W, A)

  ' usuwamy macierze W i P
  For i = 0 To n-1
    Delete [] W[i]
    Delete [] P[i]
  Next
  Delete [] W
  Delete [] P
 
End Sub

Dim As Double Ptr Ptr A
Dim As Integer n, i, j, k

Open Cons For Input As #1

' wczytujemy wykładnik k oraz stopień macierzy n
Input #1, k, n

' tworzymy macierz dynamiczną
' i wczytujemy dane wierszami
A = New Double Ptr [n]
For i = 0 To n-1
  A[i] = New Double[n]
  For j = 0 To n-1
    Input #1, A[i][j]
  Next
Next

Close #1

' obliczamy k-tą potęgę A
potega(k, n, A)

' wyświetlamy wyniki
For i = 0 To n-1
  For j = 0 To n-1
    Print Using "#######.## ";A[i][j];: Next
  Print
Next

' usuwamy macierz A
For i = 0 To n-1
  Delete [] A[i]
Next
Delete [] A

End

Python (dodatek)

# Szybkie potęgowanie macierzy
# Data: 13.04.2024
# (C)2024 mgr Jerzy Wałaszek
#-----------------------------


# procedura mnożenia macierzy
# Z = X x Y
#----------------------------
def mnoz(x, y, z):
    n = len(x)
    for i in range(n):
        for j in range(n):
            s = 0.0
            for k in range(n):
                s += x[i][k]*y[k][j]
            z[i][j] = s


# procedura przepisuje
# macierz X do macierzy Y
#------------------------
def przepisz(x, y):
    n = len(x)
    for i in range(n):
        y[i] = x[i].copy()


# procedura ustawia w macierzy X
# macierz jednostkową
#-------------------------------
def jednostkowa(x):
    n = len(x)
    for i in range(n):
        for j in range(n):
            x[i][j] = 0.0
        x[i][i] = 1.0


# procedura wylicza
# k-tą potęgę macierzy X
#-----------------------
def potega(k, x):

    # tworzymy macierze W i P
    n = len(x)
    w = []
    p = []
    for i in range(n):
        a = []
        for j in range(n):
            a.append(0.0)
        w.append(a.copy())
        p.append(a.copy())
    # w macierzy W ustawiamy
    # macierz jednostkową
    jednostkowa(w)
    # w pętli obliczamy potęgę
    # macierzy X w W
    while k:
        if k&1: # testujemy najmłodszy bit k
            mnoz(w, x, p) # jeśli ustawiony, to
            przepisz(p, w) # W = W x X
        k >>= 1 # przesuwamy bity k w prawo
        if not k: break # jeśli brak bitów 1,
                        # przerywamy
        mnoz(x, x, p) # X = X x X
        przepisz(p, x)
    # wynik potęgowania wraca do macierzy X
    przepisz(w, x)


#*** PROGRAM GŁÓWNY ***
#----------------------

# wczytujemy wykładnik k
# oraz stopień macierzy n
k = int(input())
n = int(input())
# tworzymy macierz dynamiczną
# i wczytujemy dane wierszami
a = []
for i in range(n):
    s = input().split()
    arr = [0.0] * n
    for j in range(n):
        arr[j] = float(s[j])
    a.append(arr)
# obliczamy k-tą potęgę A
potega(k, a)
# wyświetlamy wyniki
print()
for i in range(n):
    for j in range(n):
        print("%10.2f " % a[i][j], end="")
    print()

Na początek: podrozdziału strony

Materiały tylko do użytku dydaktycznego. Ich kopiowanie i powielanie jest dozwolone
pod warunkiem podania źródła oraz niepobierania za to pieniędzy.

Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl

Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.

Informacje dodatkowe.

Potęgowanie macierzy

Problem

Rozwiązanie nr 1

Algorytm potęgowania macierzy – wersja naiwna

Wejście:

Wyjście:

Elementy pomocnicze:

Lista kroków:

Przykładowe programy

Rozwiązanie nr 2

Algorytm potęgowania macierzy – wersja ulepszona

Wejście:

Wyjście:

Elementy pomocnicze:

Lista kroków:

Przykładowe programy