Bernstein basis polynomials
# 20230601 Raku programming solution
sub toBern2 { [ @_[0], @_[0]+@_[1]/2, @_[0..2].sum ] }
sub toBern3 (@a) { @a[0], @a[0]+@a[1]/3, @a[0]+@a[1]*2/3+@a[2]/3, @a[0..3].sum }
sub evalBern-N (@b, $t) { # uses de Casteljau's algorithm
my ($s, @bern) = 1 - $t, @b.Slip;
while ( @bern.elems > 2 ) {
@bern = @bern.rotor(2 => -1).map: { $s * .[0] + $t * .[1] };
}
return $s * @bern[0] + $t * @bern[1]
}
sub evaluations (@m, @b, $x) {
my $m = ([o] map { $_ + $x * * }, @m)(0); # Horner's rule
return "p({$x.fmt: '%.2f'}) = { evalBern-N @b, $x } (mono $m)";
}
sub bern2to3 (@a) { @a[0], @a[0]/3+@a[1]*2/3, @a[1]*2/3+@a[2]/3, @a[2] }
my (@pm,@qm,@rm) := ([1, 0, 0], [1, 2, 3], [1, 2, 3, 4]);
say "Subprogram(1) examples:";
my (@pb2,@qb2) := (@pm,@qm).map: { toBern2 $_ };
say "mono [{.[0]}] --> bern [{.[1]}]" for (@pm,@pb2,@qm,@qb2).rotor: 2;
say "\nSubprogram(2) examples:";
for (@pm,@pb2,@qm,@qb2).rotor(2) X (0.25,7.5) -> [[@m,@b], $x] {
say evaluations @m, @b, $x
}
say "\nSubprogram(3) examples:";
.push(0) for (@pm,@qm);
my (@pb3,@qb3,@rb3) := (@pm,@qm,@rm).map: { toBern3 $_ };
say "mono [{.[0]}] --> bern [{.[1]}]" for (@pm,@pb3,@qm,@qb3,@rm,@rb3).rotor: 2;
say "\nSubprogram(4) examples:";
for (@pm,@pb3,@qm,@qb3,@rm,@rb3).rotor(2) X (0.25,7.5) -> [[@m,@b], $x] {
say evaluations @m, @b, $x
}
say "\nSubprogram(5) examples:";
my (@pc,@qc) := (@pb2,@qb2).map: { bern2to3 $_ };
say "mono [{.[1]}] --> bern [{.[0]}]" for (@pc,@pb2,@qc,@qb2).rotor: 2;
Output:
Subprogram(1) examples:
mono [1 0 0] --> bern [1 1 1]
mono [1 2 3] --> bern [1 2 6]
Subprogram(2) examples:
p(0.25) = 1 (mono 1)
p(7.50) = 1 (mono 1)
p(0.25) = 1.6875 (mono 1.6875)
p(7.50) = 184.75 (mono 184.75)
Subprogram(3) examples:
mono [1 0 0 0] --> bern [1 1 1 1]
mono [1 2 3 0] --> bern [1 1.666667 3.333333 6]
mono [1 2 3 4] --> bern [1 1.666667 3.333333 10]
Subprogram(4) examples:
p(0.25) = 1 (mono 1)
p(7.50) = 1 (mono 1)
p(0.25) = 1.6875 (mono 1.6875)
p(7.50) = 184.75 (mono 184.75)
p(0.25) = 1.75 (mono 1.75)
p(7.50) = 1872.25 (mono 1872.25)
Subprogram(5) examples:
mono [1 1 1] --> bern [1 1 1 1]
mono [1 2 6] --> bern [1 1.666667 3.333333 6]
Last updated