Bernstein basis polynomials
You are encouraged to solve this task according to the task description, using any language you may know.
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.
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 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; 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 )
ALGOL 68
Algol 68 doesn't have procedure overloading but does have monadic and dyadic operator overloading and user defined operator symbols.
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 BERNTWO b )VOID:
BEGIN
b0 OF b := a0;
b1 OF b := a0 + ((1/2) * a1);
b2 OF b := a0 + a1 + a2
END # tobern2 # ;
# Subprogram (2): evaluate, at t, the polynomial with #
# Bernstein coefficients b0, b1, b2. Use de Casteljau's #
# algorithm #
OP EVAL = ( BERNTWO b, REAL t )REAL:
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 # EVAL # ;
# Subprogram (3): transform monomial coefficients #
# a0, a1, a2, a3 to Bernstein coefficients b0, b1, #
# b2, b3 #
PROC tobern3 = ( REAL a0, a1, a2, a3, REF BERNTHREE b )VOID:
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 # ;
# Subprogram (4): evaluate, at t, the polynomial #
# with Bernstein coefficients b0, b1, b2, b3. Use #
# de Casteljau's algorithm #
OP EVAL = ( BERNTHREE b, REAL t )REAL:
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 # EVAL #;
# 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
HEAP BERNTHREE c;
b0 OF c := b0 OF q;
b1 OF c := ((1/3) * b0 OF q) + ((2/3) * b1 OF q);
b2 OF c := ((2/3) * b1 OF q) + ((1/3) * b2 OF q);
b3 OF c := b2 OF q;
c
END # TOBERNTHREE # ;
BEGIN
# returns x as a string without trailing 0 decoimals #
PROC f = ( REAL x )STRING:
BEGIN
STRING v := fixed( x, -14, 11 );
INT end pos := UPB v;
WHILE IF end pos < LWB v THEN FALSE ELSE v[ end pos ] = "0" FI DO
end pos -:= 1
OD;
IF end pos >= LWB v THEN
IF v[ end pos ] = "." THEN end pos -:= 1 FI
FI;
INT start pos := LWB v;
WHILE IF start pos > end pos THEN FALSE ELSE v[ start pos ] = " " FI DO
start pos +:= 1
OD;
IF end pos < start pos THEN "0" ELSE v[ start pos : end pos ] FI
END # f # ;
PRIO SHOWEVAL = 5;
# prints the result of evaluating p with x #
OP SHOWEVAL = ( BERNTWO p, REAL x )VOID:
print( ( " ", label OF p, " ( ", f( x ), " ) = ", f( p EVAL x ), newline ) );
# prints the result of evaluating p with x #
OP SHOWEVAL = ( BERNTHREE p, REAL x )VOID:
print( ( " ", label OF p, " ( ", f( x ), " ) = ", f( p EVAL x ), newline ) );
# returns a string representation of the values of p #
OP TOSTRING = ( BERNTWO p )STRING:
"( " + 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 p2 := BERNTWO ( "p", 0, 0, 0 );
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;
REAL q0m = 1, q1m = 2, q2m = 3;
REAL r0m = 1, r1m = 2, r2m = 3, r3m = 4;
tobern2( p0m, p1m, p2m, p2 );
tobern2( q0m, q1m, q2m, q2 );
print( ( "Subprogram (1) examples:", newline ) );
print( ( " mono ", TOSTRING []REAL( p0m, p1m, p2m ), " --> bern ", TOSTRING p2, newline ) );
print( ( " mono ", TOSTRING []REAL( q0m, q1m, q2m ), " --> bern ", TOSTRING q2, newline ) );
print( ( "Subprogram (2) examples:", newline ) );
p2 SHOWEVAL 0.25;
p2 SHOWEVAL 7.50;
q2 SHOWEVAL 0.25;
q2 SHOWEVAL 7.50;
tobern3( p0m, p1m, p2m, 0, p3 );
tobern3( q0m, q1m, q2m, 0, q3 );
tobern3( r0m, r1m, r2m, r3m, r3 );
print( ( "Subprogram (3) examples:", newline ) );
print( ( " mono ", TOSTRING []REAL( p0m, p1m, p2m, 0 ), " --> bern ", TOSTRING p3, newline ) );
print( ( " mono ", TOSTRING []REAL( q0m, q1m, q2m, 0 ), " --> bern ", TOSTRING q3, newline ) );
print( ( " mono ", TOSTRING []REAL( r0m, r1m, r2m, r3m ), " --> bern ", TOSTRING r3, newline ) );
print( ( "Subprogram (4) examples:", newline ) );
p3 SHOWEVAL 0.25;
p3 SHOWEVAL 7.50;
q3 SHOWEVAL 0.25;
q3 SHOWEVAL 7.50;
r3 SHOWEVAL 0.25;
r3 SHOWEVAL 7.50;
print( ( "Subprogram (5) examples:", newline ) );
print( ( " bern ", TOSTRING p2, " --> bern ", TOSTRING TOBERNTHREE p2, newline ) );
print( ( " bern ", TOSTRING q2, " --> bern ", TOSTRING TOBERNTHREE q2, newline ) )
END
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 )
ATS
The following is full of pedagogical comments, which I hope make this mysterious programming language a bit less mysterious. ATS is a very complicated language, whose typechecker does things most language's typecheckers simply do not do.
(The template system also is complicated, and in fact can cause bad C code to be generated. I view this issue as one of bugginess that a determined programmer works around. It is like gfortran not warning you that your misuse of "ALLOCATABLE" is not going to succeed. The compiler writers have concentrated their efforts on other matters.)
(* The following can be compiled with "patscc -O3
source_file_name.dats", or with your preferred optimizer
setting. The ATS compiler generates C that assumes you will do
aggressive optimization, say "-O2" or better.
I like to add GCC's "-std=gnu2x" flag, as well, though I expect
this flag will at some point be changed to "-std=gnu23" or
similar. Otherwise (with the default patscc settings) you get
"-std=c99".
For adventure, try "patscc -std=gnu2x -Ofast -march=native
source_file_name.dats" *)
#include "share/atspre_staload.hats"
(* 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
macdef two_thirds = g0i2f 2 / g0i2f 3
(* We will try to achieve the type-safety of distinguishing between
monomial and Bernstein bases at compile time, but without overhead
at runtime. (Such type-safety can be partly achieved in many other
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"
dataprop BASIS (polyknd : tkind) =
| MONOMIAL_BASIS (monoknd) of ()
| BERNSTEIN_BASIS (bernknd) of ()
typedef poly_degree2 (polyknd : tkind, coefknd : tkind) =
@(BASIS polyknd | g0float coefknd, g0float coefknd,
g0float coefknd)
typedef poly_degree3 (polyknd : tkind, coefknd : tkind) =
@(BASIS polyknd | g0float coefknd, g0float coefknd,
g0float coefknd, g0float coefknd)
(* The :<> below means that the function is proven to terminate, does
not raise exceptions, does not access or write to memory it does
not "own", etc. In practice, however, such requirements are often
masked. The goal of ATS is not formal proof, but bug reduction,
easier debugging, avoidance of the need for runtime checks, etc. *)
fn {coefknd : tkind} (* Subprogram (1) *)
mono2bern_degree2 (coefs : poly_degree2 (monoknd, coefknd))
:<> poly_degree2 (bernknd, coefknd) =
(* AN EXAMPLE OF ATS2 TYPE-SAFETY: If you try to call
mono2bern_degree2 with Bernstein coefficients (instead of monomial
coefficients) as the argument, then, AT COMPILE TIME, you get
messages similar to these:
/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.
/some_path/bernstein_poly_task.dats: 1051(line=32, offs=31) -- 1052(line=32, offs=32): error(3): unsolved constraint: C3NSTRprop(C3TKmain(); S2Eeqeq(S2Eextkind(rosettacode_bernstein_poly); S2Eextkind(rosettacode_monomial_poly)))
typechecking has failed: there are some unsolved constraints: please inspect the above reported error message(s) for information.
exit(ATS): uncaught exception: _2tmp_2ATS_2dPostiats_2src_2pats_error_2esats__FatalErrorExn(1025)
AT RUNTIME, however, both kinds of polynomial coefficients are
represented by a plain unboxed tuple (that is, by a C struct) of
three numbers, each being of whichever type "g0float coefknd"
happens to be ("float", "double", "ldouble", etc.). *)
let
val+ @(_ | a0, a1, a2) = coefs
in
@(BERNSTEIN_BASIS () | a0, a0 + (one_half * a1), a0 + a1 + a2)
end
fn {coefknd : tkind} (* Subprogram (2) *)
evalbern_degree2 (coefs : poly_degree2 (bernknd, coefknd),
t : g0float coefknd)
:<> g0float coefknd =
let
val+ @(_ | b0, b1, b2) = coefs
and s = g0i2f 1 - t
val b01 = (s * b0) + (t * b1)
and b12 = (s * b1) + (t * b2)
val b012 = (s * b01) + (t * b12)
in
b012
end
fn {coefknd : tkind} (* Subprogram (3) *)
mono2bern_degree3 (coefs : poly_degree3 (monoknd, coefknd))
:<> poly_degree3 (bernknd, coefknd) =
let
val+ @(_ | a0, a1, a2, a3) = coefs
in
@(BERNSTEIN_BASIS () | a0, a0 + (one_third * a1),
a0 + (two_thirds * a1) + (one_third * a2),
a0 + a1 + a2 + a3)
end
fn {coefknd : tkind} (* Subprogram (4) *)
evalbern_degree3 (coefs : poly_degree3 (bernknd, coefknd),
t : g0float coefknd)
:<> g0float coefknd =
let
val+ @(_ | b0, b1, b2, b3) = coefs
and s = g0i2f 1 - t
val b01 = (s * b0) + (t * b1)
and b12 = (s * b1) + (t * b2)
and b23 = (s * b2) + (t * b3)
val b012 = (s * b01) + (t * b12)
and b123 = (s * b12) + (t * b23)
val b0123 = (s * b012) + (t * b123)
in
b0123
end
fn {coefknd : tkind} (* Subprogram (5) *)
bern_degree2_to_3 (coefs : poly_degree2 (bernknd, coefknd))
:<> poly_degree3 (bernknd, coefknd) =
let
val+ @(_ | q0, q1, q2) = coefs
in
@(BERNSTEIN_BASIS () | q0, (one_third * q0) + (two_thirds * q1),
(two_thirds * q1) + (one_third * q2), q2)
end
(* Let us make it easy to print out coefficients. Note that
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
val+ @(_ | a0, a1, a2) = coefs
macdef outf = stdout_ref
in
fprint_val<string> (outf, "(");
fprint_val<g0float coefknd> (outf, a0);
fprint_val<string> (outf, ", ");
fprint_val<g0float coefknd> (outf, a1);
fprint_val<string> (outf, ", ");
fprint_val<g0float coefknd> (outf, a2);
fprint_val<string> (outf, ")")
end
fn {coefknd : tkind}
print_poly_degree3 {polyknd : tkind}
(coefs : poly_degree3 (polyknd, coefknd))
: void =
let
val+ @(_ | a0, a1, a2, a3) = coefs
macdef outf = stdout_ref
in
fprint_val<string> (outf, "(");
fprint_val<g0float coefknd> (outf, a0);
fprint_val<string> (outf, ", ");
fprint_val<g0float coefknd> (outf, a1);
fprint_val<string> (outf, ", ");
fprint_val<g0float coefknd> (outf, a2);
fprint_val<string> (outf, ", ");
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
fn {coefknd : tkind}
print_bern_degree3 (coefs : poly_degree3 (bernknd, coefknd))
: void =
begin
fprint_val<string> (stdout_ref, "bern ");
print_poly_degree3 (coefs)
end
overload print with print_mono_degree2
overload print with print_bern_degree2
overload print with print_mono_degree3
overload print with print_bern_degree3
fn {coefknd : tkind}
evalmono_degree2 (coefs : poly_degree2 (monoknd, coefknd),
t : g0float coefknd)
:<> g0float coefknd =
let
val+ @(_ | a0, a1, a2) = coefs
in
a0 + (t * (a1 + (t * a2))) (* Horner form. *)
end
fn {coefknd : tkind}
evalmono_degree3 (coefs : poly_degree3 (monoknd, coefknd),
t : g0float coefknd)
:<> g0float coefknd =
let
val+ @(_ | a0, a1, a2, a3) = coefs
in
a0 + (t * (a1 + (t * (a2 + (t * a3))))) (* Horner form. *)
end
fn
do_examples () : void =
let
(* I am going to use type "float = g0float fltknd". That is, C
single precision. The notation for such a number ends in an
"F" or "f", just as in C. *)
val pmono2 = @(MONOMIAL_BASIS () | 1.0F, 0.0F, 0.0F)
val qmono2 = @(MONOMIAL_BASIS () | 1.0F, 2.0F, 3.0F)
val pbern2 = mono2bern_degree2 (pmono2)
(* ATS often lets you leave out the argument parentheses, as in
the following line. However, unlike in OCaml or Standard ML,
(where parentheses also are often left out) argument
parentheses do NOT indicate that arguments are being passed as
a tuple. *)
val qbern2 = mono2bern_degree2 qmono2
(* ATS will often (if it might not be confused with argument
parentheses, etc.) let you leave out the "@" in an unboxed
tuple. *)
val pmono3 = (MONOMIAL_BASIS () | 1.0F, 0.0F, 0.0F, 0.0F)
val qmono3 = (MONOMIAL_BASIS () | 1.0F, 2.0F, 3.0F, 0.0F)
val rmono3 = (MONOMIAL_BASIS () | 1.0F, 2.0F, 3.0F, 4.0F)
val pbern3 = mono2bern_degree3 (pmono3)
val qbern3 = mono2bern_degree3 (qmono3)
val rbern3 = mono2bern_degree3 (rmono3)
val pbern3a = bern_degree2_to_3 (pbern2)
val qbern3a = bern_degree2_to_3 (qbern2)
in
println! ("Subprogram (1) examples:");
println! (" ", pmono2, " --> ", pbern2);
println! (" ", qmono2, " --> ", qbern2);
println! ("Subprogram (2) examples:");
let
(* The following is written in ATS-isms, showing how to do
things (a) in procedural fashion, (b) without an allocator,
(c) without runtime checks, and yet (d) without risk of going
outside the bounds of an array.
Figuring out how I am doing all that is left partly as a
boring exercise for the reader. (Be warned that the use of
for* is not, as far as I know, documented anywhere but in
mailing list archives.)
xs is an array on the stack, and i is a loop variable on the
stack or in a register. The for* loop works like a C
for-loop, except it is proven to terminate, and it guarantees
that xs is never indexed with any value except 0 or 1.
An alternative would have been to use a regular "for" loop
(without the asterisk) and do runtime bounds checks (raising
an exception on check failure). The compiler will NOT let you
index xs with an out-of-bounds index. *)
var xs : array (float, 2) = @[float][2] (0.25F, 7.50F)
var i : Int
in
for* {i : nat | i <= 2} .<2 - i>. (i : int i) =>
(i := 0; i <> 2; i := i + 1)
let
val x = xs[i]
val ypm = evalmono_degree2 (pmono2, x)
and ypb = evalbern_degree2 (pbern2, x)
in
println! (" p(", x, ") = ", ypb, " (mono: ", ypm, ")")
end;
for* {i : nat | i <= 2} .<2 - i>. (i : int i) =>
(i := 0; i <> 2; i := i + 1)
let
val x = xs[i]
val yqm = evalmono_degree2 (qmono2, x)
and yqb = evalbern_degree2 (qbern2, x)
in
println! (" q(", x, ") = ", yqb, " (mono: ", yqm, ")")
end
end;
println! ("Subprogram (3) examples:");
println! (" ", pmono3, " --> ", pbern3);
println! (" ", qmono3, " --> ", qbern3);
println! (" ", rmono3, " --> ", rbern3);
println! ("Subprogram (4) examples:");
let
var xs : array (float, 2) = @[float][2] (0.25F, 7.50F)
var i : Int
in
for* {i : nat | i <= 2} .<2 - i>. (i : int i) =>
(i := 0; i <> 2; i := i + 1)
let
val x = xs[i]
val ypm = evalmono_degree3 (pmono3, x)
and ypb = evalbern_degree3 (pbern3, x)
in
println! (" p(", x, ") = ", ypb, " (mono: ", ypm, ")")
end;
for* {i : nat | i <= 2} .<2 - i>. (i : int i) =>
(i := 0; i <> 2; i := i + 1)
let
val x = xs[i]
val yqm = evalmono_degree3 (qmono3, x)
and yqb = evalbern_degree3 (qbern3, x)
in
println! (" q(", x, ") = ", yqb, " (mono: ", yqm, ")")
end;
for* {i : nat | i <= 2} .<2 - i>. (i : int i) =>
(i := 0; i <> 2; i := i + 1)
let
val x = xs[i]
val yrm = evalmono_degree3 (rmono3, x)
and yrb = evalbern_degree3 (rbern3, x)
in
println! (" r(", x, ") = ", yrb, " (mono: ", yrm, ")")
end
end;
println! ("Subprogram (5) examples:");
println! (" ", pbern2, " --> ", pbern3a);
(* Although ATS is a "semicolon-as-separator" language, the
compiler will not complain about the null statement follow this
last semicolon. The ATS language's designer usually has put in
the semicolon, whereas I, being a bit of a purist, usually
leave it out. (And, yes, I do often get errors due to missing
semicolons; in ATS one catches that at compile time. And, at my
work 30 years ago, we programmed in Turbo/Borland Pascal and we
usually put the semicolon in.) *)
println! (" ", qbern2, " --> ", qbern3a);
end
implement main0 () = do_examples ()- Output:
Subprogram (1) examples: mono (1.000000, 0.000000, 0.000000) --> bern (1.000000, 1.000000, 1.000000) mono (1.000000, 2.000000, 3.000000) --> bern (1.000000, 2.000000, 6.000000) Subprogram (2) examples: p(0.250000) = 1.000000 (mono: 1.000000) p(7.500000) = 1.000000 (mono: 1.000000) q(0.250000) = 1.687500 (mono: 1.687500) q(7.500000) = 184.750000 (mono: 184.750000) Subprogram (3) examples: mono (1.000000, 0.000000, 0.000000, 0.000000) --> bern (1.000000, 1.000000, 1.000000, 1.000000) mono (1.000000, 2.000000, 3.000000, 0.000000) --> bern (1.000000, 1.666667, 3.333333, 6.000000) mono (1.000000, 2.000000, 3.000000, 4.000000) --> bern (1.000000, 1.666667, 3.333333, 10.000000) Subprogram (4) examples: p(0.250000) = 1.000000 (mono: 1.000000) p(7.500000) = 1.000000 (mono: 1.000000) q(0.250000) = 1.687500 (mono: 1.687500) q(7.500000) = 184.749771 (mono: 184.750000) r(0.250000) = 1.750000 (mono: 1.750000) r(7.500000) = 1872.250000 (mono: 1872.250000) Subprogram (5) examples: 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)
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;
}
- 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}
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)}]");
}
}
- 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] --> 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]
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;
}
- Output:
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]
FreeBASIC
Dim Shared b0 As Double, b1 As Double, b2 As Double
Sub tobern2 (Byval a0 As Double, Byval a1 As Double, Byval a2 As Double, _
Byref b0 As Double, Byref b1 As Double, Byref b2 As Double)
' Subprogram (1): transform monomial coefficients
' a0, a1, a2 to Bernstein coefficients b0, b1, b2
b0 = a0
b1 = a0 + ((1/2) * a1)
b2 = a0 + a1 + a2
End Sub
Function evalbern2 (b0 As Double, b1 As Double, b2 As Double, t As Double) As Double
' Subprogram (2): evaluate, at t, the polynomial with
' Bernstein coefficients b0, b1, b2. Use de Casteljau's algorithm
Dim As Double s, b01, b12, b012
s = 1 - t
b01 = (s * b0) + (t * b1)
b12 = (s * b1) + (t * b2)
b012 = (s * b01) + (t * b12)
Return b012
End Function
Sub tobern3 (Byval a0 As Double, Byval a1 As Double, Byval a2 As Double, Byval a3 As Double, _
Byref b0 As Double, Byref b1 As Double, Byref b2 As Double, Byref b3 As Double)
' Subprogram (3): transform monomial coefficients
' a0, a1, a2, a3 to Bernstein coefficients b0, b1, b2, b3
b0 = a0
b1 = a0 + ((1/3) * a1)
b2 = a0 + ((2/3) * a1) + ((1/3) * a2)
b3 = a0 + a1 + a2 + a3
End Sub
Function evalbern3 (b0 As Double, b1 As Double, b2 As Double, _
b3 As Double, t As Double) As Double
' Subprogram (4): evaluate, at t, the polynomial
' with Bernstein coefficients b0, b1, b2, b3.
' Use de Casteljau's algorithm
Dim As Double 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)
Return b0123
End Function
Sub bern2to3 (Byval q0 As Double, Byval q1 As Double, Byval q2 As Double, _
Byref c0 As Double, Byref c1 As Double, Byref c2 As Double, Byref c3 As Double)
' Subprogram (5): transform the quadratic Bernstein
' coefficients q0, q1, q2 to the cubic Bernstein
' coefficients c0, c1, c2, c3
c0 = q0
c1 = ((1/3) * q0) + ((2/3) * q1)
c2 = ((2/3) * q1) + ((1/3) * q2)
c3 = q2
End Sub
Dim As Double p0b2, p1b2, p2b2
Dim As Double q0b2, q1b2, q2b2
Dim As Double p0b3, p1b3, p2b3, p3b3
Dim As Double q0b3, q1b3, q2b3, q3b3
Dim As Double r0b3, r1b3, r2b3, r3b3
Dim As Double pc0, pc1, pc2, pc3
Dim As Double qc0, qc1, qc2, qc3
Dim As Double x, y
Dim As Double p0m = 1, p1m = 0, p2m = 0
Dim As Double q0m = 1, q1m = 2, q2m = 3
Dim As Double r0m = 1, r1m = 2, r2m = 3, r3m = 4
tobern2 (p0m, p1m, p2m, p0b2, p1b2, p2b2)
tobern2 (q0m, q1m, q2m, q0b2, q1b2, q2b2)
Print "Subprogram (1) examples:"
Print Using " mono (#.##, #.##, #.##) --> bern (#.##, #.##, #.##)"; p0m; p1m; p2m; p0b2; p1b2; p2b2
Print Using " mono (#.##, #.##, #.##) --> bern (#.##, #.##, #.##)"; q0m; q1m; q2m; q0b2; q1b2; q2b2
Print "Subprogram (2) examples:"
x = 0.25
y = evalbern2 (p0b2, p1b2, p2b2, x)
Print Using " p (#.##) = ###.##"; x; y
x = 7.50
y = evalbern2 (p0b2, p1b2, p2b2, x)
Print Using " p (#.##) = ###.##"; x; y
x = 0.25
y = evalbern2 (q0b2, q1b2, q2b2, x)
Print Using " q (#.##) = ###.##"; x; y
x = 7.50
y = evalbern2 (q0b2, q1b2, q2b2, x)
Print Using " q (#.##) = ###.##"; x; y
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)
Print "Subprogram (3) examples:"
Print Using " mono (#.##, #.##, #.##, 0.00) --> bern (#.##, #.##, #.##, ##.##)"; p0m; p1m; p2m; p0b3; p1b3; p2b3; p3b3
Print Using " mono (#.##, #.##, #.##, 0.00) --> bern (#.##, #.##, #.##, ##.##)"; q0m; q1m; q2m; q0b3; q1b3; q2b3; q3b3
Print Using " mono (#.##, #.##, #.##, #.##) --> bern (#.##, #.##, #.##, ##.##)"; r0m; r1m; r2m; r3m; r0b3; r1b3; r2b3; r3b3
Print "Subprogram (4) examples:"
x = 0.25
y = evalbern3 (p0b3, p1b3, p2b3, p3b3, x)
Print Using " p (#.##) = ####.##"; x; y
x = 7.50
y = evalbern3 (p0b3, p1b3, p2b3, p3b3, x)
Print Using " p (#.##) = ####.##"; x; y
x = 0.25
y = evalbern3 (q0b3, q1b3, q2b3, q3b3, x)
Print Using " q (#.##) = ####.##"; x; y
x = 7.50
y = evalbern3 (q0b3, q1b3, q2b3, q3b3, x)
Print Using " q (#.##) = ####.##"; x; y
x = 0.25
y = evalbern3 (r0b3, r1b3, r2b3, r3b3, x)
Print Using " r (#.##) = ####.##"; x; y
x = 7.50
y = evalbern3 (r0b3, r1b3, r2b3, r3b3, x)
Print Using " r (#.##) = ####.##"; x; y
bern2to3 (p0b2, p1b2, p2b2, pc0, pc1, pc2, pc3)
bern2to3 (q0b2, q1b2, q2b2, qc0, qc1, qc2, qc3)
Print "Subprogram (5) examples:"
Print Using " bern (#.##, #.##, #.##) --> bern (#.##, #.##, #.##, #.##)"; p0b2; p1b2; p2b2; pc0; pc1; pc2; pc3
Print Using " bern (#.##, #.##, #.##) --> bern (#.##, #.##, #.##, #.##)"; q0b2; q1b2; q2b2; qc0; qc1; qc2; qc3
Sleep- 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)
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)
}
- 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]
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 ) ) ) ) );
}
}
- Output:
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]
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.
M4
This works (insofar as it does work) with every m4 I have, including the "legacy" Heirloom Development Tools m4.
divert(-1)
#
# POSIX m4 has only "signed integer arithmetic with at least 32-bit
# precision", so I use integers and scaling.
#
# Be aware that overflows might be silently ignored by m4, resulting
# in nonsense being printed. I demonstrate this below (where the value
# of the polynomial exceeds four times the cube of 7.50).
#
# (Were I trying to implement this task more reliably, it would be by
# first implementing decimal arithmetic in m4, with a large or
# arbitrary number of digits. Which is something I have not, to date,
# ever done.)
#
define(`printnum',`eval(($1)/10000).eval(($1)%10000)')
define(`printpoly2',`_$0($1)')
define(`_printpoly2',
`(printnum($1), printnum($2), printnum($3))')
define(`printpoly3',`_$0($1)')
define(`_printpoly3',
`(printnum($1), printnum($2), printnum($3), printnum($4))')
# Subprogram (1).
define(`tobern2',`_$0($1)')
define(`_tobern2',
`$1,eval(($1) + (($2) * 5)/10),eval(($1) + ($2) + ($3))')
# Subprogram (2).
define(`evalbern2',`_$0($1,$2)')
define(`_evalbern2',
`pushdef(`t',eval(($4)))`'pushdef(`s',eval(100 - t))dnl
pushdef(`b01',eval((s * ($1))/100 + (t * ($2))/100))dnl
pushdef(`b12',eval((s * ($2))/100 + (t * ($3))/100))dnl
eval((s * b01)/100 + (t * b12)/100)dnl
popdef(`s',`t',`b01',`b12')')
# Subprogram (3).
define(`tobern3',`_$0($1)')
define(`_tobern3',
`$1,eval(($1) + (($2) * 3333)/10000),dnl
eval(($1) + (($2) * 6667)/10000 + (($3) * 3333)/10000),dnl
eval(($1) + ($2) + ($3) + ($4))')
# Subprogram (4).
define(`evalbern3',`_$0($1,$2)')
define(`_evalbern3',
`pushdef(`t',eval(($5)))`'pushdef(`s',eval(100 - t))dnl
pushdef(`b01',eval((s * ($1))/100 + (t * ($2))/100))dnl
pushdef(`b12',eval((s * ($2))/100 + (t * ($3))/100))dnl
pushdef(`b23',eval((s * ($3))/100 + (t * ($4))/100))dnl
pushdef(`b012',eval((s * b01)/100 + (t * b12)/100))dnl
popdef(`b01')dnl
pushdef(`b123',eval((s * b12)/100 + (t * b23)/100))dnl
popdef(`b12',`b23')dnl
eval((s * b012)/100 + (t * b123)/100)dnl
popdef(`s',`t',`b012',`b123')')
# Subprogram (5).
define(`bern2to3',`_$0($1)')
define(`_bern2to3',
`$1,eval((($1) * 3333)/10000 + (($2) * 6667)/10000),dnl
eval((($2) * 6667)/10000 + (($3) * 3333)/10000),$3')
define(`pmono2',``10000,00000,00000'')
define(`qmono2',``10000,20000,30000'')
define(`pbern2',`tobern2(pmono2)')
define(`qbern2',`tobern2(qmono2)')
define(`pmono3',``10000,00000,00000,00000'')
define(`qmono3',``10000,20000,30000,00000'')
define(`rmono3',``10000,20000,30000,40000'')
define(`pbern3',`tobern3(pmono3)')
define(`qbern3',`tobern3(qmono3)')
define(`rbern3',`tobern3(rmono3)')
define(`pbern3a',`bern2to3(`pbern2')')
define(`qbern3a',`bern2to3(`qbern2')')
divert`'dnl
Subprogram (1) examples:
mono printpoly2(pmono2) --> bern printpoly2(`pbern2')
mono printpoly2(qmono2) --> bern printpoly2(`qbern2')
Subprogram (2) examples:
p(0.25) = printnum(evalbern2(`pbern2',25))
p(7.50) = printnum(evalbern2(`pbern2',750))
q(0.25) = printnum(evalbern2(`qbern2',25))
q(7.50) = printnum(evalbern2(`qbern2',750))
Subprogram (3) examples:
mono printpoly3(pmono3) --> bern printpoly3(`pbern3')
mono printpoly3(qmono3) --> bern printpoly3(`qbern3')
mono printpoly3(rmono3) --> bern printpoly3(`rbern3')
Subprogram (4) examples:
p(0.25) = printnum(evalbern3(`pbern3',25))
p(7.50) = printnum(evalbern3(`pbern3',750))
q(0.25) = printnum(evalbern3(`qbern3',25)) <-- rounding error
q(7.50) = printnum(evalbern3(`qbern3',750)) <-- rounding error
r(0.25) = printnum(evalbern3(`rbern3',25)) <-- rounding error
r(7.50) = printnum(evalbern3(`rbern3',750)) <-- overflow
Subprogram (5) examples:
printpoly2(`pbern2') --> printpoly3(`pbern3a')
printpoly2(`qbern2') --> printpoly3(`qbern3a')- Output:
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 p(7.50) = 1.0 q(0.25) = 1.6875 q(7.50) = 184.7500 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.6666, 3.3333, 6.0) mono (1.0, 2.0, 3.0, 4.0) --> bern (1.0, 1.6666, 3.3333, 10.0) Subprogram (4) examples: p(0.25) = 1.0 p(7.50) = 1.0 q(0.25) = 1.6872 <-- rounding error q(7.50) = 184.7306 <-- rounding error r(0.25) = 1.7497 <-- rounding error r(7.50) = -2422.-7366 <-- overflow Subprogram (5) examples: (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)
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
- Output:
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)
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 m_to_bern2(sequence monomial_coefficients) atom {a0,a1,a2} = monomial_coefficients, a01 = a0 + a1/2, a02 = a0 + a1 + a2 return {a0, a01, a02} end function function m_to_bern3(sequence monomial_coefficients) atom {a0, a1, a2, a3} = monomial_coefficients, a01 = a0 + a1/3, a02 = a0 + a1/3*2 + a2/3, a03 = a0 + a1 + a2 + a3 return {a0, a01, a02, a03} end function function bern2_to_bern3(sequence bernstein_coefficients) atom {b0, b1, b2} = bernstein_coefficients, b01 = b0/3 + b1*2/3, b12 = b1*2/3 + b2/3 return {b0, b01, b12, b2} end function function eval_bern2(sequence bernstein_coefficients, atom t) // using de Casteljau’s algorithm atom {b0,b1,b2} = bernstein_coefficients, s = 1-t, // first mid-points (obvs "mid"dle only when t=0.5) b01 = s*b0 + t*b1, b12 = s*b1 + t*b2, // second mid-point is on the curve b012 = s*b01 + t*b12 return b012 end function function eval_bern3(sequence bernstein_coefficients, atom t) // using de Casteljau’s algorithm atom {b0, b1, b2, b3} = bernstein_coefficients, s = 1-t, // first mid-points (ditto) b01 = s*b0 + t*b1, b12 = s*b1 + t*b2, b23 = s*b2 + t*b3, // second mid-points b012 = s*b01 + t*b12, b123 = s*b12 + t*b23, // third mid-point is on the curve b0123 = s*b012 + t*b123 return b0123 end function function eval_m2(sequence monomial_coefficients, atom t) atom {a0, a1, a2} = monomial_coefficients return a0 + t*(a1 + t*a2) -- Horner’s rule end function function eval_m3(sequence monomial_coefficients, atom t) atom {a0, a1, a2, a3} = monomial_coefficients return a0 + t*(a1 + t*(a2 + t*a3)) -- Horner’s rule end function constant pm2 = {1, 0, 0}, -- p(x) = 1 pm3 = {1, 0, 0, 0}, -- (ditto in degree 3) qm2 = {1, 2, 3}, -- q(x) = 1 + 2x + 3x² qm3 = {1, 2, 3, 0}, -- (ditto in degree 3) rm3 = {1, 2, 3, 4}, -- r(x) = 1 + 2x + 3x² + 4x³ pb2 = m_to_bern2(pm2), pb3 = m_to_bern3(pm3), qb2 = m_to_bern2(qm2), qb3 = m_to_bern3(qm3), rb3 = m_to_bern3(rm3) printf(1," m_to_bern2(%v) --> %v\n", {pm2,pb2}) printf(1," m_to_bern2(%v) --> %v\n", {qm2,qb2}) printf(1," m_to_bern3(%v) --> %v\n", {pm3,pb3}) printf(1," m_to_bern3(%v) --> %v\n", {qm3,qb3}) printf(1," m_to_bern3(%v) --> %v\n", {rm3,rb3}) printf(1," bern2_to_bern3(%v) --> %v\n", {pb2,bern2_to_bern3(pb2)}) printf(1," bern2_to_bern3(%v) --> %v\n", {qb2,bern2_to_bern3(qb2)}) printf(1,"Evaluating bernstein degree 2 examples:\n") for x in {0.25, 7.50} do printf(1," p(%.2f) = %8g (mono: %g)\n",{x,eval_bern2(pb2,x),eval_m2(pm2,x)}) printf(1," q(%.2f) = %8g (mono: %g)\n",{x,eval_bern2(qb2,x),eval_m2(qm2,x)}) end for printf(1,"Evaluating bernstein degree 3 examples:\n") for x in {0.25, 7.50} do printf(1," p(%.2f) = %8g (mono: %g)\n",{x,eval_bern3(pb3,x),eval_m3(pm3,x)}) printf(1," q(%.2f) = %8g (mono: %g)\n",{x,eval_bern3(qb3,x),eval_m3(qm3,x)}) printf(1," r(%.2f) = %8g (mono: %g)\n",{x,eval_bern3(rb3,x),eval_m3(rm3,x)}) end for
- Output:
m_to_bern2({1,0,0}) --> {1,1,1}
m_to_bern2({1,2,3}) --> {1,2,6}
m_to_bern3({1,0,0,0}) --> {1,1,1,1}
m_to_bern3({1,2,3,0}) --> {1,1.666666667,3.333333333,6}
m_to_bern3({1,2,3,4}) --> {1,1.666666667,3.333333333,10}
bern2_to_bern3({1,1,1}) --> {1,1,1,1}
bern2_to_bern3({1,2,6}) --> {1,1.666666667,3.333333333,6}
Evaluating bernstein degree 2 examples:
p(0.25) = 1 (mono: 1)
q(0.25) = 1.6875 (mono: 1.6875)
p(7.50) = 1 (mono: 1)
q(7.50) = 184.75 (mono: 184.75)
Evaluating bernstein degree 3 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.50) = 1 (mono: 1)
q(7.50) = 184.75 (mono: 184.75)
r(7.50) = 1872.25 (mono: 1872.25)
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)
Raku
# 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]
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)