Bernstein basis polynomials: Difference between revisions
Add C# implementation
mNo edit summary |
(Add C# implementation) |
||
| (18 intermediate revisions by 7 users not shown) | |||
Line 1:
{{
The <math>n + 1</math> [[wp:Bernstein_polynomial|''Bernstein basis polynomials'']] of degree <math>n</math> are defined as
:<math>b_{k,n}(t) = \binom{n}{k} t^{k} \left( 1 - t \right)^{n - k},\quad k = 0,\ldots,n</math>
Line 21:
[[#ALGOL_60|ALGOL 60]] and [[#Python|Python]] implementations are provided as initial examples. The latter does the optional monomial-basis evaluations.
You can use the following algorithms. They are written in unambiguous Algol 60 reference language instead of a made up pseudo-language. The [[#ALGOL_60|ALGOL 60]] example was necessary to check my work, but these reference versions are in the actual standard language designed for the printing of algorithms.
'''procedure''' ''tobern2'' (''a0'', ''a1'', ''a2'', ''b0'', ''b1'', ''b2'');
'''value''' ''a0'', ''a1'', ''a2''; '''comment''' pass by value;
Line 333:
=={{header|ALGOL 68}}==
Algol 68 doesn't have procedure overloading but does have monadic and dyadic operator overloading and user defined operator symbols.
<syntaxhighlight lang="algol68">
BEGIN # Bernstein Basis Polynomials - translated from the Algol 60 #
# Note the "Subrptoghram (n) designations are as specified in the #
# task #
MODE BERNTWO = STRUCT( STRING label, REAL b0, b1, b2 );
MODE BERNTHREE = STRUCT( STRING label, REAL b0, b1, b2, b3 );
PRIO EVAL = 5; # set the priority of the dyadic EVAL #
# Subprogram (1): transform monomial coefficients #
# a0, a1, a2 to Bernstein coefficients b0, b1, b2 #
PROC tobern2 = ( REAL a0, a1, a2, REF
BEGIN
b0 OF b := a0;
b1 OF b := a0 + ((1/2) * a1);
b2 OF b := a0 + a1 + a2
END # tobern2 # ;
Line 350 ⟶ 356:
# Bernstein coefficients b0, b1, b2. Use de Casteljau's #
# algorithm #
BEGIN
REAL s = 1 - t;
REAL b01 = (s * b0 OF b) + (t * b1 OF b);
REAL b12 = (s * b1 OF b) + (t * b2 OF b);
(s * b01) + (t * b12)
END #
# Subprogram (3): transform monomial coefficients #
# a0, a1, a2, a3 to Bernstein coefficients b0, b1, #
# b2, b3 #
PROC tobern3 = ( REAL a0, a1, a2, a3, REF
BEGIN
b0 OF b := a0;
b1 OF b := a0 + ((1/3) * a1);
b2 OF b := a0 + ((2/3) * a1) + ((1/3) * a2);
b3 OF b := a0 + a1 + a2 + a3
END # tobern3 # ;
Line 372 ⟶ 378:
# with Bernstein coefficients b0, b1, b2, b3. Use #
# de Casteljau's algorithm #
BEGIN
REAL s = 1 - t;
REAL b01 = (s * b0 OF b) + (t * b1 OF b);
REAL b12 = (s * b1 OF b) + (t * b2 OF b);
REAL b23 = (s * b2 OF b) + (t * b3 OF b);
REAL b012 = (s * b01) + (t * b12);
REAL b123 = (s * b12) + (t * b23);
(s * b012) + (t * b123)
END #
# Subprogram (5): transform the quadratic Bernstein #
# coefficients q0, q1, q2 to the cubic Bernstein #
# coefficients c0, c1, c2, c3; #
OP TOBERNTHREE = ( BERNTWO q )BERNTHREE:
BEGIN
c
END # TOBERNTHREE # ;
BEGIN
Line 414 ⟶ 422:
END # f # ;
OP SHOWEVAL = ( BERNTWO p, REAL x )VOID:
print( ( " ", label OF p, " ( ", f( x ), " ) = ", f( p EVAL x ), newline ) );
"( " + f( b0 OF p ) + ", " + f( b1 OF p ) + ", " + f( b2 OF p ) + " )";
# returns a string representation of the values of p #
OP TOSTRING = ( BERNTHREE p )STRING:
"( " + f( b0 OF p ) + ", " + f( b1 OF p ) + ", " + f( b2 OF p ) + ", " + f( b3 OF p ) + " )";
# returns a string representation of the values of p #
OP TOSTRING = ( []REAL p )STRING:
BEGIN
STRING result := "(", separator := "";
FOR i FROM LWB p TO UPB p DO
result +:= separator + " " + f( p[ i ] );
separator := ","
OD;
result + " )"
END # TOSTRING # ;
BERNTWO q2 := BERNTWO ( "q", 0, 0, 0 );
BERNTHREE p3 := BERNTHREE( "p", 0, 0, 0, 0 );
BERNTHREE q3 := BERNTHREE( "q", 0, 0, 0, 0 );
BERNTHREE r3 := BERNTHREE( "r", 0, 0, 0, 0 );
REAL p0m = 1, p1m = 0, p2m = 0;
Line 430 ⟶ 457:
REAL r0m = 1, r1m = 2, r2m = 3, r3m = 4;
tobern2( p0m, p1m, p2m,
tobern2( q0m, q1m, q2m,
print( ( "Subprogram (1) examples:", newline ) );
print( ( " mono
print( ( " mono
print( ( "Subprogram (2) examples:", newline ) );
q2 SHOWEVAL 0.25;
tobern3( p0m, p1m, p2m, 0,
tobern3( q0m, q1m, q2m, 0,
tobern3( r0m, r1m, r2m, r3m,
print( ( "Subprogram (3) examples:", newline ) );
print( ( " mono
print(
print( ( " mono ", TOSTRING []REAL( r0m,
print( ( "Subprogram (4) examples:", newline ) );
p3 SHOWEVAL 7.50;
q3 SHOWEVAL 0.25;
r3 SHOWEVAL 0.25;
r3 SHOWEVAL 7.50;
print( ( "Subprogram (5) examples:", newline ) );
print( ( " bern
print( ( " bern ", TOSTRING q2, "
END
Line 563 ⟶ 541:
(* Some macros that are convenient for type-safe floating point
without actually knowing which floating-point type it will
be. These macros will work if, in the context where they appear,
the typechecker can infer which floating-point type is meant. *)
macdef one_third = g0i2f 1 / g0i2f 3
macdef one_half = g0i2f 1 / g0i2f 2
Line 572 ⟶ 550:
(* We will try to achieve the type-safety of distinguishing between
monomial and Bernstein bases at compile time, but without overhead
at runtime.
languages simply by giving the types different names but the same
fields. But then you cannot easily have TYPE-SAFE subprograms that
work on EITHER type. You could have that with a type hierarchy, but
then, typically, you start getting overhead from runtime
polymorphism. You could do it by nesting record types within record
types, but then you start needing the overhead of pointers to the
different nested types. Furthermore, these approaches are all
essentially WORKAROUNDS, rather than direct solutions to the
problem, which is "how to distinguish between a tuple representing
coefficients in ONE basis from the exact same kind of tuple
representing coefficients in A DIFFERENT basis, while also treating
both as coefficients in SOME basis, such as for addition or scalar
multiplication". Here we do that by attaching an object to the
tuple AT TYPECHECKING TIME that is entirely left out after
typechecking is complete.) *)
tkindef monoknd = "rosettacode_monomial_poly"
tkindef bernknd = "rosettacode_bernstein_poly"
Line 597 ⟶ 590:
mono2bern_degree2 with Bernstein coefficients (instead of monomial
coefficients) as the argument, then, AT COMPILE TIME, you get
/some_path/bernstein_poly_task.dats: 1051(line=32, offs=31) -- 1052(line=32, offs=32): warning(3): the constraint [S2Eeqeq(S2Eextkind(rosettacode_bernstein_poly); S2Eextkind(rosettacode_monomial_poly))] cannot be translated into a form accepted by the constraint solver.
Line 671 ⟶ 664:
end
(* Let us make it easy to print out coefficients.
print_poly_degree2 operates TYPE-SAFELY on ANY object of type
"poly", without regard to whether it represents monomial
coefficients or Bernstein coefficients. In C we might do such a
thing by casting a pointer of one type to a pointer of another
type, but that is not a TYPE-SAFE operation. *)
fn {coefknd : tkind}
print_poly_degree2 {polyknd : tkind}
(coefs : poly_degree2 (polyknd, coefknd))
: void =
let
Line 679 ⟶ 678:
macdef outf = stdout_ref
in
fprint_val<string> (outf, "
fprint_val<g0float coefknd> (outf, a0);
fprint_val<string> (outf, ", ");
Line 689 ⟶ 688:
fn {coefknd : tkind}
print_poly_degree3 {polyknd : tkind}
(coefs : poly_degree3 (polyknd, coefknd))
: void =
let
Line 711 ⟶ 695:
macdef outf = stdout_ref
in
fprint_val<string> (outf, "
fprint_val<g0float coefknd> (outf, a0);
fprint_val<string> (outf, ", ");
Line 720 ⟶ 704:
fprint_val<g0float coefknd> (outf, a3);
fprint_val<string> (outf, ")")
end
fn {coefknd : tkind}
print_mono_degree2 (coefs : poly_degree2 (monoknd, coefknd))
: void =
begin
fprint_val<string> (stdout_ref, "mono ");
print_poly_degree2 (coefs)
end
fn {coefknd : tkind}
print_bern_degree2 (coefs : poly_degree2 (bernknd, coefknd))
: void =
begin
fprint_val<string> (stdout_ref, "bern ");
print_poly_degree2 (coefs)
end
fn {coefknd : tkind}
print_mono_degree3 (coefs : poly_degree3 (monoknd, coefknd))
: void =
begin
fprint_val<string> (stdout_ref, "mono ");
print_poly_degree3 (coefs)
end
Line 725 ⟶ 733:
print_bern_degree3 (coefs : poly_degree3 (bernknd, coefknd))
: void =
fprint_val<string> (stdout_ref, "bern ");
print_poly_degree3 (coefs)
end
Line 923 ⟶ 921:
bern (1.000000, 1.000000, 1.000000) --> bern (1.000000, 1.000000, 1.000000, 1.000000)
bern (1.000000, 2.000000, 6.000000) --> bern (1.000000, 1.666667, 3.333333, 6.000000)</pre>
=={{header|C}}==
{{trans|Wren}}
<syntaxhighlight lang="c">#include <stdio.h>
typedef double dbl;
void tobern2(dbl a0, dbl a1, dbl a2, dbl *b0, dbl *b1, dbl *b2) {
*b0 = a0;
*b1 = a0 + a1 / 2;
*b2 = a0 + a1 + a2;
}
/* uses de Casteljau’s algorithm */
dbl evalbern2(dbl b0, dbl b1, dbl b2, dbl t) {
dbl s = 1.0 - t;
dbl b01 = s * b0 + t * b1;
dbl b12 = s * b1 + t * b2;
return s * b01 + t * b12;
}
void tobern3(dbl a0, dbl a1, dbl a2, dbl a3, dbl *b0, dbl *b1, dbl *b2, dbl *b3) {
*b0 = a0;
*b1 = a0 + a1 / 3;
*b2 = a0 + a1 * 2/3 + a2 / 3;
*b3 = a0 + a1 + a2 + a3;
}
/* uses de Casteljau’s algorithm */
dbl evalbern3(dbl b0, dbl b1, dbl b2, dbl b3, dbl t) {
dbl s = 1 - t;
dbl b01 = s * b0 + t * b1;
dbl b12 = s * b1 + t * b2;
dbl b23 = s * b2 + t * b3;
dbl b012 = s * b01 + t * b12;
dbl b123 = s * b12 + t * b23;
return s * b012 + t * b123;
}
void bern2to3(dbl q0, dbl q1, dbl q2, dbl *c0, dbl *c1, dbl *c2, dbl *c3) {
*c0 = q0;
*c1 = q0 / 3 + q1 * 2/3;
*c2 = q1 * 2/3 + q2 / 3;
*c3 = q2;
}
/* uses Horner's rule */
dbl evalmono2(dbl a0, dbl a1, dbl a2, dbl t) {
return a0 + (t * (a1 + (t * a2)));
}
/* uses Horner's rule */
dbl evalmono3(dbl a0, dbl a1, dbl a2, dbl a3, dbl t) {
return a0 + (t * (a1 + (t * (a2 + (t * a3)))));
}
int main() {
dbl pm0 = 1, pm1 = 0, pm2 = 0, pm3 = 0, pb0, pb1, pb2, pb3, pc0, pc1, pc2, pc3;
dbl qm0 = 1, qm1 = 2, qm2 = 3, qm3 = 0, qb0, qb1, qb2, qb3, qc0, qc1, qc2, qc3;
dbl rm0 = 1, rm1 = 2, rm2 = 3, rm3 = 4, rb0, rb1, rb2, rb3;
dbl x, y, m;
const char *fmt;
printf("Subprogram(1) examples:\n");
tobern2(pm0, pm1, pm2, &pb0, &pb1, &pb2);
tobern2(qm0, qm1, qm2, &qb0, &qb1, &qb2);
fmt = "mono {%g, %g, %g} --> bern {%g, %g, %g}\n";
printf(fmt, pm0, pm1, pm2, pb0, pb1, pb2);
printf(fmt, qm0, qm1, qm2, qb0, qb1, qb2);
printf("\nSubprogram(2) examples:\n");
x = 0.25;
y = evalbern2(pb0, pb1, pb2, x);
m = evalmono2(pm0, pm1, pm2, x);
printf("p(%4.2f) = %g (mono %g)\n", x, y, m);
x = 7.5;
y = evalbern2(pb0, pb1, pb2, x);
m = evalmono2(pm0, pm1, pm2, x);
printf("p(%4.2f) = %g (mono %g)\n", x, y, m);
x = 0.25;
y = evalbern2(qb0, qb1, qb2, x);
m = evalmono2(qm0, qm1, qm2, x);
printf("q(%4.2f) = %g (mono %g)\n", x, y, m);
x = 7.5;
y = evalbern2(qb0, qb1, qb2, x);
m = evalmono2(qm0, qm1, qm2, x);
printf("q(%4.2f) = %g (mono %g)\n", x, y, m);
printf("\nSubprogram(3) examples:\n");
tobern3(pm0, pm1, pm2, pm3, &pb0, &pb1, &pb2, &pb3);
tobern3(qm0, qm1, qm2, qm3, &qb0, &qb1, &qb2, &qb3);
tobern3(rm0, rm1, rm2, rm3, &rb0, &rb1, &rb2, &rb3);
fmt = "mono {%g, %g, %g, %g} --> bern {%0.14g, %0.14g, %0.14g, %0.14g}\n";
printf(fmt, pm0, pm1, pm2, pm3, pb0, pb1, pb2, pb3);
printf(fmt, qm0, qm1, qm2, qm3, qb0, qb1, qb2, qb3);
printf(fmt, rm0, rm1, rm2, rm3, rb0, rb1, rb2, rb3);
printf("\nSubprogram(4) examples:\n");
x = 0.25;
y = evalbern3(pb0, pb1, pb2, pb3, x);
m = evalmono3(pm0, pm1, pm2, pm3, x);
printf("p(%4.2f) = %g (mono %g)\n", x, y, m);
x = 7.5;
y = evalbern3(pb0, pb1, pb2, pb3, x);
m = evalmono3(pm0, pm1, pm2, pm3, x);
printf("p(%4.2f) = %g (mono %g)\n", x, y, m);
x = 0.25;
y = evalbern3(qb0, qb1, qb2, qb3, x);
m = evalmono3(qm0, qm1, qm2, qm3, x);
printf("q(%4.2f) = %g (mono %g)\n", x, y, m);
x = 7.5;
y = evalbern3(qb0, qb1, qb2, qb3, x);
m = evalmono3(qm0, qm1, qm2, qm3, x);
printf("q(%4.2f) = %g (mono %g)\n", x, y, m);
x = 0.25;
y = evalbern3(rb0, rb1, rb2, rb3, x);
m = evalmono3(rm0, rm1, rm2, rm3, x);
printf("r(%4.2f) = %g (mono %g)\n", x, y, m);
x = 7.5;
y = evalbern3(rb0, rb1, rb2, rb3, x);
m = evalmono3(rm0, rm1, rm2, rm3, x);
printf("r(%4.2f) = %g (mono %g)\n", x, y, m);
printf("\nSubprogram(5) examples:\n");
tobern2(pm0, pm1, pm2, &pb0, &pb1, &pb2);
tobern2(qm0, qm1, qm2, &qb0, &qb1, &qb2);
bern2to3(pb0, pb1, pb2, &pc0, &pc1, &pc2, &pc3);
bern2to3(qb0, qb1, qb2, &qc0, &qc1, &qc2, &qc3);
fmt = "mono {%g, %g, %g} --> bern {%0.14g, %0.14g, %0.14g, %0.14g}\n";
printf(fmt, pb0, pb1, pb2, pc0, pc1, pc2, pc3);
printf(fmt, qb0, qb1, qb2, qc0, qc1, qc2, qc3);
return 0;
}</syntaxhighlight>
{{out}}
<pre>
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)
q(0.25) = 1.6875 (mono 1.6875)
q(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.6666666666667, 3.3333333333333, 6}
mono {1, 2, 3, 4} --> bern {1, 1.6666666666667, 3.3333333333333, 10}
Subprogram(4) examples:
p(0.25) = 1 (mono 1)
p(7.50) = 1 (mono 1)
q(0.25) = 1.6875 (mono 1.6875)
q(7.50) = 184.75 (mono 184.75)
r(0.25) = 1.75 (mono 1.75)
r(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.6666666666667, 3.3333333333333, 6}
</pre>
=={{header|C#}}==
{{trans|Go}}
<syntaxhighlight lang="C#">
using System;
class Program
{
static double[] ToBern2(double[] a)
{
return new double[] { a[0], a[0] + a[1] / 2, a[0] + a[1] + a[2] };
}
static double EvalBern2(double[] b, double t)
{
double s = 1.0 - t;
double b01 = s * b[0] + t * b[1];
double b12 = s * b[1] + t * b[2];
return s * b01 + t * b12;
}
static double[] ToBern3(double[] a)
{
return new double[] { a[0], a[0] + a[1] / 3, a[0] + a[1] * 2 / 3 + a[2] / 3, a[0] + a[1] + a[2] + a[3] };
}
static double EvalBern3(double[] b, double t)
{
double s = 1.0 - t;
double b01 = s * b[0] + t * b[1];
double b12 = s * b[1] + t * b[2];
double b23 = s * b[2] + t * b[3];
double b012 = s * b01 + t * b12;
double b123 = s * b12 + t * b23;
return s * b012 + t * b123;
}
static double[] Bern2To3(double[] q)
{
return new double[] { q[0], q[0] / 3 + q[1] * 2 / 3, q[1] * 2 / 3 + q[2] / 3, q[2] };
}
static double EvalMono2(double[] a, double t)
{
return a[0] + (t * (a[1] + (t * a[2])));
}
static double EvalMono3(double[] a, double t)
{
return a[0] + (t * (a[1] + (t * (a[2] + (t * a[3])))));
}
static void Main(string[] args)
{
double[] pm = { 1, 0, 0 };
double[] qm = { 1, 2, 3 };
double[] rm = { 1, 2, 3, 4 };
double x, y, m;
Console.WriteLine("Subprogram(1) examples:");
var pb2 = ToBern2(pm);
var qb2 = ToBern2(qm);
Console.WriteLine($"mono [{string.Join(", ", pm)}] --> bern [{string.Join(", ", pb2)}]");
Console.WriteLine($"mono [{string.Join(", ", qm)}] --> bern [{string.Join(", ", qb2)}]");
Console.WriteLine("\nSubprogram(2) examples:");
x = 0.25;
y = EvalBern2(pb2, x);
m = EvalMono2(pm, x);
Console.WriteLine($"p({x:F2}) = {y:G14} (mono {m:G14})");
x = 7.5;
y = EvalBern2(pb2, x);
m = EvalMono2(pm, x);
Console.WriteLine($"p({x:F2}) = {y:G14} (mono {m:G14})");
x = 0.25;
y = EvalBern2(qb2, x);
m = EvalMono2(qm, x);
Console.WriteLine($"q({x:F2}) = {y:G14} (mono {m:G14})");
x = 7.5;
y = EvalBern2(qb2, x);
m = EvalMono2(qm, x);
Console.WriteLine($"q({x:F2}) = {y:G14} (mono {m:G14})");
Console.WriteLine("\nSubprogram(3) examples:");
var pb3 = ToBern3(new double[] { pm[0], pm[1], pm[2], 0 });
var qb3 = ToBern3(new double[] { qm[0], qm[1], qm[2], 0 });
var rb3 = ToBern3(rm);
Console.WriteLine($"mono [{string.Join(", ", pm)}] --> bern [{string.Join(", ", pb3)}]");
Console.WriteLine($"mono [{string.Join(", ", qm)}] --> bern [{string.Join(", ", qb3)}]");
Console.WriteLine($"mono [{string.Join(", ", rm)}] --> bern [{string.Join(", ", rb3)}]");
Console.WriteLine("\nSubprogram(4) examples:");
x = 0.25;
y = EvalBern3(pb3, x);
m = EvalMono3(new double[] { pm[0], pm[1], pm[2], 0 }, x);
Console.WriteLine($"p({x:F2}) = {y:G14} (mono {m:G14})");
x = 7.5;
y = EvalBern3(pb3, x);
m = EvalMono3(new double[] { pm[0], pm[1], pm[2], 0 }, x);
Console.WriteLine($"p({x:F2}) = {y:G14} (mono {m:G14})");
x = 0.25;
y = EvalBern3(qb3, x);
m = EvalMono3(new double[] { qm[0], qm[1], qm[2], 0 }, x);
Console.WriteLine($"q({x:F2}) = {y:G14} (mono {m:G14})");
x = 7.5;
y = EvalBern3(qb3, x);
m = EvalMono3(new double[] { qm[0], qm[1], qm[2], 0 }, x);
Console.WriteLine($"q({x:F2}) = {y:G14} (mono {m:G14})");
x = 0.25;
y = EvalBern3(rb3, x);
m = EvalMono3(rm, x);
Console.WriteLine($"r({x:F2}) = {y:G14} (mono {m:G14})");
x = 7.5;
y = EvalBern3(rb3, x);
m = EvalMono3(rm, x);
Console.WriteLine($"r({x:F2}) = {y:G14} (mono {m:G14})");
Console.WriteLine("\nSubprogram(5) examples:");
var pc = Bern2To3(pb2);
var qc = Bern2To3(qb2);
Console.WriteLine($"bern [{string.Join(", ", pb2)}] --> bern3 [{string.Join(", ", pc)}]");
Console.WriteLine($"bern [{string.Join(", ", qb2)}] --> bern3 [{string.Join(", ", qc)}]");
}
}
</syntaxhighlight>
{{out}}
<pre>
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)
q(0.25) = 1.6875 (mono 1.6875)
q(7.50) = 184.75 (mono 184.75)
Subprogram(3) examples:
mono [1, 0, 0] --> bern [1, 1, 1, 1]
mono [1, 2, 3] --> bern [1, 1.66666666666667, 3.33333333333333, 6]
mono [1, 2, 3, 4] --> bern [1, 1.66666666666667, 3.33333333333333, 10]
Subprogram(4) examples:
p(0.25) = 1 (mono 1)
p(7.50) = 1 (mono 1)
q(0.25) = 1.6875 (mono 1.6875)
q(7.50) = 184.75 (mono 184.75)
r(0.25) = 1.75 (mono 1.75)
r(7.50) = 1872.25 (mono 1872.25)
Subprogram(5) examples:
bern [1, 1, 1] --> bern3 [1, 1, 1, 1]
bern [1, 2, 6] --> bern3 [1, 1.66666666666667, 3.33333333333333, 6]
</pre>
=={{header|C++}}==
<syntaxhighlight lang="c++">
#include <cstdint>
#include <iostream>
#include <vector>
std::string to_string(const std::vector<double>& list) {
std::string result = "[";
for ( uint64_t i = 0; i < list.size() - 1; ++i ) {
result += std::to_string(list[i]) + ", ";
}
result += std::to_string(list.back()) + "]";
return result;
}
// Subprogram (1)
std::vector<double> monomial_to_bernstein_degree2(const std::vector<double>& monomial) {
return std::vector<double>{ monomial[0],
monomial[0] + ( monomial[1] / 2.0 ),
monomial[0] + monomial[1] + monomial[2] };
}
// Subprogram (2)
double evaluate_bernstein_degree2(const std::vector<double>& bernstein, const double& t) {
// de Casteljau’s algorithm
const double s = 1 - t;
const double b01 = ( s * bernstein[0] ) + ( t * bernstein[1] );
const double b12 = ( s * bernstein[1] ) + ( t * bernstein[2] );
return ( s * b01 ) + ( t * b12 );
}
// Subprogram (3)
std::vector<double> monomial_to_bernstein_degree3(const std::vector<double>& monomial) {
return std::vector<double>{ monomial[0],
monomial[0] + ( monomial[1] / 3.0 ),
monomial[0] + ( 2.0 * monomial[1] / 3.0 ) + ( monomial[2] / 3.0 ),
monomial[0] + monomial[1] + monomial[2] + monomial[3] };
}
// Subprogram (4)
double evaluate_bernstein_degree3(const std::vector<double>& bernstein, const double& t) {
// de Casteljau’s algorithm
const double s = 1 - t;
const double b01 = ( s * bernstein[0] ) + ( t * bernstein[1] );
const double b12 = ( s * bernstein[1] ) + ( t * bernstein[2] );
const double b23 = ( s * bernstein[2] ) + ( t * bernstein[3] );
const double b012 = ( s * b01 ) + ( t * b12 );
const double b123 = ( s * b12 ) + ( t * b23 );
return ( s * b012 ) + ( t * b123 );
}
// Subprogram (5)
std::vector<double> bernstein_degree2_to_degree3(const std::vector<double>& bernstein) {
return std::vector<double>{ bernstein[0],
( bernstein[0] / 3.0 ) + ( 2.0 * bernstein[1] / 3.0 ),
( 2.0 * bernstein[1] / 3.0 ) + ( bernstein[2] / 3.0 ),
bernstein[2] };
}
double evaluate_monomial_degree2(const std::vector<double>& monomial, const double& t) {
// Horner’s rule
return monomial[0] + ( t * ( monomial[1] + ( t * monomial[2] ) ) );
}
double evaluate_monomial_degree3(const std::vector<double>& monomial, const double& t) {
// Horner’s rule
return monomial[0] + ( t * ( monomial[1] + ( t * ( monomial[2] + ( t * monomial[3] ) ) ) ) );
}
int main() {
/**
* For the following polynomials, use Subprogram (1) to find coefficients in the degree-2 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
*/
std::vector<double> pMonomial2 = { 1.0, 0.0, 0.0 };
std::vector<double> qMonomial2 = { 1.0, 2.0, 3.0 };
std::vector<double> pBernstein2 = monomial_to_bernstein_degree2(pMonomial2);
std::vector<double> qBernstein2 = monomial_to_bernstein_degree2(qMonomial2);
std::cout << "Subprogram (1) examples:" << std::endl;
std::cout << " monomial " + to_string(pMonomial2) + " --> bernstein " + to_string(pBernstein2) << std::endl;
std::cout << " monomial " + to_string(qMonomial2) + " --> bernstein " + to_string(qBernstein2) << std::endl;
/**
* Use Subprogram (2) to evaluate p(x) and q(x) at x = 0.25, 7.50. Display the results.
* Optionally also display results from evaluating in the original monomial basis.
*/
std::cout << "Subprogram (2) examples:" << std::endl;
for ( const double& x : { 0.25, 7.50 } ) {
std::cout << " p(" << x << ") = " << evaluate_bernstein_degree2(pBernstein2, x)
<< " ( mono: " << evaluate_monomial_degree2(pMonomial2, x) << " )" << std::endl;
}
for ( const double& x : { 0.25, 7.50 } ) {
std::cout << " q(" << x << ") = " << evaluate_bernstein_degree2(qBernstein2, x)
<< " ( mono: " << evaluate_monomial_degree2(qMonomial2, x) << " )" << std::endl;
}
/**
* For the following polynomials, use Subprogram (3) to find coefficients in the degree-3 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
* r(x) = 1 + 2x + 3x² + 4x³
*
* Display the results.
*/
std::vector<double> pMonomial3 = { 1.0, 0.0, 0.0, 0.0 };
std::vector<double> qMonomial3 = { 1.0, 2.0, 3.0, 0.0 };
std::vector<double> rMonomial3 = { 1.0, 2.0, 3.0, 4.0 };
std::vector<double> pBernstein3 = monomial_to_bernstein_degree3(pMonomial3);
std::vector<double> qBernstein3 = monomial_to_bernstein_degree3(qMonomial3);
std::vector<double> rBernstein3 = monomial_to_bernstein_degree3(rMonomial3);
std::cout << "Subprogram (3) examples:" << std::endl;
std::cout << " monomial " + to_string(pMonomial3) + " --> bernstein " + to_string(pBernstein3) << std::endl;
std::cout << " monomial " + to_string(qMonomial3) + " --> bernstein " + to_string(qBernstein3) << std::endl;
std::cout << " monomial " + to_string(rMonomial3) + " --> bernstein " + to_string(rBernstein3) << std::endl;
/**
* Use Subprogram (4) to evaluate p(x), q(x), and r(x) at x = 0.25, 7.50. Display the results.
* Optionally also display results from evaluating in the original monomial basis.
*/
std::cout << "Subprogram (4) examples:" << std::endl;
for ( const double& x : { 0.25, 7.50 } ) {
std::cout << " p(" << x << ") = " << evaluate_bernstein_degree3(pBernstein3, x)
<< " ( mono: " << evaluate_monomial_degree3(pMonomial3, x) << " )" << std::endl;
}
for ( const double& x : { 0.25, 7.50 } ) {
std::cout << " q(" << x << ") = " << evaluate_bernstein_degree3(qBernstein3, x)
<< " ( mono: " << evaluate_monomial_degree3(qMonomial3, x) << " )" << std::endl;
}
for ( const double& x : { 0.25, 7.50 } ) {
std::cout << " r(" << x << ") = " << evaluate_bernstein_degree3(rBernstein3, x)
<< " ( mono: " << evaluate_monomial_degree3(rMonomial3, x) << " )" << std::endl;
}
/**
* For the following polynomials, using the result of Subprogram (1) applied to the polynomial,
* use Subprogram (5) to get coefficients for the degree-3 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
*
* Display the results.
*/
std::cout << "Subprogram (5) examples:" << std::endl;
std::vector<double> pBernstein3a = bernstein_degree2_to_degree3(pBernstein2);
std::vector<double> qBernstein3a = bernstein_degree2_to_degree3(qBernstein2);
std::cout << " bernstein " + to_string(pBernstein2) + " --> bernstein " + to_string(pBernstein3a) << std::endl;
std::cout << " bernstein " + to_string(qBernstein2) + " --> bernstein " + to_string(qBernstein3a) << std::endl;
}
</syntaxhighlight>
{{ out }}
<pre>
Subprogram (1) examples:
monomial [1.000000, 0.000000, 0.000000] --> bernstein [1.000000, 1.000000, 1.000000]
monomial [1.000000, 2.000000, 3.000000] --> bernstein [1.000000, 2.000000, 6.000000]
Subprogram (2) examples:
p(0.25) = 1 ( mono: 1 )
p(7.5) = 1 ( mono: 1 )
q(0.25) = 1.6875 ( mono: 1.6875 )
q(7.5) = 184.75 ( mono: 184.75 )
Subprogram (3) examples:
monomial [1.000000, 0.000000, 0.000000, 0.000000] --> bernstein [1.000000, 1.000000, 1.000000, 1.000000]
monomial [1.000000, 2.000000, 3.000000, 0.000000] --> bernstein [1.000000, 1.666667, 3.333333, 6.000000]
monomial [1.000000, 2.000000, 3.000000, 4.000000] --> bernstein [1.000000, 1.666667, 3.333333, 10.000000]
Subprogram (4) examples:
p(0.25) = 1 ( mono: 1 )
p(7.5) = 1 ( mono: 1 )
q(0.25) = 1.6875 ( mono: 1.6875 )
q(7.5) = 184.75 ( mono: 184.75 )
r(0.25) = 1.75 ( mono: 1.75 )
r(7.5) = 1872.25 ( mono: 1872.25 )
Subprogram (5) examples:
bernstein [1.000000, 1.000000, 1.000000] --> bernstein [1.000000, 1.000000, 1.000000, 1.000000]
bernstein [1.000000, 2.000000, 6.000000] --> bernstein [1.000000, 1.666667, 3.333333, 6.000000]
</pre>
=={{header|FreeBASIC}}==
Line 1,086 ⟶ 1,584:
bern (1.00, 1.00, 1.00) --> bern (1.00, 1.00, 1.00, 1.00)
bern (1.00, 2.00, 6.00) --> bern (1.00, 1.67, 3.33, 6.00)</pre>
=={{header|Go}}==
{{trans|Wren}}
<syntaxhighlight lang="go">package main
import "fmt"
func toBern2(a []float64) []float64 {
return []float64{a[0], a[0] + a[1]/2, a[0] + a[1] + a[2]}
}
// uses de Casteljau's algorithm
func evalBern2(b []float64, t float64) float64 {
s := 1.0 - t
b01 := s*b[0] + t*b[1]
b12 := s*b[1] + t*b[2]
return s*b01 + t*b12
}
func toBern3(a []float64) []float64 {
b := make([]float64, 4)
b[0] = a[0]
b[1] = a[0] + a[1]/3
b[2] = a[0] + a[1]*2/3 + a[2]/3
b[3] = a[0] + a[1] + a[2] + a[3]
return b
}
// uses de Casteljau's algorithm
func evalBern3(b []float64, t float64) float64 {
s := 1.0 - t
b01 := s*b[0] + t*b[1]
b12 := s*b[1] + t*b[2]
b23 := s*b[2] + t*b[3]
b012 := s*b01 + t*b12
b123 := s*b12 + t*b23
return s*b012 + t*b123
}
func bern2to3(q []float64) []float64 {
c := make([]float64, 4)
c[0] = q[0]
c[1] = q[0]/3 + q[1]*2/3
c[2] = q[1]*2/3 + q[2]/3
c[3] = q[2]
return c
}
// uses Horner's rule
func evalMono2(a []float64, t float64) float64 {
return a[0] + (t * (a[1] + (t * a[2])))
}
// uses Horner's rule
func evalMono3(a []float64, t float64) float64 {
return a[0] + (t * (a[1] + (t * (a[2] + (t * a[3])))))
}
func main() {
pm := []float64{1, 0, 0}
qm := []float64{1, 2, 3}
rm := []float64{1, 2, 3, 4}
var x, y, m float64
fmt.Println("Subprogram(1) examples:")
pb2 := toBern2(pm)
qb2 := toBern2(qm)
fmt.Printf("mono %v --> bern %v\n", pm, pb2)
fmt.Printf("mono %v --> bern %v\n", qm, qb2)
fmt.Println("\nSubprogram(2) examples:")
x = 0.25
y = evalBern2(pb2, x)
m = evalMono2(pm, x)
fmt.Printf("p(%4.2f) = %0.14g (mono %0.14g)\n", x, y, m)
x = 7.5
y = evalBern2(pb2, x)
m = evalMono2(pm, x)
fmt.Printf("p(%4.2f) = %0.14g (mono %0.14g)\n", x, y, m)
x = 0.25
y = evalBern2(qb2, x)
m = evalMono2(qm, x)
fmt.Printf("q(%4.2f) = %0.14g (mono %0.14g)\n", x, y, m)
x = 7.5
y = evalBern2(qb2, x)
m = evalMono2(qm, x)
fmt.Printf("q(%4.2f) = %0.14g (mono %0.14g)\n", x, y, m)
fmt.Println("\nSubprogram(3) examples:")
pm = append(pm, 0)
qm = append(qm, 0)
pb3 := toBern3(pm)
qb3 := toBern3(qm)
rb3 := toBern3(rm)
f := "mono $%v --> bern %0.14v\n"
fmt.Printf(f, pm, pb3)
fmt.Printf(f, qm, qb3)
fmt.Printf(f, rm, rb3)
fmt.Println("\nSubprogram(4) examples:")
x = 0.25
y = evalBern3(pb3, x)
m = evalMono3(pm, x)
fmt.Printf("p(%4.2f) = %0.14g (mono %0.14g)\n", x, y, m)
x = 7.5
y = evalBern3(pb3, x)
m = evalMono3(pm, x)
fmt.Printf("p(%4.2f) = %0.14g (mono %0.14g)\n", x, y, m)
x = 0.25
y = evalBern3(qb3, x)
m = evalMono3(qm, x)
fmt.Printf("q(%4.2f) = %0.14g (mono %0.14g)\n", x, y, m)
x = 7.5
y = evalBern3(qb3, x)
m = evalMono3(qm, x)
fmt.Printf("q(%4.2f) = %0.14g (mono %0.14g)\n", x, y, m)
x = 0.25
y = evalBern3(rb3, x)
m = evalMono3(rm, x)
fmt.Printf("r(%4.2f) = %0.14g (mono %0.14g)\n", x, y, m)
x = 7.5
y = evalBern3(rb3, x)
m = evalMono3(rm, x)
fmt.Printf("r(%4.2f) = %0.14g (mono %0.14g)\n", x, y, m)
fmt.Println("\nSubprogram(5) examples:")
pc := bern2to3(pb2)
qc := bern2to3(qb2)
fmt.Printf("mono %v --> bern %0.14v\n", pb2, pc)
fmt.Printf("mono %v --> bern %0.14v\n", qb2, qc)
}</syntaxhighlight>
{{out}}
<pre>
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)
q(0.25) = 1.6875 (mono 1.6875)
q(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.6666666666667 3.3333333333333 6]
mono $[1 2 3 4] --> bern [1 1.6666666666667 3.3333333333333 10]
Subprogram(4) examples:
p(0.25) = 1 (mono 1)
p(7.50) = 1 (mono 1)
q(0.25) = 1.6875 (mono 1.6875)
q(7.50) = 184.75 (mono 184.75)
r(0.25) = 1.75 (mono 1.75)
r(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.6666666666667 3.3333333333333 6]
</pre>
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.util.Arrays;
import java.util.List;
public final class BernsteinBasisPolynomials {
public static void main(String[] args) {
/**
* For the following polynomials, use Subprogram (1) to find coefficients in the degree-2 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
*/
QuadraticCoefficients pMonomial2 = new QuadraticCoefficients(1.0, 0.0, 0.0);
QuadraticCoefficients qMonomial2 = new QuadraticCoefficients(1.0, 2.0, 3.0);
QuadraticCoefficients pBernstein2 = monomialToBernsteinDegree2(pMonomial2);
QuadraticCoefficients qBernstein2 = monomialToBernsteinDegree2(qMonomial2);
System.out.println("Subprogram (1) examples:");
System.out.println(" monomial " + pMonomial2 + " --> bernstein " + pBernstein2);
System.out.println(" monomial " + qMonomial2 + " --> bernstein " + qBernstein2);
/**
* Use Subprogram (2) to evaluate p(x) and q(x) at x = 0.25, 7.50. Display the results.
* Optionally also display results from evaluating in the original monomial basis.
*/
System.out.println("Subprogram (2) examples:");
for ( double x : List.of( 0.25, 7.50 ) ) {
System.out.println(" p(" + x + ") = " + evaluateBernsteinDegree2(pBernstein2, x) +
" ( mono: " + evaluateMonomialDegree2(pMonomial2, x) + " )");
}
for ( double x : List.of( 0.25, 7.50 ) ) {
System.out.println(" q(" + x + ") = " + evaluateBernsteinDegree2(qBernstein2, x) +
" ( mono: " + evaluateMonomialDegree2(qMonomial2, x) + " )");
}
/**
* For the following polynomials, use Subprogram (3) to find coefficients in the degree-3 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
* r(x) = 1 + 2x + 3x² + 4x³
*
* Display the results.
*/
CubicCoefficients pMonomial3 = new CubicCoefficients(1.0, 0.0, 0.0, 0.0);
CubicCoefficients qMonomial3 = new CubicCoefficients(1.0, 2.0, 3.0, 0.0);
CubicCoefficients rMonomial3 = new CubicCoefficients(1.0, 2.0, 3.0, 4.0);
CubicCoefficients pBernstein3 = monomialToBernsteinDegree3(pMonomial3);
CubicCoefficients qBernstein3 = monomialToBernsteinDegree3(qMonomial3);
CubicCoefficients rBernstein3 = monomialToBernsteinDegree3(rMonomial3);
System.out.println("Subprogram (3) examples:");
System.out.println(" monomial " + pMonomial3 + " --> bernstein " + pBernstein3);
System.out.println(" monomial " + qMonomial3 + " --> bernstein " + qBernstein3);
System.out.println(" monomial " + rMonomial3 + " --> bernstein " + rBernstein3);
/**
* Use Subprogram (4) to evaluate p(x), q(x), and r(x) at x = 0.25, 7.50. Display the results.
* Optionally also display results from evaluating in the original monomial basis.
*/
System.out.println("Subprogram (4) examples:");
for ( double x : List.of( 0.25, 7.50 ) ) {
System.out.println(" p(" + x + ") = " + evaluateBernsteinDegree3(pBernstein3, x) +
" ( mono: " + evaluateMonomialDegree3(pMonomial3, x) + " )");
}
for ( double x : List.of( 0.25, 7.50 ) ) {
System.out.println(" q(" + x + ") = " + evaluateBernsteinDegree3(qBernstein3, x) +
" ( mono: " + evaluateMonomialDegree3(qMonomial3, x) + " )");
}
for ( double x : List.of( 0.25, 7.50 ) ) {
System.out.println(" r(" + x + ") = " + evaluateBernsteinDegree3(rBernstein3, x) +
" ( mono: " + evaluateMonomialDegree3(rMonomial3, x) + " )");
}
/**
* For the following polynomials, using the result of Subprogram (1) applied to the polynomial,
* use Subprogram (5) to get coefficients for the degree-3 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
*
* Display the results.
*/
System.out.println("Subprogram (5) examples:");
CubicCoefficients pBernstein3a = bernsteinDegree2ToDegree3(pBernstein2);
CubicCoefficients qBernstein3a = bernsteinDegree2ToDegree3(qBernstein2);
System.out.println(" bernstein " + pBernstein2 + " --> bernstein " + pBernstein3a);
System.out.println(" bernstein " + qBernstein2 + " --> bernstein " + qBernstein3a);
}
private static record QuadraticCoefficients(double q0, double q1, double q2) {
@Override
public String toString() {
return Arrays.asList(q0, q1, q2).toString();
}
}
private static record CubicCoefficients(double c0, double c1, double c2, double c3) {
@Override
public String toString() {
return Arrays.asList(c0, c1, c2, c3).toString();
}
}
// Subprogram (1)
private static QuadraticCoefficients monomialToBernsteinDegree2(QuadraticCoefficients monomial) {
return new QuadraticCoefficients( monomial.q0,
monomial.q0 + ( monomial.q1 / 2.0 ),
monomial.q0 + monomial.q1 + monomial.q2 );
}
// Subprogram (2)
private static double evaluateBernsteinDegree2(QuadraticCoefficients bernstein, double t) {
// de Casteljau’s algorithm
final double s = 1 - t;
final double b01 = ( s * bernstein.q0 ) + ( t * bernstein.q1 );
final double b12 = ( s * bernstein.q1 ) + ( t * bernstein.q2 );
return ( s * b01 ) + ( t * b12 );
}
// Subprogram (3)
private static CubicCoefficients monomialToBernsteinDegree3(CubicCoefficients monomial) {
return new CubicCoefficients( monomial.c0,
monomial.c0 + ( monomial.c1 / 3.0 ),
monomial.c0 + ( 2.0 * monomial.c1 / 3.0 ) + ( monomial.c2 / 3.0 ),
monomial.c0 + monomial.c1 + monomial.c2 + monomial.c3 );
}
// Subprogram (4)
private static double evaluateBernsteinDegree3(CubicCoefficients bernstein, double t) {
// de Casteljau’s algorithm
final double s = 1 - t;
final double b01 = ( s * bernstein.c0 ) + ( t * bernstein.c1 );
final double b12 = ( s * bernstein.c1 ) + ( t * bernstein.c2 );
final double b23 = ( s * bernstein.c2 ) + ( t * bernstein.c3 );
final double b012 = ( s * b01 ) + ( t * b12 );
final double b123 = ( s * b12 ) + ( t * b23 );
return ( s * b012 ) + ( t * b123 );
}
// Subprogram (5)
private static CubicCoefficients bernsteinDegree2ToDegree3(QuadraticCoefficients bernstein) {
return new CubicCoefficients( bernstein.q0,
( bernstein.q0 / 3.0 ) + ( 2.0 * bernstein.q1 / 3.0 ),
( 2.0 * bernstein.q1 / 3.0 ) + ( bernstein.q2 / 3.0 ),
bernstein.q2 );
}
private static double evaluateMonomialDegree2(QuadraticCoefficients monomial, double t) {
// Horner’s rule
return monomial.q0 + ( t * ( monomial.q1 + ( t * monomial.q2 ) ) );
}
private static double evaluateMonomialDegree3(CubicCoefficients monomial, double t) {
// Horner’s rule
return monomial.c0 + ( t * ( monomial.c1 + ( t * ( monomial.c2 + ( t * monomial.c3 ) ) ) ) );
}
}
</syntaxhighlight>
{{ out }}
<pre>
Subprogram (1) examples:
monomial [1.0, 0.0, 0.0] --> bernstein [1.0, 1.0, 1.0]
monomial [1.0, 2.0, 3.0] --> bernstein [1.0, 2.0, 6.0]
Subprogram (2) examples:
p(0.25) = 1.0 ( mono: 1.0 )
p(7.5) = 1.0 ( mono: 1.0 )
q(0.25) = 1.6875 ( mono: 1.6875 )
q(7.5) = 184.75 ( mono: 184.75 )
Subprogram (3) examples:
monomial [1.0, 0.0, 0.0, 0.0] --> bernstein [1.0, 1.0, 1.0, 1.0]
monomial [1.0, 2.0, 3.0, 0.0] --> bernstein [1.0, 1.6666666666666665, 3.333333333333333, 6.0]
monomial [1.0, 2.0, 3.0, 4.0] --> bernstein [1.0, 1.6666666666666665, 3.333333333333333, 10.0]
Subprogram (4) examples:
p(0.25) = 1.0 ( mono: 1.0 )
p(7.5) = 1.0 ( mono: 1.0 )
q(0.25) = 1.6874999999999998 ( mono: 1.6875 )
q(7.5) = 184.75000000000034 ( mono: 184.75 )
r(0.25) = 1.7499999999999998 ( mono: 1.75 )
r(7.5) = 1872.25 ( mono: 1872.25 )
Subprogram (5) examples:
bernstein [1.0, 1.0, 1.0] --> bernstein [1.0, 1.0, 1.0, 1.0]
bernstein [1.0, 2.0, 6.0] --> bernstein [1.0, 1.6666666666666665, 3.333333333333333, 6.0]
</pre>
=={{header|jq}}==
Line 1,492 ⟶ 2,340:
(1.0, 1.0, 1.0) --> (1.0, 1.0, 1.0, 1.0)
(1.0, 2.0, 6.0) --> (1.0, 1.6667, 3.3332, 6.0)</pre>
=={{header|Nim}}==
{{trans|Python}}
<syntaxhighlight lang="Nim">type
Coeffs2 = (float, float, float)
Coeffs3 = (float, float, float, float)
# Subprogram (1).
func monomialToBernsteinDegree2(monomialCoefficients: Coeffs2): Coeffs2 =
let (a0, a1, a2) = monomialCoefficients
result = (a0, a0 + ((1/2) * a1), a0 + a1 + a2)
# Subprogram (2).
func evaluateBernsteinDegree2(bernsteinCoefficients: Coeffs2; t: float): float =
let (b0, b1, b2) = bernsteinCoefficients
# de Casteljau’s algorithm.
let s = 1 - t
let b01 = (s * b0) + (t * b1)
let b12 = (s * b1) + (t * b2)
result = (s * b01) + (t * b12)
# Subprogram (3).
func monomialToBernsteinDegree3(monomialCoefficients: Coeffs3): Coeffs3 =
let (a0, a1, a2, a3) = monomialCoefficients
result = (a0, a0 + ((1/3) * a1), a0 + ((2/3) * a1) + ((1/3) * a2), a0 + a1 + a2 + a3)
# Subprogram (4).
func evaluateBernsteinDegree3(bernsteinCoefficients: Coeffs3; t: float): float =
let (b0, b1, b2, b3) = bernsteinCoefficients
# de Casteljau’s algorithm.
let s = 1 - t
let b01 = (s * b0) + (t * b1)
let b12 = (s * b1) + (t * b2)
let b23 = (s * b2) + (t * b3)
let b012 = (s * b01) + (t * b12)
let b123 = (s * b12) + (t * b23)
result = (s * b012) + (t * b123)
# Subprogram (5).
func bernsteinDegree2ToDegree3(bernsteinCoefficients: Coeffs2): Coeffs3 =
let (b0, b1, b2) = bernstein_coefficients
result = (b0, ((1/3) * b0) + ((2/3) * b1), ((2/3) * b1) + ((1/3) * b2), b2)
func evaluateMonomialDegree2(monomialCoefficients: Coeffs2; t: float): float =
let (a0, a1, a2) = monomialCoefficients
# Horner’s rule.
result = a0 + (t * (a1 + (t * a2)))
func evaluateMonomialDegree3(monomialCoefficients: Coeffs3; t: float): float =
let (a0, a1, a2, a3) = monomialCoefficients
# Horner’s rule.
result = a0 + (t * (a1 + (t * (a2 + (t * a3)))))
#
# For the following polynomials, use Subprogram (1) to find
# coefficients in the degree-2 Bernstein basis:
#
# p(x) = 1
# q(x) = 1 + 2x + 3x²
#
# Display the results.
#
let pmono2 = (1.0, 0.0, 0.0)
let qmono2 = (1.0, 2.0, 3.0)
let pbern2 = monomialToBernsteinDegree2(pmono2)
let qbern2 = monomialToBernsteinDegree2(qmono2)
echo "Subprogram (1) examples:"
echo " mono ", pmono2, " --> bern ", pbern2
echo " mono ", qmono2, " --> bern ", qbern2
#
# Use Subprogram (2) to evaluate p(x) and q(x) at x = 0.25, 7.50.
# Display the results. Optionally also display results from evaluating
# in the original monomial basis.
#
echo "Subprogram (2) examples:"
for x in [0.25, 7.50]:
echo " p(", x, ") = ", evaluateBernsteinDegree2(pbern2, x),
" (mono: ", evaluateMonomialDegree2(pmono2, x), ')'
for x in [0.25, 7.50]:
echo " q(", x, ") = ", evaluateBernsteinDegree2(qbern2, x),
" (mono: ", evaluateMonomialDegree2(qmono2, x), ')'
#
# For the following polynomials, use Subprogram (3) to find
# coefficients in the degree-3 Bernstein basis:
#
# p(x) = 1
# q(x) = 1 + 2x + 3x²
# r(x) = 1 + 2x + 3x² + 4x³
#
# Display the results.
#
let pmono3 = (1.0, 0.0, 0.0, 0.0)
let qmono3 = (1.0, 2.0, 3.0, 0.0)
let rmono3 = (1.0, 2.0, 3.0, 4.0)
let pbern3 = monomialToBernsteinDegree3(pmono3)
let qbern3 = monomialToBernsteinDegree3(qmono3)
let rbern3 = monomialToBernsteinDegree3(rmono3)
echo "Subprogram (3) examples:"
echo " mono ", pmono3, " --> bern ", pbern3
echo " mono ", qmono3, " --> bern ", qbern3
echo " mono ", rmono3, " --> bern ", rbern3
#
# Use Subprogram (4) to evaluate p(x), q(x), and r(x) at x = 0.25,
# 7.50. Display the results. Optionally also display results from
# evaluating in the original monomial basis.
#
echo "Subprogram (4) examples:"
for x in [0.25, 7.50]:
echo " p(", x, ") = ", evaluateBernsteinDegree3(pbern3, x),
" (mono: ", evaluateMonomialDegree3(pmono3, x), ')'
for x in [0.25, 7.50]:
echo " q(", x, ") = ", evaluateBernsteinDegree3(qbern3, x),
" (mono: ", evaluateMonomialDegree3(qmono3, x), ')'
for x in [0.25, 7.50]:
echo " r(", x, ") = ", evaluateBernsteinDegree3(rbern3, x),
" (mono: ", evaluateMonomialDegree3(rmono3, x), ')'
#
# For the following polynomials, using the result of Subprogram (1)
# applied to the polynomial, use Subprogram (5) to get coefficients
# for the degree-3 Bernstein basis:
#
# p(x) = 1
# q(x) = 1 + 2x + 3x²
#
# Display the results.
#
echo "Subprogram (5) examples:"
let pbern3a = bernsteinDegree2ToDegree3(pbern2)
let qbern3a = bernsteinDegree2ToDegree3(qbern2)
echo " bern ", pbern2, "--> bern ", pbern3a
echo " bern ", qbern2, "--> bern ", qbern3a
</syntaxhighlight>
{{out}}
<pre>Subprogram (1) examples:
mono (1.0, 0.0, 0.0) --> bern (1.0, 1.0, 1.0)
mono (1.0, 2.0, 3.0) --> bern (1.0, 2.0, 6.0)
Subprogram (2) examples:
p(0.25) = 1.0 (mono: 1.0)
p(7.5) = 1.0 (mono: 1.0)
q(0.25) = 1.6875 (mono: 1.6875)
q(7.5) = 184.75 (mono: 184.75)
Subprogram (3) examples:
mono (1.0, 0.0, 0.0, 0.0) --> bern (1.0, 1.0, 1.0, 1.0)
mono (1.0, 2.0, 3.0, 0.0) --> bern (1.0, 1.666666666666667, 3.333333333333333, 6.0)
mono (1.0, 2.0, 3.0, 4.0) --> bern (1.0, 1.666666666666667, 3.333333333333333, 10.0)
Subprogram (4) examples:
p(0.25) = 1.0 (mono: 1.0)
p(7.5) = 1.0 (mono: 1.0)
q(0.25) = 1.6875 (mono: 1.6875)
q(7.5) = 184.7500000000003 (mono: 184.75)
r(0.25) = 1.75 (mono: 1.75)
r(7.5) = 1872.25 (mono: 1872.25)
Subprogram (5) examples:
bern (1.0, 1.0, 1.0)--> bern (1.0, 1.0, 1.0, 1.0)
bern (1.0, 2.0, 6.0)--> bern (1.0, 1.666666666666667, 3.333333333333333, 6.0)
</pre>
=={{header|ObjectIcon}}==
Line 2,003 ⟶ 3,012:
bern (1, 1.0, 1) --> bern (1, 1.0, 1.0, 1)
bern (1, 2.0, 6) --> bern (1, 1.6666666666666665, 3.333333333333333, 6)</pre>
=={{header|Raku}}==
{{trans|Wren}}
<syntaxhighlight lang="raku" line># 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;
</syntaxhighlight>
{{out}}
<pre>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]</pre>
=={{header|Wren}}==
Line 2,009 ⟶ 3,099:
{{libheader|Wren-math}}
Note that the library method, ''Math.evalPoly'', evaluates polynomials of any degree using Horner's rule but requires the coefficients to be presented in highest to lowest degree order.
<syntaxhighlight lang="
import "./math" for Math
| |||