Rosetta Code

Matrix-exponentiation operator: Difference between revisions

Page Discussion

Matrix-exponentiation operator (view source)

Revision as of 15:52, 14 March 2020

19 bytes added ,  4 years ago
Rename Perl 6 -> Raku, alphabetize, minor clean-up
(Merge omitted languages at bottom and add Processing)
(Rename Perl 6 -> Raku, alphabetize, minor clean-up)
Line 534:
| 0.00000 0.00000 1.00000 |
</pre>
 
=={{header|C sharp}}==
<lang csharp>using System;
using System.Collections;
using System.Collections.Generic;
using static System.Linq.Enumerable;
 
public static class MatrixExponentation
{
public static double[,] Identity(int size) {
double[,] matrix = new double[size, size];
for (int i = 0; i < size; i++) matrix[i, i] = 1;
return matrix;
}
 
public static double[,] Multiply(this double[,] left, double[,] right) {
if (left.ColumnCount() != right.RowCount()) throw new ArgumentException();
double[,] m = new double[left.RowCount(), right.ColumnCount()];
foreach (var (row, column) in from r in Range(0, m.RowCount()) from c in Range(0, m.ColumnCount()) select (r, c)) {
m[row, column] = Range(0, m.RowCount()).Sum(i => left[row, i] * right[i, column]);
}
return m;
}
 
public static double[,] Pow(this double[,] matrix, int exp) {
if (matrix.RowCount() != matrix.ColumnCount()) throw new ArgumentException("Matrix must be square.");
double[,] accumulator = Identity(matrix.RowCount());
for (int i = 0; i < exp; i++) {
accumulator = accumulator.Multiply(matrix);
}
return accumulator;
}
 
private static int RowCount(this double[,] matrix) => matrix.GetLength(0);
private static int ColumnCount(this double[,] matrix) => matrix.GetLength(1);
 
private static void Print(this double[,] m) {
foreach (var row in Rows()) {
Console.WriteLine("[ " + string.Join(" ", row) + " ]");
}
Console.WriteLine();
 
IEnumerable<IEnumerable<double>> Rows() =>
Range(0, m.RowCount()).Select(row => Range(0, m.ColumnCount()).Select(column => m[row, column]));
}
 
public static void Main() {
var matrix = new double[,] {
{ 3, 2 },
{ 2, 1 }
};
matrix.Pow(0).Print();
matrix.Pow(1).Print();
matrix.Pow(2).Print();
matrix.Pow(3).Print();
matrix.Pow(4).Print();
matrix.Pow(50).Print();
}
 
}</lang>
{{out}}
<pre style="height:30ex;overflow:scroll">
[ 1 0 ]
[ 0 1 ]
 
[ 3 2 ]
[ 2 1 ]
 
[ 13 8 ]
[ 8 5 ]
 
[ 55 34 ]
[ 34 21 ]
 
[ 233 144 ]
[ 144 89 ]
 
[ 1.61305314249046E+31 9.9692166771893E+30 ]
[ 9.9692166771893E+30 6.16131474771528E+30 ]</pre>
 
