Gaussian elimination

Gaussian elimination results in a matrix in row echelon form. Gaussian elimination with back-substitution (also known as Gauss-Jordan elimination) results in a matrix in reduced row echelon form. That being the case, we can reuse much of the code from the Reduced row echelon form task. Raku stores and does calculations on decimal numbers within its limit of precision using Rational numbers by default, meaning the calculations are exact.

sub gauss-jordan-solve (@a, @b) {
    @b.kv.map: { @a[$^k].append: $^v };
    @a.&rref[*]»[*-1];
}

# reduced row echelon form
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;
    return $num.narrow if $num.narrow ~~ Int;
    $num.nude.join: '/';
}

sub say-it ($message, @array, $fmt = " %8s") {
    say "\n$message";
    $_».&rat-or-int.fmt($fmt).put for @array;
}

my @a = (
    [ 1.00, 0.00, 0.00,  0.00,  0.00,   0.00 ],
    [ 1.00, 0.63, 0.39,  0.25,  0.16,   0.10 ],
    [ 1.00, 1.26, 1.58,  1.98,  2.49,   3.13 ],
    [ 1.00, 1.88, 3.55,  6.70, 12.62,  23.80 ],
    [ 1.00, 2.51, 6.32, 15.88, 39.90, 100.28 ],
    [ 1.00, 3.14, 9.87, 31.01, 97.41, 306.02 ],
);
my @b = ( -0.01, 0.61, 0.91, 0.99, 0.60, 0.02 );

say-it 'A matrix:', @a, "%6.2f";
say-it 'or, A in exact rationals:', @a;
say-it 'B matrix:', @b, "%6.2f";
say-it 'or, B in exact rationals:', @b;
say-it 'x matrix:', (my @gj = gauss-jordan-solve @a, @b), "%16.12f";
say-it 'or, x in exact rationals:', @gj, "%28s";

Output:

A matrix:
  1.00   0.00   0.00   0.00   0.00   0.00
  1.00   0.63   0.39   0.25   0.16   0.10
  1.00   1.26   1.58   1.98   2.49   3.13
  1.00   1.88   3.55   6.70  12.62  23.80
  1.00   2.51   6.32  15.88  39.90 100.28
  1.00   3.14   9.87  31.01  97.41 306.02

or, A in exact rationals:
        1         0         0         0         0         0
        1    63/100    39/100       1/4      4/25      1/10
        1     63/50     79/50     99/50   249/100   313/100
        1     47/25     71/20     67/10    631/50     119/5
        1   251/100    158/25    397/25    399/10   2507/25
        1    157/50   987/100  3101/100  9741/100  15301/50

B matrix:
 -0.01
  0.61
  0.91
  0.99
  0.60
  0.02

or, B in exact rationals:
   -1/100
   61/100
   91/100
   99/100
      3/5
     1/50

x matrix:
 -0.010000000000
  1.602790394502
 -1.613203059906
  1.245494121371
 -0.490989719585
  0.065760696175

or, x in exact rationals:
                      -1/100
   655870882787/409205648497
  -660131804286/409205648497
   509663229635/409205648497
  -200915766608/409205648497
    26909648324/409205648497

Last updated