Gauss-Jordan matrix inversion
Last updated
Was this helpful?
Last updated
Was this helpful?
Uses bits and pieces from other tasks, primarily.
sub gauss-jordan-invert (@m where &is-square) {
^@m .map: { @m[$_].append: identity(+@m)[$_] };
@m.&rref[*]»[+@m .. *];
}
sub is-square (@m) { so @m == all @m[*] }
sub identity ($n) { [ 1, |(0 xx $n-1) ], *.rotate(-1).Array ... *.tail }
# reduced row echelon form (from 'Gauss-Jordan elimination' task)
sub rref (@m) {
my ($lead, $rows, $cols) = 0, @m, @m[0];
for ^$rows -> $r {
$lead < $cols or return @m;
my $i = $r;
until @m[$i;$lead] {
++$i == $rows or next;
$i = $r;
++$lead == $cols and return @m;
}
@m[$i, $r] = @m[$r, $i] if $r != $i;
@m[$r] »/=» $ = @m[$r;$lead];
for ^$rows -> $n {
next if $n == $r;
@m[$n] »-=» @m[$r] »×» (@m[$n;$lead] // 0);
}
++$lead;
}
@m
}
sub rat-or-int ($num) {
return $num unless $num ~~ Rat|FatRat;
return $num.narrow if $num.narrow ~~ Int;
$num.nude.join: '/';
}
sub say_it ($message, @array) {
my $max;
@array.map: {$max max= max $_».&rat-or-int.comb(/\S+/)».chars};
say "\n$message";
$_».&rat-or-int.fmt(" %{$max}s").put for @array;
}
multi to-matrix ($str) { [$str.split(';').map(*.words.Array)] }
multi to-matrix (@array) { @array }
sub hilbert-matrix ($h) {
[ (1..$h).map( -> $n { [ ($n ..^ $n + $h).map: { (1/$_).FatRat } ] } ) ]
}
my @tests =
'1 2 3; 4 1 6; 7 8 9',
'2 -1 0; -1 2 -1; 0 -1 2',
'-1 -2 3 2; -4 -1 6 2; 7 -8 9 1; 1 -2 1 3',
'1 2 3 4; 5 6 7 8; 9 33 11 12; 13 14 15 17',
'3 1 8 9 6; 6 2 8 10 1; 5 7 2 10 3; 3 2 7 7 9; 3 5 6 1 1',
# Test with a Hilbert matrix
hilbert-matrix 10;
@tests.map: {
my @matrix = .&to-matrix;
say_it( " {'=' x 20} Original Matrix: {'=' x 20}", @matrix );
say_it( ' Gauss-Jordan Inverted:', my @invert = gauss-jordan-invert @matrix );
say_it( ' Re-inverted:', gauss-jordan-invert @invert».Array );
} ==================== Original Matrix: ====================
1 2 3
4 1 6
7 8 9
Gauss-Jordan Inverted:
-13/16 1/8 3/16
1/8 -1/4 1/8
25/48 1/8 -7/48
Re-inverted:
1 2 3
4 1 6
7 8 9
==================== Original Matrix: ====================
2 -1 0
-1 2 -1
0 -1 2
Gauss-Jordan Inverted:
3/4 1/2 1/4
1/2 1 1/2
1/4 1/2 3/4
Re-inverted:
2 -1 0
-1 2 -1
0 -1 2
==================== Original Matrix: ====================
-1 -2 3 2
-4 -1 6 2
7 -8 9 1
1 -2 1 3
Gauss-Jordan Inverted:
-21/23 17/69 13/138 19/46
-38/23 15/23 1/23 15/23
-16/23 25/69 11/138 9/46
-13/23 16/69 -2/69 13/23
Re-inverted:
-1 -2 3 2
-4 -1 6 2
7 -8 9 1
1 -2 1 3
==================== Original Matrix: ====================
1 2 3 4
5 6 7 8
9 33 11 12
13 14 15 17
Gauss-Jordan Inverted:
19/184 -199/184 -1/46 1/2
1/23 -2/23 1/23 0
-441/184 813/184 -1/46 -3/2
2 -3 0 1
Re-inverted:
1 2 3 4
5 6 7 8
9 33 11 12
13 14 15 17
==================== Original Matrix: ====================
3 1 8 9 6
6 2 8 10 1
5 7 2 10 3
3 2 7 7 9
3 5 6 1 1
Gauss-Jordan Inverted:
-4525/6238 2529/6238 -233/3119 1481/3119 -639/6238
1033/6238 -1075/6238 342/3119 -447/3119 871/6238
1299/6238 -289/6238 -204/3119 -390/3119 739/6238
782/3119 -222/3119 237/3119 -556/3119 -177/3119
-474/3119 -17/3119 -24/3119 688/3119 -140/3119
Re-inverted:
3 1 8 9 6
6 2 8 10 1
5 7 2 10 3
3 2 7 7 9
3 5 6 1 1
==================== Original Matrix: ====================
1 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10
1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11
1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12
1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13
1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14
1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15
1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16
1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17
1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18
1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19
Gauss-Jordan Inverted:
100 -4950 79200 -600600 2522520 -6306300 9609600 -8751600 4375800 -923780
-4950 326700 -5880600 47567520 -208107900 535134600 -832431600 770140800 -389883780 83140200
79200 -5880600 112907520 -951350400 4281076800 -11237826600 17758540800 -16635041280 8506555200 -1829084400
-600600 47567520 -951350400 8245036800 -37875637800 101001700800 -161602721280 152907955200 -78843164400 17071454400
2522520 -208107900 4281076800 -37875637800 176752976400 -477233036280 771285715200 -735869534400 382086104400 -83223340200
-6306300 535134600 -11237826600 101001700800 -477233036280 1301544644400 -2121035716800 2037792556800 -1064382719400 233025352560
9609600 -832431600 17758540800 -161602721280 771285715200 -2121035716800 3480673996800 -3363975014400 1766086882560 -388375587600
-8751600 770140800 -16635041280 152907955200 -735869534400 2037792556800 -3363975014400 3267861442560 -1723286307600 380449555200
4375800 -389883780 8506555200 -78843164400 382086104400 -1064382719400 1766086882560 -1723286307600 912328045200 -202113826200
-923780 83140200 -1829084400 17071454400 -83223340200 233025352560 -388375587600 380449555200 -202113826200 44914183600
Re-inverted:
1 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10
1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11
1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12
1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13
1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14
1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15
1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16
1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17
1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18
1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19