Gauss-Jordan matrix inversion

Uses bits and pieces from other tasks, Reduced row echelon form 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 );
}

Output:

 ==================== 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

Last updated