Matrix-exponentiation operator: Difference between revisions
m
→{{header|Wren}}: Minor tidy
m (→{{header|Wren}}: Minor tidy) |
|||
| (4 intermediate revisions by 2 users not shown) | |||
Line 226:
begin
for K in RL'Range (1) loop
Sum := Sum + X (
end loop;
R (I, J) := Sum;
Line 237:
begin
for I in A'Range (1) loop
for J in A'Range (
Put (A (I, J));
end loop;
Line 251:
end Test_Matrix;</syntaxhighlight>
This solution is not tested, because the available version of GNAT GPL Ada compiler (20070405-41) does not provide an implementation of the standard library.
''(Another person is talking here:)'' I have made small corrections and tested this in 2023, and it did not work as I expected. However, I have questions about the mathematical libraries. I tried both GCC 12 and GCC 13. (I also tried the last GNAT Community Edition, but it no longer functions on my system.) What might be needed here is one's own eigensystem routine.
On the other hand, I did get a version working to raise a real matrix to a natural number power, thus demonstrating the correctness of the approach:
<syntaxhighlight lang="ada">
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Float_Text_IO; use Ada.Float_Text_IO;
with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;
procedure Test_Matrix is
procedure Put (A : Real_Matrix) is
begin
for I in A'Range (1) loop
for J in A'Range (2) loop
Put (" ");
Put (A (I, J));
end loop;
New_Line;
end loop;
end Put;
function "**" (A : Real_Matrix; Power : Integer) return Real_Matrix is
L : Real_Vector (A'Range (1));
X : Real_Matrix (A'Range (1), A'Range (2));
R : Real_Matrix (A'Range (1), A'Range (2));
RL : Real_Vector (A'Range (1));
begin
Eigensystem (A, L, X);
for I in L'Range loop
RL (I) := L (I) ** Power;
end loop;
for I in R'Range (1) loop
for J in R'Range (2) loop
declare
Sum : Float := 0.0;
begin
for K in RL'Range loop
Sum := Sum + X (I, K) * RL (K) * X (J, K);
end loop;
R (I, J) := Sum;
end;
end loop;
end loop;
return R;
end "**";
M : Real_Matrix (1..2, 1..2) := ((3.0, 2.0), (2.0, 1.0));
begin
Put_Line ("M ="); Put (M);
Put_Line ("M**0 ="); Put (M**0);
Put_Line ("M**1 ="); Put (M**1);
Put_Line ("M**2 ="); Put (M**2);
Put_Line ("M**3 ="); Put (M**3);
Put_Line ("M**50 ="); Put (M**50);
end Test_Matrix;
</syntaxhighlight>
{{out}}
<pre>
M =
3.00000E+00 2.00000E+00
2.00000E+00 1.00000E+00
M**0 =
1.00000E+00 0.00000E+00
0.00000E+00 1.00000E+00
M**1 =
3.00000E+00 2.00000E+00
2.00000E+00 1.00000E+00
M**2 =
1.30000E+01 8.00000E+00
8.00000E+00 5.00000E+00
M**3 =
5.50000E+01 3.40000E+01
3.40000E+01 2.10000E+01
M**50 =
1.61305E+31 9.96919E+30
9.96919E+30 6.16130E+30
</pre>
=={{header|ALGOL 68}}==
Line 1,868 ⟶ 1,943:
<syntaxhighlight lang="j">pow=: pow1=: 4 : 'mp/ mp~^:(I.|.#:y) x'</syntaxhighlight>
This implements an optimization where the exponent is represented in base 2, and repeated squaring is used to create a list of relevant powers of the base matrix, which are then combined using matrix multiplication. Note, however, that these two definitions treat a zero exponent differently (<code>m pow0 0</code> gives an identity matrix whose shape matches m, while <code>m pow1 0</code> gives a scalar 1).
Example use:
<syntaxhighlight lang=J>
(3 2,:2 1) pow 3
55 34
34 21
</syntaxhighlight>
=={{header|JavaScript}}==
Line 3,980 ⟶ 4,057:
The Matrix class in the above module also has a 'pow' method but, as an alternative, overloads the otherwise unused '^' operator to provide the same functionality.
<syntaxhighlight lang="
import "./fmt" for Fmt
var m = Matrix.new([[0, 1], [1, 1]])
| |||