for (  [1, 3, 5], # Test Matrices
       [2, 4, 7],
       [1, 1, 0]
    ),
    (  [11,  9, 24,  2],
       [ 1,  5,  2,  6],
       [ 3, 17, 18,  1],
       [ 2,  5,  7,  1]
    )
    -> @test {
    say-it 'A Matrix', @test;
    say-it( .[0], @(.[1]) ) for 'P Matrix', 'Aʼ Matrix', 'L Matrix', 'U Matrix' Z, lu @test;
}

sub lu (@a) {
    die unless @a.&is-square;
    my $n  = @a;
    my @P  = pivotize @a;
    my @Aʼ = mmult @P, @a;
    my @L  = matrix-ident $n;
    my @U  = matrix-zero  $n;
    for ^$n X ^$n -> ($i,$j) {
        if $j ≥ $i { @U[$i;$j] =  @Aʼ[$i;$j] - [+] map { @U[$_;$j] × @L[$i;$_] }, ^$i              }
        else       { @L[$i;$j] = (@Aʼ[$i;$j] - [+] map { @U[$_;$j] × @L[$i;$_] }, ^$j) / @U[$j;$j] }
    }
    @P, @Aʼ, @L, @U;
}

sub pivotize (@m) {
    my $size = @m;
    my @id   = matrix-ident $size;
    for ^$size -> $i {
        my $max = @m[$i;$i];
        my $row = $i;
        for $i ..^ $size -> $j {
            if @m[$j;$i] > $max {
                $max = @m[$j;$i];
                $row = $j;
            }
        }
        @id[$row, $i] = @id[$i, $row] if $row != $i;
    }
    @id
}

sub is-square (@m) { so @m == all @m }

sub matrix-zero ($n, $m = $n) { map { [ flat 0 xx $n ] }, ^$m }

sub matrix-ident ($n) { map { [ flat 0 xx $_, 1, 0 xx $n - 1 - $_ ] }, ^$n }

sub mmult(@a,@b) {
    my @p;
    for ^@a X ^@b[0] -> ($r, $c) {
        @p[$r;$c] += @a[$r;$_] × @b[$_;$c] for ^@b;
    }
    @p
}

sub rat-int ($num) {
    return $num unless $num ~~ Rat;
    return $num.narrow if $num.narrow ~~ Int;
    $num.nude.join: '/';

}

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