Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Hilbert matrix cross product its inversion is a good demonstration of floating point versus rational arbitrary precision
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
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
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;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.