Bernstein basis polynomials
The Bernstein basis polynomials of degree are defined as
Any real polynomial can written as a linear combination of such Bernstein basis polynomials. Let us call the coefficients in such linear combinations Bernstein coefficients.
The goal of this task is to write subprograms for working with degree-2 and degree-3 Bernstein coefficients. A programmer is likely to have to deal with such representations. For example, OpenType fonts store glyph outline data as as either degree-2 or degree-3 Bernstein coefficients.
The task is as follows:
- Write a subprogram (which we will call Subprogram (1)) to transform the ordinary monomial coefficients of a real polynomial, of degree 2 or less, to Bernstein coefficients for the degree-2 basis.
- Write Subprogram (2) to evaluate, at a given point, a real polynomial written in degree-2 Bernstein coefficients. For this use de Casteljau's algorithm.
- Write Subprogram (3) to transform the monomial coefficients of a real polynomial, of degree 3 or less, to Bernstein coefficients for degree 3.
- Write Subprogram (4) to evaluate, at a given point, a real polynomial written in degree-3 Bernstein coefficients. Use de Casteljau's algorithm.
- Write Subprogram (5) to transform Bernstein coefficients for degree 2 to Bernstein coefficients for degree 3.
- For the following polynomials, use Subprogram (1) to find Bernstein coefficients for the degree-2 basis: , . Display the results.
- Use Subprogram (2) to evaluate and at and at . Display the results. Optionally also evaluate the polynomials in the monomial basis and display the results. (For the monomial basis it is best to use a Horner scheme.)
- Use Subprogram (3) to find degree-3 Bernstein coefficients for and , and additionally also for . Display the results.
- Use Subprogram (4) to evaluate , , and at and . Display the results. Optionally also evaluate them in the monomial basis and display the results. (You may have done that already for and .)
- Use Subprogram (5) to get degree-3 Bernstein coefficients for and from their degree-2 Bernstein coefficients. Display the results.
ALGOL 60 and Python implementations are provided as initial examples. The latter does the optional monomial-basis evaluations.
Here is ALGOL 60 more nicely printed:
procedure tobern2 (a0, a1, a2, b0, b1, b2); value a0, a1, a2; comment pass by value; real a0, a1, a2; real b0, b1, b2; comment Subprogram (1): transform monomial coefficients a0, a1, a2 to Bernstein coefficients b0, b1, b2; begin b0 := a0; b1 := a0 + ((1/2) × a1); b2 := a0 + a1 + a2 end tobern2;
real procedure evalbern2 (b0, b1, b2, t); value b0, b1, b2, t; real b0, b1, b2, t; comment Subprogram (2): evaluate, at t, the polynomial with Bernstein coefficients b0, b1, b2. Use de Casteljau’s algorithm; begin real s, b01, b12, b012; s := 1 - t; b01 := (s × b0) + (t × b1); b12 := (s × b1) + (t × b2); b012 := (s × b01) + (t × b12); evalbern2 := b012 end evalbern2;
procedure tobern3 (a0, a1, a2, a3, b0, b1, b2, b3); value a0, a1, a2, a3; real a0, a1, a2, a3; real b0, b1, b2, b3; comment Subprogram (3): transform monomial coefficients a0, a1, a2, a3 to Bernstein coefficients b0, b1, b2, b3; begin b0 := a0; b1 := a0 + ((1/3) × a1); b2 := a0 + ((2/3) × a1) + ((1/3) × a2); b3 := a0 + a1 + a2 + a3 end tobern3;
real procedure evalbern3 (b0, b1, b2, b3, t); value b0, b1, b2, b3, t; real b0, b1, b2, b3, t; comment Subprogram (4): evaluate, at t, the polynomial with Bernstein coefficients b0, b1, b2, b3. Use de Casteljau's algorithm; begin real s, b01, b12, b23, b012, b123, b0123; s := 1 - t; b01 := (s × b0) + (t × b1); b12 := (s × b1) + (t × b2); b23 := (s × b2) + (t × b3); b012 := (s × b01) + (t × b12); b123 := (s × b12) + (t × b23); b0123 := (s × b012) + (t × b123); evalbern3 := b0123 end evalbern3;
procedure bern2to3 (q0, q1, q2, c0, c1, c2, c3); value q0, q1, q2; real q0, q1, q2; real c0, c1, c2, c3; comment Subprogram (5): transform the quadratic Bernstein coefficients q0, q1, q2 to the cubic Bernstein coefficients c0, c1, c2, c3; begin c0 := q0; c1 := ((1/3) × q0) + ((2/3) × q1); c2 := ((2/3) × q1) + ((1/3) × q2); c3 := q2 end bern2to3;
ALGOL 60
procedure tobern2 (a0, a1, a2, b0, b1, b2);
value a0, a1, a2; comment pass by value;
real a0, a1, a2;
real b0, b1, b2;
comment Subprogram (1): transform monomial coefficients
a0, a1, a2 to Bernstein coefficients b0, b1, b2;
begin
b0 := a0;
b1 := a0 + ((1/2) * a1);
b2 := a0 + a1 + a2
end tobern2;
real procedure evalbern2 (b0, b1, b2, t);
value b0, b1, b2, t;
real b0, b1, b2, t;
comment Subprogram (2): evaluate, at t, the polynomial with
Bernstein coefficients b0, b1, b2. Use de Casteljau's
algorithm;
begin
real s, b01, b12, b012;
s := 1 - t;
b01 := (s * b0) + (t * b1);
b12 := (s * b1) + (t * b2);
b012 := (s * b01) + (t * b12);
evalbern2 := b012
end evalbern2;
procedure tobern3 (a0, a1, a2, a3, b0, b1, b2, b3);
value a0, a1, a2, a3;
real a0, a1, a2, a3;
real b0, b1, b2, b3;
comment Subprogram (3): transform monomial coefficients
a0, a1, a2, a3 to Bernstein coefficients b0, b1,
b2, b3;
begin
b0 := a0;
b1 := a0 + ((1/3) * a1);
b2 := a0 + ((2/3) * a1) + ((1/3) * a2);
b3 := a0 + a1 + a2 + a3
end tobern3;
real procedure evalbern3 (b0, b1, b2, b3, t);
value b0, b1, b2, b3, t;
real b0, b1, b2, b3, t;
comment Subprogram (4): evaluate, at t, the polynomial
with Bernstein coefficients b0, b1, b2, b3. Use
de Casteljau's algorithm;
begin
real s, b01, b12, b23, b012, b123, b0123;
s := 1 - t;
b01 := (s * b0) + (t * b1);
b12 := (s * b1) + (t * b2);
b23 := (s * b2) + (t * b3);
b012 := (s * b01) + (t * b12);
b123 := (s * b12) + (t * b23);
b0123 := (s * b012) + (t * b123);
evalbern3 := b0123
end evalbern3;
procedure bern2to3 (q0, q1, q2, c0, c1, c2, c3);
value q0, q1, q2;
real q0, q1, q2;
real c0, c1, c2, c3;
comment Subprogram (5): transform the quadratic Bernstein
coefficients q0, q1, q2 to the cubic Bernstein
coefficients c0, c1, c2, c3;
begin
c0 := q0;
c1 := ((1/3) * q0) + ((2/3) * q1);
c2 := ((2/3) * q1) + ((1/3) * q2);
c3 := q2
end bern2to3;
begin
real p0m, p1m, p2m;
real q0m, q1m, q2m;
real r0m, r1m, r2m, r3m;
real p0b2, p1b2, p2b2;
real q0b2, q1b2, q2b2;
real p0b3, p1b3, p2b3, p3b3;
real q0b3, q1b3, q2b3, q3b3;
real r0b3, r1b3, r2b3, r3b3;
real pc0, pc1, pc2, pc3;
real qc0, qc1, qc2, qc3;
real x, y;
p0m := 1; p1m := 0; p2m := 0;
q0m := 1; q1m := 2; q2m := 3;
r0m := 1; r1m := 2; r2m := 3; r3m := 4;
tobern2 (p0m, p1m, p2m, p0b2, p1b2, p2b2);
tobern2 (q0m, q1m, q2m, q0b2, q1b2, q2b2);
outstring (1, "Subprogram (1) examples:\n");
outstring (1, " mono ( ");
outreal (1, p0m); outstring (1, ", ");
outreal (1, p1m); outstring (1, ", ");
outreal (1, p2m); outstring (1, ") --> bern ( ");
outreal (1, p0b2); outstring (1, ", ");
outreal (1, p1b2); outstring (1, ", ");
outreal (1, p2b2); outstring (1, ")\n");
outstring (1, " mono ( ");
outreal (1, q0m); outstring (1, ", ");
outreal (1, q1m); outstring (1, ", ");
outreal (1, q2m); outstring (1, ") --> bern ( ");
outreal (1, q0b2); outstring (1, ", ");
outreal (1, q1b2); outstring (1, ", ");
outreal (1, q2b2); outstring (1, ")\n");
outstring (1, "Subprogram (2) examples:\n");
x := 0.25;
y := evalbern2 (p0b2, p1b2, p2b2, x);
outstring (1, " p ( "); outreal (1, x);
outstring (1, ") = "); outreal (1, y);
outstring (1, "\n");
x := 7.50;
y := evalbern2 (p0b2, p1b2, p2b2, x);
outstring (1, " p ( "); outreal (1, x);
outstring (1, ") = "); outreal (1, y);
outstring (1, "\n");
x := 0.25;
y := evalbern2 (q0b2, q1b2, q2b2, x);
outstring (1, " q ( "); outreal (1, x);
outstring (1, ") = "); outreal (1, y);
outstring (1, "\n");
x := 7.50;
y := evalbern2 (q0b2, q1b2, q2b2, x);
outstring (1, " q ( "); outreal (1, x);
outstring (1, ") = "); outreal (1, y);
outstring (1, "\n");
tobern3 (p0m, p1m, p2m, 0, p0b3, p1b3, p2b3, p3b3);
tobern3 (q0m, q1m, q2m, 0, q0b3, q1b3, q2b3, q3b3);
tobern3 (r0m, r1m, r2m, r3m, r0b3, r1b3, r2b3, r3b3);
outstring (1, "Subprogram (3) examples:\n");
outstring (1, " mono ( ");
outreal (1, p0m); outstring (1, ", ");
outreal (1, p1m); outstring (1, ", ");
outreal (1, p2m); outstring (1, ", ");
outreal (1, 0); outstring (1, ") --> bern ( ");
outreal (1, p0b3); outstring (1, ", ");
outreal (1, p1b3); outstring (1, ", ");
outreal (1, p2b3); outstring (1, ", ");
outreal (1, p3b3); outstring (1, ")\n");
outstring (1, " mono ( ");
outreal (1, q0m); outstring (1, ", ");
outreal (1, q1m); outstring (1, ", ");
outreal (1, q2m); outstring (1, ", ");
outreal (1, 0); outstring (1, ") --> bern ( ");
outreal (1, q0b3); outstring (1, ", ");
outreal (1, q1b3); outstring (1, ", ");
outreal (1, q2b3); outstring (1, ", ");
outreal (1, q3b3); outstring (1, ")\n");
outstring (1, " mono ( ");
outreal (1, r0m); outstring (1, ", ");
outreal (1, r1m); outstring (1, ", ");
outreal (1, r2m); outstring (1, ", ");
outreal (1, r3m); outstring (1, ") --> bern ( ");
outreal (1, r0b3); outstring (1, ", ");
outreal (1, r1b3); outstring (1, ", ");
outreal (1, r2b3); outstring (1, ", ");
outreal (1, r3b3); outstring (1, ")\n");
outstring (1, "Subprogram (4) examples:\n");
x := 0.25;
y := evalbern3 (p0b3, p1b3, p2b3, p3b3, x);
outstring (1, " p ( "); outreal (1, x);
outstring (1, ") = "); outreal (1, y);
outstring (1, "\n");
x := 7.50;
y := evalbern3 (p0b3, p1b3, p2b3, p3b3, x);
outstring (1, " p ( "); outreal (1, x);
outstring (1, ") = "); outreal (1, y);
outstring (1, "\n");
x := 0.25;
y := evalbern3 (q0b3, q1b3, q2b3, q3b3, x);
outstring (1, " q ( "); outreal (1, x);
outstring (1, ") = "); outreal (1, y);
outstring (1, "\n");
x := 7.50;
y := evalbern3 (q0b3, q1b3, q2b3, q3b3, x);
outstring (1, " q ( "); outreal (1, x);
outstring (1, ") = "); outreal (1, y);
outstring (1, "\n");
x := 0.25;
y := evalbern3 (r0b3, r1b3, r2b3, r3b3, x);
outstring (1, " r ( "); outreal (1, x);
outstring (1, ") = "); outreal (1, y);
outstring (1, "\n");
x := 7.50;
y := evalbern3 (r0b3, r1b3, r2b3, r3b3, x);
outstring (1, " r ( "); outreal (1, x);
outstring (1, ") = "); outreal (1, y);
outstring (1, "\n");
bern2to3 (p0b2, p1b2, p2b2, pc0, pc1, pc2, pc3);
bern2to3 (q0b2, q1b2, q2b2, qc0, qc1, qc2, qc3);
outstring (1, "Subprogram (5) examples:\n");
outstring (1, " bern ( ");
outreal (1, p0b2); outstring (1, ", ");
outreal (1, p1b2); outstring (1, ", ");
outreal (1, p2b2); outstring (1, ") --> bern ( ");
outreal (1, pc0); outstring (1, ", ");
outreal (1, pc1); outstring (1, ", ");
outreal (1, pc2); outstring (1, ", ");
outreal (1, pc3); outstring (1, ")\n");
outstring (1, " bern ( ");
outreal (1, q0b2); outstring (1, ", ");
outreal (1, q1b2); outstring (1, ", ");
outreal (1, q2b2); outstring (1, ") --> bern ( ");
outreal (1, qc0); outstring (1, ", ");
outreal (1, qc1); outstring (1, ", ");
outreal (1, qc2); outstring (1, ", ");
outreal (1, qc3); outstring (1, ")\n")
end
- 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 p ( 7.5 ) = 1 q ( 0.25 ) = 1.6875 q ( 7.5 ) = 184.75 Subprogram (3) examples: mono ( 1 , 0 , 0 , 0 ) --> bern ( 1 , 1 , 1 , 1 ) mono ( 1 , 2 , 3 , 0 ) --> bern ( 1 , 1.66666666667 , 3.33333333333 , 6 ) mono ( 1 , 2 , 3 , 4 ) --> bern ( 1 , 1.66666666667 , 3.33333333333 , 10 ) Subprogram (4) examples: p ( 0.25 ) = 1 p ( 7.5 ) = 1 q ( 0.25 ) = 1.6875 q ( 7.5 ) = 184.75 r ( 0.25 ) = 1.75 r ( 7.5 ) = 1872.25 Subprogram (5) examples: bern ( 1 , 1 , 1 ) --> bern ( 1 , 1 , 1 , 1 ) bern ( 1 , 2 , 6 ) --> bern ( 1 , 1.66666666667 , 3.33333333333 , 6 )
jq
Adapted from Wren
Works with jq and gojq, the C and Go implementations of jq
def toBern2(a):
[a[0], a[0] + a[1] / 2,
a[0] + a[1] + a[2]];
# uses de Casteljau's algorithm
def evalBern2($b; $t):
(1 - $t) as $s
| ($s * $b[0] + $t * $b[1]) as $b01
| ($s * $b[1] + $t * $b[2]) as $b12
| $s * $b01 + $t * $b12;
def toBern3(a):
[a[0],
a[0] + a[1] / 3,
a[0] + a[1] * 2/3 + a[2] / 3,
a[0] + a[1] + a[2] + a[3] ];
# uses de Casteljau's algorithm
def evalBern3($b; $t):
(1 - $t) as $s
| ($s * $b[0] + $t * $b[1]) as $b01
| ($s * $b[1] + $t * $b[2]) as $b12
| ($s * $b[2] + $t * $b[3]) as $b23
| ($s * $b01 + $t * $b12) as $b012
| ($s * $b12 + $t * $b23) as $b123
| $s * $b012 + $t * $b123;
def bern2to3(q):
[q[0],
q[0] / 3 + q[1] * 2/3,
q[1] * 2/3 + q[2] / 3,
q[2]] ;
def pm: [1, 0, 0];
def qm: [1, 2, 3];
def rm: [1, 2, 3, 4];
def pb2: toBern2(pm);
def qb2: toBern2(qm);
"Subprogram(1) examples:",
"mono\(pm) --> bern\(pb2)",
"mono\(qm) --> bern\(qb2)",
"\nSubprogram(2) examples:",
({x: 0.25}
| .y = evalBern2(pb2; .x)
| "p(\(.x)) = \(.y)" ),
({x: 7.5}
| .y = evalBern2(pb2; .x)
| "p(\(.x)) = \(.y)" ),
({x: 0.25}
| .y = evalBern2(qb2; .x)
| "q(\(.x)) = \(.y)" ),
({x: 7.5}
| .y = evalBern2(qb2; .x)
| "q(\(.x)) = \(.y)" ),
"\nSubprogram(3) examples:",
({}
| .pm0 = pm + [0]
| .qm0 = qm + [0]
| .pb3 = toBern3(.pm0)
| .qb3 = toBern3(.qm0)
| .rb3 = toBern3(rm)
| "mono\(.pm0) --> bern\(.pb3)",
"mono\(.qm0) --> bern\(.qb3)",
"mono\(rm) --> bern\(.rb3)",
"\nSubprogram(4) examples:",
( .x = 0.25
| .y = evalBern3(.pb3; .x)
| "p(\(.x)) = \(.y)" ),
( .x = 7.5
| .y = evalBern3(.pb3; .x)
| "p(\(.x)) = \(.y)" ),
( .x = 0.25
| .y = evalBern3(.qb3; .x)
| "q(\(.x)) = \(.y)" ),
( .x = 7.5
| .y = evalBern3(.qb3; .x)
| "q(\(.x)) = \(.y)" ),
( .x = 0.25
| .y = evalBern3(.rb3; .x)
| "r(\(.x)) = \(.y)" ),
( .x = 7.5
| .y = evalBern3(.rb3; .x)
| "r(\(.x)) = \(.y)" ),
"\nSubprogram(5) examples:",
"mono\(pb2) --> bern\( bern2to3(pb2))",
"mono\(qb2) --> bern\( bern2to3(qb2) )"
)
- 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 p(7.5) = 1 q(0.25) = 1.6875 q(7.5) = 184.75 Subprogram(3) examples: mono[1,0,0,0] --> bern[1,1,1,1] mono[1,2,3,0] --> bern[1,1.6666666666666665,3.333333333333333,6] mono[1,2,3,4] --> bern[1,1.6666666666666665,3.333333333333333,10] Subprogram(4) examples: p(0.25) = 1 p(7.5) = 1 q(0.25) = 1.6874999999999998 q(7.5) = 184.75000000000034 r(0.25) = 1.7499999999999998 r(7.5) = 1872.25 Subprogram(5) examples: mono[1,1,1] --> bern[1,1,1,1] mono[1,2,6] --> bern[1,1.6666666666666665,3.333333333333333,6]
Julia
""" rosettacode.org/wiki/Bernstein_basis_polynomials """
""" Subprogram (1). """
monomial_to_bernstein_degree2(mc) = ((a0, a1, a2) = mc; (a0, a0 + ((1/2) * a1), a0 + a1 + a2))
""" Subprogram (2). """
function evaluate_bernstein_degree2(bernstein_coefficients, t)
(b0, b1, b2) = bernstein_coefficients
# de Casteljau’s algorithm.
s = 1 - t
b01 = (s * b0) + (t * b1)
b12 = (s * b1) + (t * b2)
b012 = (s * b01) + (t * b12)
return b012
end
""" Subprogram (3). """
function monomial_to_bernstein_degree3(monomial_coefficients)
(a0, a1, a2, a3) = monomial_coefficients
return (a0, a0 + ((1/3) * a1),
a0 + ((2/3) * a1) + ((1/3) * a2),
a0 + a1 + a2 + a3)
end
""" Subprogram (4). """
function evaluate_bernstein_degree3(bernstein_coefficients, t)
(b0, b1, b2, b3) = bernstein_coefficients
# de Casteljau’s algorithm.
s = 1 - t
b01 = (s * b0) + (t * b1)
b12 = (s * b1) + (t * b2)
b23 = (s * b2) + (t * b3)
b012 = (s * b01) + (t * b12)
b123 = (s * b12) + (t * b23)
b0123 = (s * b012) + (t * b123)
return b0123
end
""" Subprogram (5). """
function bernstein_degree2_to_degree3(bernstein_coefficients)
(b0, b1, b2) = bernstein_coefficients
return (b0, ((1/3) * b0) + ((2/3) * b1),
((2/3) * b1) + ((1/3) * b2), b2)
end
evaluate_monomial_degree2(mc, t) = ((a0, a1, a2) = mc; a0 + (t * (a1 + (t * a2))))
evaluate_monomial_degree3(mc, t) = ((a0, a1, a2, a3) = mc; 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.
#
pmono2 = (1, 0, 0)
qmono2 = (1, 2, 3)
pbern2 = monomial_to_bernstein_degree2(pmono2)
qbern2 = monomial_to_bernstein_degree2(qmono2)
println("Subprogram (1) examples:")
println(" mono ", pmono2, " --> bern ", pbern2)
println(" 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.
println("Subprogram (2) examples:")
for x in (0.25, 7.50)
println(" p(", x, ") = ", evaluate_bernstein_degree2(pbern2, x),
" ( mono: ", evaluate_monomial_degree2(pmono2, x), ")")
end
for x in (0.25, 7.50)
println(" q(", x, ") = ", evaluate_bernstein_degree2(qbern2, x),
" ( mono: ", evaluate_monomial_degree2(qmono2, x), ")")
end
#
# 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.
#
pmono3 = (1, 0, 0, 0)
qmono3 = (1, 2, 3, 0)
rmono3 = (1, 2, 3, 4)
pbern3 = monomial_to_bernstein_degree3(pmono3)
qbern3 = monomial_to_bernstein_degree3(qmono3)
rbern3 = monomial_to_bernstein_degree3(rmono3)
println("Subprogram (3) examples:")
println(" mono ", pmono3, " --> bern ", pbern3)
println(" mono ", qmono3, " --> bern ", qbern3)
println(" 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.
println("Subprogram (4) examples:")
for x in (0.25, 7.50)
println(" p(", x, ") = " , evaluate_bernstein_degree3(pbern3, x),
" ( mono: " , evaluate_monomial_degree3(pmono3, x), ")")
end
for x in (0.25, 7.50)
println(" q(", x, ") = ", evaluate_bernstein_degree3(qbern3, x),
" ( mono: ", evaluate_monomial_degree3(qmono3, x), ")")
end
for x in (0.25, 7.50)
println(" r(", x, ") = ", evaluate_bernstein_degree3(rbern3, x),
" ( mono: ", evaluate_monomial_degree3(rmono3, x), ")")
end
#
# 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.
#
println("Subprogram (5) examples:")
pbern3a = bernstein_degree2_to_degree3(pbern2)
qbern3a = bernstein_degree2_to_degree3(qbern2)
println(" bern ", pbern2, " --> bern ", pbern3a)
println(" bern ", qbern2, " --> bern ", qbern3a)
- Output:
Same as Python example.
ObjectIcon
#!/bin/env oiscript
import
io(),
ipl.lists(list2str)
final abstract class Bernstein ()
public static from_monomial_degree2 (a0_a1_a2)
# Subprogram (1): transform monomial coefficients [a0, a1, a2] to
# Bernstein coefficients [b0, b1, b2].
local a0, a1, a2
a0 := a0_a1_a2[1]; a1 := a0_a1_a2[2]; a2 := a0_a1_a2[3]
return [a0, a0 + (0.5 * a1), a0 + a1 + a2]
end
public static evaluate_degree2 (b0_b1_b2, t)
# Subprogram (2): evaluate, at t, the polynomial with Bernstein
# coefficients [b0, b1, b2]. Use de Casteljau's
# algorithm.
local b0, b1, b2
local s, b01, b12, b012
b0 := b0_b1_b2[1]; b1 := b0_b1_b2[2]; b2 := b0_b1_b2[3]
s := 1 - t
b01 := (s * b0) + (t * b1)
b12 := (s * b1) + (t * b2)
b012 := (s * b01) + (t * b12)
return b012
end
public static from_monomial_degree3 (a0_a1_a2_a3)
# Subprogram (3): transform monomial coefficients [a0, a1, a2, a3]
# to Bernstein coefficients [b0, b1, b2, b3].
local a0, a1, a2, a3
a0 := a0_a1_a2_a3[1]; a1 := a0_a1_a2_a3[2]
a2 := a0_a1_a2_a3[3]; a3 := a0_a1_a2_a3[4]
#
# In Object Icon, the same symbol / is used for both floating
# point and integer division. Thus you will get the WRONG result
# if you write 1/3 instead of 1.0/3.0, etc. (This is
# error-encouraging language design but true of many languages. I
# fell victim while writing this code.)
#
return [a0, a0 + ((1.0 / 3.0) * a1),
a0 + ((2.0 / 3.0) * a1) + ((1.0 / 3.0) * a2),
a0 + a1 + a2 + a3]
end
public static evaluate_degree3 (b0_b1_b2_b3, t)
# Subprogram (4): evaluate, at t, the polynomial with Bernstein
# coefficients [b0, b1, b2, b3]. Use de Casteljau's
# algorithm.
local b0, b1, b2, b3
local s, b01, b12, b23, b012, b123, b0123
b0 := b0_b1_b2_b3[1]; b1 := b0_b1_b2_b3[2]
b2 := b0_b1_b2_b3[3]; b3 := b0_b1_b2_b3[4]
s := 1 - t
b01 := (s * b0) + (t * b1)
b12 := (s * b1) + (t * b2)
b23 := (s * b2) + (t * b3)
b012 := (s * b01) + (t * b12)
b123 := (s * b12) + (t * b23)
b0123 := (s * b012) + (t * b123)
return b0123
end
public static degree2_to_degree3 (q0_q1_q2)
# Subprogram (5): transform the quadratic Bernstein coefficients
# [q0, q1, q2] to the cubic Bernstein coefficients
# [c0, c1, c2, c3].
local q0, q1, q2
q0 := q0_q1_q2[1]; q1 := q0_q1_q2[2]; q2 := q0_q1_q2[3]
return [q0, ((1.0 / 3.0) * q0) + ((2.0 / 3.0) * q1),
((2.0 / 3.0) * q1) + ((1.0 / 3.0) * q2), q2]
end
end # final abstract class Bernstein
final abstract class Monomial ()
public static evaluate_degree2 (a0_a1_a2, t)
# Evaluate, at t, the polynomial with monomial coefficients [a0,
# a1, a2].
local a0, a1, a2
a0 := a0_a1_a2[1]; a1 := a0_a1_a2[2]; a2 := a0_a1_a2[3]
return a0 + (t * (a1 + (t * a2))) # Horner form.
end
public static evaluate_degree3 (a0_a1_a2_a3, t)
# Evaluate, at t, the polynomial with monomial coefficients [a0,
# a1, a2, a3].
local a0, a1, a2, a3
a0 := a0_a1_a2_a3[1]; a1 := a0_a1_a2_a3[2]
a2 := a0_a1_a2_a3[3]; a3 := a0_a1_a2_a3[4]
return a0 + (t * (a1 + (t * (a2 + (t * a3))))) # Horner form.
end
end # final abstract class Monomial
procedure lstr (lst)
return "[" || list2str (lst, ", ") || "]"
end
procedure main ()
local x
local pmono2, qmono2, pbern2, qbern2
local pmono3, qmono3, rmono3, pbern3, qbern3, rbern3
local pbern3a, qbern3a
#
# 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.
#
pmono2 := [1, 0, 0]
qmono2 := [1, 2, 3]
pbern2 := Bernstein.from_monomial_degree2 (pmono2)
qbern2 := Bernstein.from_monomial_degree2 (qmono2)
io.write ("Subprogram (1) examples:")
io.write (" mono ", lstr (pmono2), " --> bern ", lstr (pbern2))
io.write (" mono ", lstr (qmono2), " --> bern ", lstr (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.
#
io.write ("Subprogram (2) examples:")
every x := 0.25 | 7.50 do
io.write (" p(", x, ") = ",
Bernstein.evaluate_degree2 (pbern2, x),
" (mono: ",
Monomial.evaluate_degree2 (pmono2, x), ")")
every x := 0.25 | 7.50 do
io.write (" q(", x, ") = ",
Bernstein.evaluate_degree2 (qbern2, x),
" (mono: ",
Monomial.evaluate_degree2 (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.
#
pmono3 := [1, 0, 0, 0]
qmono3 := [1, 2, 3, 0]
rmono3 := [1, 2, 3, 4]
pbern3 := Bernstein.from_monomial_degree3 (pmono3)
qbern3 := Bernstein.from_monomial_degree3 (qmono3)
rbern3 := Bernstein.from_monomial_degree3 (rmono3)
io.write ("Subprogram (3) examples:")
io.write (" mono ", lstr (pmono3), " --> bern ", lstr (pbern3))
io.write (" mono ", lstr (qmono3), " --> bern ", lstr (qbern3))
io.write (" mono ", lstr (rmono3), " --> bern ", lstr (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.
#
io.write ("Subprogram (4) examples:")
every x := 0.25 | 7.50 do
io.write (" p(", x, ") = ",
Bernstein.evaluate_degree3 (pbern3, x),
" (mono: ",
Monomial.evaluate_degree3 (pmono3, x), ")")
every x := 0.25 | 7.50 do
io.write (" q(", x, ") = ",
Bernstein.evaluate_degree3 (qbern3, x),
" (mono: ",
Monomial.evaluate_degree3 (qmono3, x), ")")
every x := 0.25 | 7.50 do
io.write (" r(", x, ") = ",
Bernstein.evaluate_degree3 (rbern3, x),
" (mono: ",
Monomial.evaluate_degree3 (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.
#
io.write ("Subprogram (5) examples:")
pbern3a := Bernstein.degree2_to_degree3 (pbern2)
qbern3a := Bernstein.degree2_to_degree3 (qbern2)
io.write (" bern ", lstr (pbern2), " --> bern ", lstr (pbern3a))
io.write (" bern ", lstr (qbern2), " --> bern ", lstr (qbern3a))
end
- Output:
Subprogram (1) examples: mono [1, 0, 0] --> bern [1, 1.0, 1] mono [1, 2, 3] --> bern [1, 2.0, 6] 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] --> bern [1, 1.0, 1.0, 1] mono [1, 2, 3, 0] --> bern [1, 1.666666667, 3.333333333, 6] mono [1, 2, 3, 4] --> bern [1, 1.666666667, 3.333333333, 10] 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.75 (mono: 184.75) r(0.25) = 1.75 (mono: 1.75) r(7.5) = 1872.25 (mono: 1872.25) Subprogram (5) examples: bern [1, 1.0, 1] --> bern [1, 1.0, 1.0, 1] bern [1, 2.0, 6] --> bern [1, 1.666666667, 3.333333333, 6]
Phix
with javascript_semantics function monomial_to_bernstein_degree2(sequence monomial_coefficients) atom {a0, a1, a2} = monomial_coefficients return {a0, a0 + ((1/2) * a1), a0 + a1 + a2} end function function evaluate_bernstein_degree2(sequence bernstein_coefficients, atom t) atom {b0, b1, b2} = bernstein_coefficients -- de Casteljau’s algorithm. atom s = 1 - t, b01 = (s * b0) + (t * b1), b12 = (s * b1) + (t * b2), b012 = (s * b01) + (t * b12) return b012 end function function monomial_to_bernstein_degree3(sequence monomial_coefficients) atom {a0, a1, a2, a3} = monomial_coefficients return {a0, a0 + ((1/3) * a1), a0 + ((2/3) * a1) + ((1/3) * a2), a0 + a1 + a2 + a3} end function function evaluate_bernstein_degree3(sequence bernstein_coefficients, atom t) atom {b0, b1, b2, b3} = bernstein_coefficients // de Casteljau’s algorithm. atom s = 1 - t, b01 = (s * b0) + (t * b1), b12 = (s * b1) + (t * b2), b23 = (s * b2) + (t * b3), b012 = (s * b01) + (t * b12), b123 = (s * b12) + (t * b23), b0123 = (s * b012) + (t * b123) return b0123 end function function bernstein_degree2_to_degree3(sequence bernstein_coefficients) atom {b0, b1, b2} = bernstein_coefficients return {b0, ((1/3) * b0) + ((2/3) * b1), ((2/3) * b1) + ((1/3) * b2), b2} end function function evaluate_monomial_degree2(sequence monomial_coefficients, atom t) atom {a0, a1, a2} = monomial_coefficients /* Horner’s rule. */ return a0 + (t * (a1 + (t * a2))) end function function evaluate_monomial_degree3(sequence monomial_coefficients, atom t) atom {a0, a1, a2, a3} = monomial_coefficients --/* Horner’s rule. --*/ return a0 + (t * (a1 + (t * (a2 + (t * a3))))) end function constant pmono2 = {1, 0, 0}, -- p(x) = 1 qmono2 = {1, 2, 3}, -- q(x) = 1 + 2x + 3x² pbern2 = monomial_to_bernstein_degree2(pmono2), qbern2 = monomial_to_bernstein_degree2(qmono2) printf(1,"Subprogram (1) examples:\n") printf(1," mono %v --> bern %v\n", {pmono2,pbern2}) printf(1," mono %v --> bern %v\n", {qmono2,qbern2}) printf(1,"Subprogram (2) examples:\n") for x in {0.25, 7.50} do printf(1," p(%g) = %g (mono %g)\n",{x,evaluate_bernstein_degree2(pbern2,x), evaluate_monomial_degree2(pmono2,x),}) printf(1," q(%g) = %g (mono %g)\n",{x,evaluate_bernstein_degree2(qbern2,x), evaluate_monomial_degree2(qmono2,x)}) end for constant pmono3 = {1, 0, 0, 0}, -- p(x) = 1 qmono3 = {1, 2, 3, 0}, -- q(x) = 1 + 2x + 3x² rmono3 = {1, 2, 3, 4}, -- r(x) = 1 + 2x + 3x² + 4x³ pbern3 = monomial_to_bernstein_degree3(pmono3), qbern3 = monomial_to_bernstein_degree3(qmono3), rbern3 = monomial_to_bernstein_degree3(rmono3) printf(1,"Subprogram (3) examples:\n") printf(1," mono %v --> bern %v\n", {pmono3,pbern3}) printf(1," mono %v --> bern %v\n", {qmono3,qbern3}) printf(1," mono %v --> bern %v\n", {rmono3,rbern3}) printf(1,"Subprogram (4) examples:\n") for x in {0.25, 7.50} do printf(1," p(%g) = %g (mono %g)\n",{x,evaluate_bernstein_degree3(pbern3,x), evaluate_monomial_degree3(pmono3,x)}) printf(1," q(%g) = %g (mono %g)\n",{x,evaluate_bernstein_degree3(qbern3,x), evaluate_monomial_degree3(qmono3,x)}) printf(1," r(%g) = %g (mono %g)\n",{x,evaluate_bernstein_degree3(rbern3,x), evaluate_monomial_degree3(rmono3,x)}) end for printf(1,"Subprogram (5) examples:\n") constant pbern3a = bernstein_degree2_to_degree3(pbern2), qbern3a = bernstein_degree2_to_degree3(qbern2) printf(1," bern %v --> bern %v\n", {pbern2,pbern3a}) printf(1," bern %v --> bern %v\n", {qbern2,qbern3a})
- 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) q(0.25) = 1.6875 (mono 1.6875) p(7.5) = 1 (mono 1) q(7.5) = 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.666666667,3.333333333,6} mono {1,2,3,4} --> bern {1,1.666666667,3.333333333,10} Subprogram (4) examples: p(0.25) = 1 (mono 1) q(0.25) = 1.6875 (mono 1.6875) r(0.25) = 1.75 (mono 1.75) p(7.5) = 1 (mono 1) q(7.5) = 184.75 (mono 184.75) r(7.5) = 1872.25 (mono 1872.25) Subprogram (5) examples: bern {1,1,1} --> bern {1,1,1,1} bern {1,2,6} --> bern {1,1.666666667,3.333333333,6}
Python
#!/bin/env python3
# Subprogram (1).
def monomial_to_bernstein_degree2(monomial_coefficients):
(a0, a1, a2) = monomial_coefficients
return (a0, a0 + ((1/2) * a1), a0 + a1 + a2)
# Subprogram (2).
def evaluate_bernstein_degree2(bernstein_coefficients, t):
(b0, b1, b2) = bernstein_coefficients
# de Casteljau’s algorithm.
s = 1 - t
b01 = (s * b0) + (t * b1)
b12 = (s * b1) + (t * b2)
b012 = (s * b01) + (t * b12)
return b012
# Subprogram (3).
def monomial_to_bernstein_degree3(monomial_coefficients):
(a0, a1, a2, a3) = monomial_coefficients
return (a0, a0 + ((1/3) * a1),
a0 + ((2/3) * a1) + ((1/3) * a2),
a0 + a1 + a2 + a3)
# Subprogram (4).
def evaluate_bernstein_degree3(bernstein_coefficients, t):
(b0, b1, b2, b3) = bernstein_coefficients
# de Casteljau’s algorithm.
s = 1 - t
b01 = (s * b0) + (t * b1)
b12 = (s * b1) + (t * b2)
b23 = (s * b2) + (t * b3)
b012 = (s * b01) + (t * b12)
b123 = (s * b12) + (t * b23)
b0123 = (s * b012) + (t * b123)
return b0123
# Subprogram (5).
def bernstein_degree2_to_degree3(bernstein_coefficients):
(b0, b1, b2) = bernstein_coefficients
return (b0, ((1/3) * b0) + ((2/3) * b1),
((2/3) * b1) + ((1/3) * b2), b2)
def evaluate_monomial_degree2(monomial_coefficients, t):
(a0, a1, a2) = monomial_coefficients
# Horner’s rule.
return a0 + (t * (a1 + (t * a2)))
def evaluate_monomial_degree3(monomial_coefficients, t):
(a0, a1, a2, a3) = monomial_coefficients
# Horner’s rule.
return 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.
#
pmono2 = (1, 0, 0)
qmono2 = (1, 2, 3)
pbern2 = monomial_to_bernstein_degree2(pmono2)
qbern2 = monomial_to_bernstein_degree2(qmono2)
print("Subprogram (1) examples:")
print(" mono", pmono2, " --> bern", pbern2)
print(" 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.
#
print("Subprogram (2) examples:")
for x in (0.25, 7.50):
print(" p(", x, ") = ", evaluate_bernstein_degree2(pbern2, x),
" ( mono: ", evaluate_monomial_degree2(pmono2, x), ")")
for x in (0.25, 7.50):
print(" q(", x, ") = ", evaluate_bernstein_degree2(qbern2, x),
" ( mono: ", evaluate_monomial_degree2(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.
#
pmono3 = (1, 0, 0, 0)
qmono3 = (1, 2, 3, 0)
rmono3 = (1, 2, 3, 4)
pbern3 = monomial_to_bernstein_degree3(pmono3)
qbern3 = monomial_to_bernstein_degree3(qmono3)
rbern3 = monomial_to_bernstein_degree3(rmono3)
print("Subprogram (3) examples:")
print(" mono", pmono3, " --> bern", pbern3)
print(" mono", qmono3, " --> bern", qbern3)
print(" 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.
#
print("Subprogram (4) examples:")
for x in (0.25, 7.50):
print(" p(", x, ") = ", evaluate_bernstein_degree3(pbern3, x),
" ( mono: ", evaluate_monomial_degree3(pmono3, x), ")")
for x in (0.25, 7.50):
print(" q(", x, ") = ", evaluate_bernstein_degree3(qbern3, x),
" ( mono: ", evaluate_monomial_degree3(qmono3, x), ")")
for x in (0.25, 7.50):
print(" r(", x, ") = ", evaluate_bernstein_degree3(rbern3, x),
" ( mono: ", evaluate_monomial_degree3(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.
#
print("Subprogram (5) examples:")
pbern3a = bernstein_degree2_to_degree3(pbern2)
qbern3a = bernstein_degree2_to_degree3(qbern2)
print(" bern", pbern2, " --> bern", pbern3a)
print(" bern", qbern2, " --> bern", qbern3a)
- Output:
Subprogram (1) examples: mono (1, 0, 0) --> bern (1, 1.0, 1) mono (1, 2, 3) --> bern (1, 2.0, 6) 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) --> bern (1, 1.0, 1.0, 1) mono (1, 2, 3, 0) --> bern (1, 1.6666666666666665, 3.333333333333333, 6) mono (1, 2, 3, 4) --> bern (1, 1.6666666666666665, 3.333333333333333, 10) 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: bern (1, 1.0, 1) --> bern (1, 1.0, 1.0, 1) bern (1, 2.0, 6) --> bern (1, 1.6666666666666665, 3.333333333333333, 6)
Wren
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.
import "./fmt" for Fmt
import "./math" for Math
var toBern2 = Fn.new { |a| [a[0], a[0] + a[1] / 2, a[0] + a[1] + a[2]] }
// uses de Casteljau's algorithm
var evalBern2 = Fn.new { |b, t|
var s = 1 - t
var b01 = s * b[0] + t * b[1]
var b12 = s * b[1] + t * b[2]
return s * b01 + t * b12
}
var toBern3 = Fn.new { |a|
var b = List.filled(4, 0)
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
var evalBern3 = Fn.new { |b, t|
var s = 1 - t
var b01 = s * b[0] + t * b[1]
var b12 = s * b[1] + t * b[2]
var b23 = s * b[2] + t * b[3]
var b012 = s * b01 + t * b12
var b123 = s * b12 + t * b23
return s * b012 + t * b123
}
var bern2to3 = Fn.new { |q|
var c = List.filled(4, 0)
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
}
var pm = [1, 0, 0]
var qm = [1, 2, 3]
var rm = [1, 2, 3, 4]
var x
var y
var m
System.print("Subprogram(1) examples:")
var pb2 = toBern2.call(pm)
var qb2 = toBern2.call(qm)
Fmt.print("mono $n --> bern $n", pm, pb2)
Fmt.print("mono $n --> bern $n", qm, qb2)
System.print("\nSubprogram(2) examples:")
x = 0.25
y = evalBern2.call(pb2, x)
m = Math.evalPoly(pm[-1..0], x)
Fmt.print("p($4.2f) = $j (mono $j)", x, y, m)
x = 7.5
y = evalBern2.call(pb2, x)
m = Math.evalPoly(pm[-1..0], x)
Fmt.print("p($4.2f) = $j (mono $j)", x, y, m)
x = 0.25
y = evalBern2.call(qb2, x)
m = Math.evalPoly(qm[-1..0], x)
Fmt.print("q($4.2f) = $j (mono $j)", x, y, m)
x = 7.5
y = evalBern2.call(qb2, x)
m = Math.evalPoly(qm[-1..0], x)
Fmt.print("q($4.2f) = $j (mono $j)", x, y, m)
System.print("\nSubprogram(3) examples:")
pm.add(0)
qm.add(0)
var pb3 = toBern3.call(pm)
var qb3 = toBern3.call(qm)
var rb3 = toBern3.call(rm)
Fmt.print("mono $n --> bern $n", pm, pb3)
Fmt.print("mono $n --> bern $n", qm, qb3)
Fmt.print("mono $n --> bern $n", rm, rb3)
System.print("\nSubprogram(4) examples:")
x = 0.25
y = evalBern3.call(pb3, x)
m = Math.evalPoly(pm[-1..0], x)
Fmt.print("p($4.2f) = $j (mono $j)", x, y, m)
x = 7.5
y = evalBern3.call(pb3, x)
m = Math.evalPoly(pm[-1..0], x)
Fmt.print("p($4.2f) = $j (mono $j)", x, y, m)
x = 0.25
y = evalBern3.call(qb3, x)
m = Math.evalPoly(qm[-1..0], x)
Fmt.print("q($4.2f) = $j (mono $j)", x, y, m)
x = 7.5
y = evalBern3.call(qb3, x)
m = Math.evalPoly(qm[-1..0], x)
Fmt.print("q($4.2f) = $j (mono $j)", x, y, m)
x = 0.25
y = evalBern3.call(rb3, x)
m = Math.evalPoly(rm[-1..0], x)
Fmt.print("r($4.2f) = $j (mono $j)", x, y, m)
x = 7.5
y = evalBern3.call(rb3, x)
m = Math.evalPoly(rm[-1..0], x)
Fmt.print("r($4.2f) = $j (mono $j)", x, y, m)
System.print("\nSubprogram(5) examples:")
var pc = bern2to3.call(pb2)
var qc = bern2to3.call(qb2)
Fmt.print("mono $n --> bern $n", pb2, pc)
Fmt.print("mono $n --> bern $n", qb2, qc)
- 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) 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]
XPL0
procedure ToBern2 (A0, A1, A2, B0, B1, B2);
real A0, A1, A2; \pass by value
real B0, B1, B2; \pass by reference
\Subprogram (1): transform monomial coefficients
\ A0, A1, A2 to Bernstein coefficients B0, B1, B2;
begin
B0(0) := A0;
B1(0) := A0 + ((1./2.) * A1);
B2(0) := A0 + A1 + A2
end \ToBern2\;
function real EvalBern2 (B0, B1, B2, T);
real B0, B1, B2, T;
\Subprogram (2): evaluate, at T, the polynomial with
\ Bernstein coefficients B0, B1, B2. Use de Casteljau's
\ algorithm;
real S, B01, B12, B012;
begin
S := 1. - T;
B01 := (S * B0) + (T * B1);
B12 := (S * B1) + (T * B2);
B012 := (S * B01) + (T * B12);
return B012
end \EvalBern2\;
procedure ToBern3 (A0, A1, A2, A3, B0, B1, B2, B3);
real A0, A1, A2, A3;
real B0, B1, B2, B3;
\Subprogram (3): transform monomial coefficients
\ A0, A1, A2, A3 to Bernstein coefficients B0, B1,
\ B2, B3;
begin
B0(0) := A0;
B1(0) := A0 + ((1./3.) * A1);
B2(0) := A0 + ((2./3.) * A1) + ((1./3.) * A2);
B3(0) := A0 + A1 + A2 + A3
end \ToBern3\;
function real EvalBern3 (B0, B1, B2, B3, T);
real B0, B1, B2, B3, T;
\Subprogram (4): evaluate, at t, the polynomial
\ with Bernstein coefficients b0, b1, b2, b3. Use
\ de Casteljau's algorithm;
real S, B01, B12, B23, B012, B123, B0123;
begin
S := 1. - T;
B01 := (S * B0) + (T * B1);
B12 := (S * B1) + (T * B2);
B23 := (S * B2) + (T * B3);
B012 := (S * B01) + (T * B12);
B123 := (S * B12) + (T * B23);
B0123 := (S * B012) + (T * B123);
return B0123
end \EvalBern3\;
procedure Bern2To3 (Q0, Q1, Q2, C0, C1, C2, C3);
real Q0, Q1, Q2;
real C0, C1, C2, C3;
\Subprogram (5): transform the quadratic Bernstein
\ coefficients q0, q1, q2 to the cubic Bernstein
\ coefficients c0, c1, c2, c3;
begin
C0(0) := Q0;
C1(0) := ((1./3.) * Q0) + ((2./3.) * Q1);
C2(0) := ((2./3.) * Q1) + ((1./3.) * Q2);
C3(0) := Q2
end \Bern2To3\;
real P0M, P1M, P2M;
real Q0M, Q1M, Q2M;
real R0M, R1M, R2M, R3M;
real P0B2, P1B2, P2B2;
real Q0B2, Q1B2, Q2B2;
real P0B3, P1B3, P2B3, P3B3;
real Q0B3, Q1B3, Q2B3, Q3B3;
real R0B3, R1B3, R2B3, R3B3;
real PC0, PC1, PC2, PC3;
real QC0, QC1, QC2, QC3;
real X, Y;
begin
P0M := 1.; P1M := 0.; P2M := 0.;
Q0M := 1.; Q1M := 2.; Q2M := 3.;
R0M := 1.; R1M := 2.; R2M := 3.; R3M := 4.;
Format(1, 2);
ToBern2 (P0M, P1M, P2M, @P0B2, @P1B2, @P2B2);
ToBern2 (Q0M, Q1M, Q2M, @Q0B2, @Q1B2, @Q2B2);
Text (0, "Subprogram (1) examples:^m^j");
Text (0, " mono ( ");
RlOut (0, P0M); Text (0, ", ");
RlOut (0, P1M); Text (0, ", ");
RlOut (0, P2M); Text (0, ") --> bern ( ");
RlOut (0, P0B2); Text (0, ", ");
RlOut (0, P1B2); Text (0, ", ");
RlOut (0, P2B2); Text (0, ")^m^j");
Text (0, " mono ( ");
RlOut (0, Q0M); Text (0, ", ");
RlOut (0, Q1M); Text (0, ", ");
RlOut (0, Q2M); Text (0, ") --> bern ( ");
RlOut (0, Q0B2); Text (0, ", ");
RlOut (0, Q1B2); Text (0, ", ");
RlOut (0, Q2B2); Text (0, ")^m^j");
Text (0, "Subprogram (2) examples:^m^j");
X := 0.25;
Y := EvalBern2 (P0B2, P1B2, P2B2, X);
Text (0, " p ( "); RlOut (0, X);
Text (0, ") = "); RlOut (0, Y);
Text (0, "^m^j");
X := 7.50;
Y := EvalBern2 (P0B2, P1B2, P2B2, X);
Text (0, " p ( "); RlOut (0, X);
Text (0, ") = "); RlOut (0, Y);
Text (0, "^m^j");
X := 0.25;
Y := EvalBern2 (Q0B2, Q1B2, Q2B2, X);
Text (0, " q ( "); RlOut (0, X);
Text (0, ") = "); RlOut (0, Y);
Text (0, "^m^j");
X := 7.50;
Y := EvalBern2 (Q0B2, Q1B2, Q2B2, X);
Text (0, " q ( "); RlOut (0, X);
Text (0, ") = "); RlOut (0, Y);
Text (0, "^m^j");
ToBern3 (P0M, P1M, P2M, 0., @P0B3, @P1B3, @P2B3, @P3B3);
ToBern3 (Q0M, Q1M, Q2M, 0., @Q0B3, @Q1B3, @Q2B3, @Q3B3);
ToBern3 (R0M, R1M, R2M, R3M, @R0B3, @R1B3, @R2B3, @R3B3);
Text (0, "Subprogram (3) examples:^m^j");
Text (0, " mono ( ");
RlOut (0, P0M); Text (0, ", ");
RlOut (0, P1M); Text (0, ", ");
RlOut (0, P2M); Text (0, ", ");
RlOut (0, 0.); Text (0, ") --> bern ( ");
RlOut (0, P0B3); Text (0, ", ");
RlOut (0, P1B3); Text (0, ", ");
RlOut (0, P2B3); Text (0, ", ");
RlOut (0, P3B3); Text (0, ")^m^j");
Text (0, " mono ( ");
RlOut (0, Q0M); Text (0, ", ");
RlOut (0, Q1M); Text (0, ", ");
RlOut (0, Q2M); Text (0, ", ");
RlOut (0, 0.); Text (0, ") --> bern ( ");
RlOut (0, Q0B3); Text (0, ", ");
RlOut (0, Q1B3); Text (0, ", ");
RlOut (0, Q2B3); Text (0, ", ");
RlOut (0, Q3B3); Text (0, ")^m^j");
Text (0, " mono ( ");
RlOut (0, R0M); Text (0, ", ");
RlOut (0, R1M); Text (0, ", ");
RlOut (0, R2M); Text (0, ", ");
RlOut (0, R3M); Text (0, ") --> bern ( ");
RlOut (0, R0B3); Text (0, ", ");
RlOut (0, R1B3); Text (0, ", ");
RlOut (0, R2B3); Text (0, ", ");
RlOut (0, R3B3); Text (0, ")^m^j");
Text (0, "Subprogram (4) examples:^m^j");
X := 0.25;
Y := EvalBern3 (P0B3, P1B3, P2B3, P3B3, X);
Text (0, " p ( "); RlOut (0, X);
Text (0, ") = "); RlOut (0, Y);
Text (0, "^m^j");
X := 7.50;
Y := EvalBern3 (P0B3, P1B3, P2B3, P3B3, X);
Text (0, " p ( "); RlOut (0, X);
Text (0, ") = "); RlOut (0, Y);
Text (0, "^m^j");
X := 0.25;
Y := EvalBern3 (Q0B3, Q1B3, Q2B3, Q3B3, X);
Text (0, " q ( "); RlOut (0, X);
Text (0, ") = "); RlOut (0, Y);
Text (0, "^m^j");
X := 7.50;
Y := EvalBern3 (Q0B3, Q1B3, Q2B3, Q3B3, X);
Text (0, " q ( "); RlOut (0, X);
Text (0, ") = "); RlOut (0, Y);
Text (0, "^m^j");
X := 0.25;
Y := EvalBern3 (R0B3, R1B3, R2B3, R3B3, X);
Text (0, " r ( "); RlOut (0, X);
Text (0, ") = "); RlOut (0, Y);
Text (0, "^m^j");
X := 7.50;
Y := EvalBern3 (R0B3, R1B3, R2B3, R3B3, X);
Text (0, " r ( "); RlOut (0, X);
Text (0, ") = "); RlOut (0, Y);
Text (0, "^m^j");
Bern2To3 (P0B2, P1B2, P2B2, @PC0, @PC1, @PC2, @PC3);
Bern2To3 (Q0B2, Q1B2, Q2B2, @QC0, @QC1, @QC2, @QC3);
Text (0, "Subprogram (5) examples:^m^j");
Text (0, " bern ( ");
RlOut (0, P0B2); Text (0, ", ");
RlOut (0, P1B2); Text (0, ", ");
RlOut (0, P2B2); Text (0, ") --> bern ( ");
RlOut (0, PC0); Text (0, ", ");
RlOut (0, PC1); Text (0, ", ");
RlOut (0, PC2); Text (0, ", ");
RlOut (0, PC3); Text (0, ")^m^j");
Text (0, " bern ( ");
RlOut (0, Q0B2); Text (0, ", ");
RlOut (0, Q1B2); Text (0, ", ");
RlOut (0, Q2B2); Text (0, ") --> bern ( ");
RlOut (0, QC0); Text (0, ", ");
RlOut (0, QC1); Text (0, ", ");
RlOut (0, QC2); Text (0, ", ");
RlOut (0, QC3); Text (0, ")^m^j")
end
- Output:
Subprogram (1) examples: mono ( 1.00, 0.00, 0.00) --> bern ( 1.00, 1.00, 1.00) mono ( 1.00, 2.00, 3.00) --> bern ( 1.00, 2.00, 6.00) Subprogram (2) examples: p ( 0.25) = 1.00 p ( 7.50) = 1.00 q ( 0.25) = 1.69 q ( 7.50) = 184.75 Subprogram (3) examples: mono ( 1.00, 0.00, 0.00, 0.00) --> bern ( 1.00, 1.00, 1.00, 1.00) mono ( 1.00, 2.00, 3.00, 0.00) --> bern ( 1.00, 1.67, 3.33, 6.00) mono ( 1.00, 2.00, 3.00, 4.00) --> bern ( 1.00, 1.67, 3.33, 10.00) Subprogram (4) examples: p ( 0.25) = 1.00 p ( 7.50) = 1.00 q ( 0.25) = 1.69 q ( 7.50) = 184.75 r ( 0.25) = 1.75 r ( 7.50) = 1872.25 Subprogram (5) examples: 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)