=={{header|C++}}==
Line 626 ⟶ 706:
An alternative way would be to implement <tt>operator*=</tt> and conversion from number (giving multiples of the identity matrix) for the matrix and use the generic code from [[Exponentiation operator#C++]] with support for negative exponents removed (or alternatively, implement matrix inversion as well, implement /= in terms of it, and use the generic code unchanged). Note that the algorithm used there is much faster as well.
 
=={{header|C sharpChapel}}==
<lang csharp>using System;
using System.Collections;
using System.Collections.Generic;
using static System.Linq.Enumerable;
 
This uses the '*' operator for arrays as defined in [[Matrix_multiplication#Chapel]]
public static class MatrixExponentation
<lang chapel>proc **(a, e) {
{
// create result matrix of same dimensions
public static double[,] Identity(int size) {
var r:[a.domain] a.eltType;
double[,] matrix = new double[size, size];
// and initialize to for (int i = 0; i < size; i++)identity matrix[i, i] = 1;
forall ij in r.domain do
return matrix;
r(ij) = if ij(1) == ij(2) then 1 else 0;
}
 
for 1..e do
public static double[,] Multiply(this double[,] left, double[,] right) {
r *= a;
if (left.ColumnCount() != right.RowCount()) throw new ArgumentException();
double[,] m = new double[left.RowCount(), right.ColumnCount()];
foreach (var (row, column) in from r in Range(0, m.RowCount()) from c in Range(0, m.ColumnCount()) select (r, c)) {
m[row, column] = Range(0, m.RowCount()).Sum(i => left[row, i] * right[i, column]);
}
return m;
}
 
return r;
public static double[,] Pow(this double[,] matrix, int exp) {
}</lang>
if (matrix.RowCount() != matrix.ColumnCount()) throw new ArgumentException("Matrix must be square.");
double[,] accumulator = Identity(matrix.RowCount());
for (int i = 0; i < exp; i++) {
accumulator = accumulator.Multiply(matrix);
}
return accumulator;
}
 
Usage example (like Perl):
private static int RowCount(this double[,] matrix) => matrix.GetLength(0);
<lang chapel>var m:[1..3, 1..3] int;
private static int ColumnCount(this double[,] matrix) => matrix.GetLength(1);
m(1,1) = 1; m(1,2) = 2; m(1,3) = 0;
m(2,1) = 0; m(2,2) = 3; m(2,3) = 1;
m(3,1) = 1; m(3,2) = 0; m(3,3) = 0;
 
config param n = 10;
private static void Print(this double[,] m) {
foreach (var row in Rows()) {
Console.WriteLine("[ " + string.Join(" ", row) + " ]");
}
Console.WriteLine();
 
IEnumerable<IEnumerable<double>> Rows() =>
Range(0, m.RowCount()).Select(row => Range(0, m.ColumnCount()).Select(column => m[row, column]));
}
 
public static void Main() {
var matrix = new double[,] {
{ 3, 2 },
{ 2, 1 }
};
matrix.Pow(0).Print();
matrix.Pow(1).Print();
matrix.Pow(2).Print();
matrix.Pow(3).Print();
matrix.Pow(4).Print();
matrix.Pow(50).Print();
}
 
for i in 0..n do {
writeln("Order ", i);
writeln(m ** i, "\n");
}</lang>
 
{{out}}
Order 0
<pre style="height:30ex;overflow:scroll">
[ 1 0 0 ]
[ 0 1 ]0
0 0 1
 
[ 3 2 ]
Order 1
[ 2 1 ]
1 2 0
 
0 3 1
[ 13 8 ]
1 0 0
[ 8 5 ]
 
Order 2
[ 55 34 ]
1 8 2
[ 34 21 ]
1 9 3
 
1 2 0
[ 233 144 ]
[ 144 89 ]
Order 3
 
3 26 8
[ 1.61305314249046E+31 9.9692166771893E+30 ]
4 29 9
[ 9.9692166771893E+30 6.16131474771528E+30 ]</pre>
1 8 2
Order 4
11 84 26
13 95 29
3 26 8
Order 5
37 274 84
42 311 95
11 84 26
Order 6
121 896 274
137 1017 311
37 274 84
Order 7
395 2930 896
448 3325 1017
121 896 274
Order 8
1291 9580 2930
1465 10871 3325
395 2930 896
Order 9
4221 31322 9580
4790 35543 10871
1291 9580 2930
Order 10
13801 102408 31322
15661 116209 35543
4221 31322 9580
 
=={{header|Common Lisp}}==
Line 800 ⟶ 885:
(-5315/9 66493/45 90883/135 -54445/36)
(37033/144 -27374/45 -15515/54 12109/18))
 
=={{header|Chapel}}==
 
This uses the '*' operator for arrays as defined in [[Matrix_multiplication#Chapel]]
<lang chapel>proc **(a, e) {
// create result matrix of same dimensions
var r:[a.domain] a.eltType;
// and initialize to identity matrix
forall ij in r.domain do
r(ij) = if ij(1) == ij(2) then 1 else 0;
 
for 1..e do
r *= a;
 
return r;
}</lang>
 
Usage example (like Perl):
<lang chapel>var m:[1..3, 1..3] int;
m(1,1) = 1; m(1,2) = 2; m(1,3) = 0;
m(2,1) = 0; m(2,2) = 3; m(2,3) = 1;
m(3,1) = 1; m(3,2) = 0; m(3,3) = 0;
 
config param n = 10;
 
for i in 0..n do {
writeln("Order ", i);
writeln(m ** i, "\n");
}</lang>
 
{{out}}
Order 0
1 0 0
0 1 0
0 0 1
Order 1
1 2 0
0 3 1
1 0 0
Order 2
1 8 2
1 9 3
1 2 0
Order 3
3 26 8
4 29 9
1 8 2
Order 4
11 84 26
13 95 29
3 26 8
Order 5
37 274 84
42 311 95
11 84 26
Order 6
121 896 274
137 1017 311
37 274 84
Order 7
395 2930 896
448 3325 1017
121 896 274
Order 8
1291 9580 2930
1465 10871 3325
395 2930 896
Order 9
4221 31322 9580
4790 35543 10871
1291 9580 2930
Order 10
13801 102408 31322
15661 116209 35543
4221 31322 9580
 
=={{header|D}}==
Line 1,739:
 
Matrix.Show( n )</lang>
 
=={{header|M2000 Interpreter}}==
<lang M2000 Interpreter>
Line 1,840 ⟶ 1,841:
196418 121393
</pre >
 
 
=={{header|Maple}}==
Line 2,183:
print "\n### But identity matrix can handle that\n",
$m->identity ** 1_000_000_000_000;</lang>
=={{header|Perl 6}}==
{{works with|rakudo|2015.11}}
<lang perl6>subset SqMat of Array where { .elems == all(.[]».elems) }
 
multi infix:<*>(SqMat $a, SqMat $b) {[
for ^$a -> $r {[
for ^$b[0] -> $c {
[+] ($a[$r][] Z* $b[].map: *[$c])
}
]}
]}
 
multi infix:<**> (SqMat $m, Int $n is copy where { $_ >= 0 }) {
my $tmp = $m;
my $out = [for ^$m -> $i { [ for ^$m -> $j { +($i == $j) } ] } ];
loop {
$out = $out * $tmp if $n +& 1;
last unless $n +>= 1;
$tmp = $tmp * $tmp;
}
 
$out;
}
 
multi show (SqMat $m) {
my $size = $m.flatmap( *.list».chars ).max;
say .fmt("%{$size}s", ' ') for $m.list;
}
 
my @m = [1, 2, 0],
[0, 3, 1],
[1, 0, 0];
 
for 0 .. 10 -> $order {
say "### Order $order";
show @m ** $order;
}</lang>
{{out}}
<pre>### Order 0
1 0 0
0 1 0
0 0 1
### Order 1
1 2 0
0 3 1
1 0 0
### Order 2
1 8 2
1 9 3
1 2 0
### Order 3
3 26 8
4 29 9
1 8 2
### Order 4
11 84 26
13 95 29
3 26 8
### Order 5
37 274 84
42 311 95
11 84 26
### Order 6
121 896 274
137 1017 311
37 274 84
### Order 7
395 2930 896
448 3325 1017
121 896 274
### Order 8
1291 9580 2930
1465 10871 3325
395 2930 896
### Order 9
4221 31322 9580
4790 35543 10871
1291 9580 2930
### Order 10
13801 102408 31322
15661 116209 35543
4221 31322 9580</pre>
 
=={{header|Phix}}==
Line 2,536 ⟶ 2,454:
(printf "a^~a = ~s\n" i (matrix-expt a i)))
</lang>
 
=={{header|Raku}}==
(formerly Perl 6)
{{works with|rakudo|2015.11}}
<lang perl6>subset SqMat of Array where { .elems == all(.[]».elems) }
 
multi infix:<*>(SqMat $a, SqMat $b) {[
for ^$a -> $r {[
for ^$b[0] -> $c {
[+] ($a[$r][] Z* $b[].map: *[$c])
}
]}
]}
 
multi infix:<**> (SqMat $m, Int $n is copy where { $_ >= 0 }) {
my $tmp = $m;
my $out = [for ^$m -> $i { [ for ^$m -> $j { +($i == $j) } ] } ];
loop {
$out = $out * $tmp if $n +& 1;
last unless $n +>= 1;
$tmp = $tmp * $tmp;
}
 
$out;
}
 
multi show (SqMat $m) {
my $size = $m.flatmap( *.list».chars ).max;
say .fmt("%{$size}s", ' ') for $m.list;
}
 
my @m = [1, 2, 0],
[0, 3, 1],
[1, 0, 0];
 
for 0 .. 10 -> $order {
say "### Order $order";
show @m ** $order;
}</lang>
{{out}}
<pre>### Order 0
1 0 0
0 1 0
0 0 1
### Order 1
1 2 0
0 3 1
1 0 0
### Order 2
1 8 2
1 9 3
1 2 0
### Order 3
3 26 8
4 29 9
1 8 2
### Order 4
11 84 26
13 95 29
3 26 8
### Order 5
37 274 84
42 311 95
11 84 26
### Order 6
121 896 274
137 1017 311
37 274 84
### Order 7
395 2930 896
448 3325 1017
121 896 274
### Order 8
1291 9580 2930
1465 10871 3325
395 2930 896
### Order 9
4221 31322 9580
4790 35543 10871
1291 9580 2930
### Order 10
13801 102408 31322
15661 116209 35543
4221 31322 9580</pre>
 
=={{header|Ruby}}==
Line 2,706 ⟶ 2,708:
4221 31322 9580
</pre>
 
=={{header|Scala}}==
<lang scala>class Matrix[T](matrix:Array[Array[T]])(implicit n: Numeric[T], m: ClassManifest[T])
Line 3,150 ⟶ 3,153:
<1.346269e+06,8.320400e+05>,
<8.320400e+05,5.142290e+05>>></pre>
 
=={{header|VBA}}==
No operator overloading in VBA. Implemented as a function. Can not handle scalars. Requires matrix size greater than one. Does allow for negative exponents.
Thundergnat
10,333

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/59170"