Hilbert matrix cross product its inversion is a good demonstration of floating point versus rational arbitrary precision

raku rationale-matrix.p6 --timings-only 8

Created 'Num' Hilbert of 8x8: 0s
Created 'Rat' Hilbert of 8x8: 0.0469023s
Starting 'Num' Hilbert test at 2020-08-07T17:32:09.217803+02:00
'Num' Hilbert inverse calculation: 8.0857941s
'Num' Hilbert (is identity? False): 8.0857941s
Starting 'Rat' Hilbert test at 2020-08-07T17:32:17.303597+02:00
'Rat' Hilbert inverse calculation: 5.3155128s
'Rat' Hilbert (is identity? True): 5.331197s

raku rationale-matrix.p6 5

raku rationale-matrix.p6 5


Num
         1       0.5  0.333333      0.25       0.2
       0.5  0.333333      0.25       0.2  0.166667
  0.333333      0.25       0.2  0.166667  0.142857
      0.25       0.2  0.166667  0.142857     0.125
       0.2  0.166667  0.142857     0.125  0.111111

            1             0            0             0            0
            0             1            0  -7.27596e-12  1.81899e-12
  2.84217e-14  -6.82121e-13            1  -3.63798e-12            0
  1.42109e-14  -2.27374e-13  2.72848e-12             1            0
            0  -2.27374e-13  9.09495e-13  -1.81899e-12            1


Rat
         1       0.5  0.333333      0.25       0.2
       0.5  0.333333      0.25       0.2  0.166667
  0.333333      0.25       0.2  0.166667  0.142857
      0.25       0.2  0.166667  0.142857     0.125
       0.2  0.166667  0.142857     0.125  0.111111

  1  0  0  0  0
  0  1  0  0  0
  0  0  1  0  0
  0  0  0  1  0
  0  0  0  0  1

rationale-matrique.raku

use Math::Matrix :ALL;

my %TYPES = :Num(1e0), :Rat(1);
subset FractionalRepresentation of Str where {%TYPES{$^t}:exists};
sub generate-hilbert($n, $type) {
    my @hilbert = [ [] xx $n ];
    for 1..$n -> $i {
        for 1..$n -> $j {
            @hilbert[$i-1;$j-1] = %TYPES{$type} / ($i + $j - 1);
        }
    }
    @hilbert
}
sub inner-hilbert(Int $n, FractionalRepresentation $type, $show-timings) {
    my $start = DateTime.now;
    my $hilbert = Math::Matrix.new: generate-hilbert($n,$type);
    say "Created '$type' Hilbert of {$n}x{$n}: {DateTime.now - $start}s"
        if $show-timings;
    $hilbert
}

sub MAIN(Int $size, Bool :$timings-only = False, Bool :$timings = False) {
    my $show-timings = $timings || $timings-only;
    my $num-bert = inner-hilbert $size, <Num>, $show-timings;
    my $rat-bert = inner-hilbert $size, <Rat>, $show-timings;

    for [$num-bert, <Num>], [$rat-bert, <Rat>] -> [$hilbert, $type] {
        say "Starting '$type' Hilbert test at " ~ my $start = DateTime.now
            if $show-timings;
        my $inverse = $hilbert ** -1;
        say "'$type' Hilbert inverse calculation: {DateTime.now - $start}s"
            if $show-timings;
        my $identity = $hilbert.dot-product($inverse);
        say "'$type' Hilbert (is identity? {$identity.is-identity}): {DateTime.now - $start}s"
            if $show-timings;
        say ("\n" xx 2) ~ "$type\n" ~ $hilbert.gist;
        say ("\n" xx 1) ~ $identity.gist
            unless $timings-only;
    }
}