Rosetta Code

Bernstein basis polynomials: Difference between revisions

Page Discussion

Bernstein basis polynomials (view source)

Revision as of 09:26, 8 February 2024

28,229 bytes added ,  3 months ago
Add C# implementation
(→‎{{header|ALGOL 68}}: Simplify a little using features of Algol 68 not in Algol 60)
(Add C# implementation)
 
(12 intermediate revisions by 6 users not shown)
Line 1:
{{draft task}}
The <math>n + 1</math> [[wp:Bernstein_polynomial|''Bernstein basis polynomials'']] of degree <math>n</math> are defined as
:<math>b_{k,n}(t) = \binom{n}{k} t^{k} \left( 1 - t \right)^{n - k},\quad k = 0,\ldots,n</math>
Line 333:
 
=={{header|ALGOL 68}}==
Algol 68 doesn't have procedure overloading but does have monadic and dyadic operator overloading and user defined operator symbols.
<syntaxhighlight lang="algol68">
BEGIN # Bernstein Basis Polynomials - translated from the Algol 60 #
Line 391 ⟶ 392:
# coefficients q0, q1, q2 to the cubic Bernstein #
# coefficients c0, c1, c2, c3; #
PROCOP bern2to3 TOBERNTHREE = ( BERNTWO q, REF REAL c0, c1, c2, c3 )VOIDBERNTHREE:
BEGIN
c0HEAP :=BERNTHREE b0 OF qc;
c1 := ((1/3) * b0 OF q)c +:= ((2/3) * b1b0 OF q);
c2b1 OF c := ((21/3) * b1b0 OF q) + ((12/3) * b2b1 OF q);
c3b2 OF c := ((2/3) * b1 OF q) + ((1/3) * b2 OF q);
END # bern2to3 #b3 OF c := b2 OF q;
c
END # TOBERNTHREE # ;
 
BEGIN
Line 426 ⟶ 429:
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 );
Line 433 ⟶ 452:
BERNTHREE q3 := BERNTHREE( "q", 0, 0, 0, 0 );
BERNTHREE r3 := BERNTHREE( "r", 0, 0, 0, 0 );
 
REAL pc0 := 0, pc1 := 0, pc2 := 0, pc3 := 0;
REAL qc0 := 0, qc1 := 0, qc2 := 0, qc3 := 0;
 
REAL p0m = 1, p1m = 0, p2m = 0;
Line 444 ⟶ 460:
tobern2( q0m, q1m, q2m, q2 );
print( ( "Subprogram (1) examples:", newline ) );
print( ( " mono ( ", TOSTRING f[]REAL( p0m, p1m, p2m ), ", ",--> f( p1m ),bern ", "TOSTRING p2, f( p2mnewline ) );
print( ( " mono ", " ) --> bernTOSTRING []REAL( "q0m, f(q1m, b0 OF p2q2m ), ", ",--> f(bern b1 OF p2 )", ",TOSTRING "q2, f(newline b2 OF p2) );
, " )", newline
)
);
print( ( " mono ( ", f( q0m ), ", ", f( q1m ), ", ", f( q2m )
, " ) --> bern ( ", f( b0 OF q2 ), ", ", f( b1 OF q2 ), ", ", f( b2 OF q2 )
, " )", newline
)
);
 
print( ( "Subprogram (2) examples:", newline ) );
Line 465 ⟶ 473:
tobern3( r0m, r1m, r2m, r3m, r3 );
print( ( "Subprogram (3) examples:", newline ) );
print( ( " mono ( ", TOSTRING f[]REAL( p0m ), "p1m, "p2m, f(0 p1m ), ", ",--> f(bern p2m )", ",TOSTRING "p3, f(newline 0) );
print( ,( " ) --> bern (mono ", fTOSTRING []REAL( b0 OF p3 )q0m, "q1m, "q2m, f(0 b1 OF p3 ), ", ",--> f( b2 OF p3 ),bern ", "TOSTRING q3, f(newline b3 OF p3) );
print( ( " mono ", TOSTRING []REAL( r0m, "r1m, r2m, r3m ), " --> bern ", TOSTRING r3, newline ) );
)
);
print( ( " mono ( ", f( q0m ), ", ", f( q1m ), ", ", f( q2m ), ", ", f( 0 )
, " ) --> bern ( ", f( b0 OF q3 ), ", ", f( b1 OF q3 ), ", ", f( b2 OF q3 ), ", ", f( b3 OF q3 )
, " )", newline
)
);
print( ( " mono ( ", f( r0m ), ", ", f( r1m ), ", ", f( r2m ), ", ", f( r3m )
, " ) --> bern ( ", f( b0 OF r3 ), ", ", f( b1 OF r3 ), ", ", f( b2 OF r3 ), ", ", f( b3 OF r3 )
, " )", newline
)
);
 
print( ( "Subprogram (4) examples:", newline ) );
Line 489 ⟶ 485:
r3 SHOWEVAL 7.50;
 
bern2to3( p2, pc0, pc1, pc2, pc3 );
bern2to3( q2, qc0, qc1, qc2, qc3 );
print( ( "Subprogram (5) examples:", newline ) );
print( ( " bern ( ", f( b0 OFTOSTRING p2 ), ", --> bern ", f(TOSTRING b1 OFTOBERNTHREE p2 ), ",newline ", f( b2 OF p2) );
print( ( " bern ", TOSTRING q2, " ) --> bern ( ", f(TOSTRING pc0TOBERNTHREE )q2, ", ", f( pc1 newline ), ", ", f( pc2 ), ", ", f( pc3 )
 
, " )", newline
)
);
print( ( " bern ( ", f( b0 OF q2 ), ", ", f( b1 OF q2 ), ", ", f( b2 OF q2 )
, " ) --> bern ( ", f( qc0 ), ", ", f( qc1 ), ", ", f( qc2 ), ", ", f( qc3 )
, " )", newline
)
)
END
 
Line 1,097 ⟶ 1,084:
mono {1, 1, 1} --> bern {1, 1, 1, 1}
mono {1, 2, 6} --> bern {1, 1.6666666666667, 3.3333333333333, 6}
</pre>
 
=={{header|C#}}==
{{trans|Go}}
<syntaxhighlight lang="C#">
using System;
 
class Program
{
static double[] ToBern2(double[] a)
{
return new double[] { a[0], a[0] + a[1] / 2, a[0] + a[1] + a[2] };
}
 
static double EvalBern2(double[] b, double t)
{
double s = 1.0 - t;
double b01 = s * b[0] + t * b[1];
double b12 = s * b[1] + t * b[2];
return s * b01 + t * b12;
}
 
static double[] ToBern3(double[] a)
{
return new double[] { a[0], a[0] + a[1] / 3, a[0] + a[1] * 2 / 3 + a[2] / 3, a[0] + a[1] + a[2] + a[3] };
}
 
static double EvalBern3(double[] b, double t)
{
double s = 1.0 - t;
double b01 = s * b[0] + t * b[1];
double b12 = s * b[1] + t * b[2];
double b23 = s * b[2] + t * b[3];
double b012 = s * b01 + t * b12;
double b123 = s * b12 + t * b23;
return s * b012 + t * b123;
}
 
static double[] Bern2To3(double[] q)
{
return new double[] { q[0], q[0] / 3 + q[1] * 2 / 3, q[1] * 2 / 3 + q[2] / 3, q[2] };
}
 
static double EvalMono2(double[] a, double t)
{
return a[0] + (t * (a[1] + (t * a[2])));
}
 
static double EvalMono3(double[] a, double t)
{
return a[0] + (t * (a[1] + (t * (a[2] + (t * a[3])))));
}
 
static void Main(string[] args)
{
double[] pm = { 1, 0, 0 };
double[] qm = { 1, 2, 3 };
double[] rm = { 1, 2, 3, 4 };
double x, y, m;
 
Console.WriteLine("Subprogram(1) examples:");
var pb2 = ToBern2(pm);
var qb2 = ToBern2(qm);
Console.WriteLine($"mono [{string.Join(", ", pm)}] --> bern [{string.Join(", ", pb2)}]");
Console.WriteLine($"mono [{string.Join(", ", qm)}] --> bern [{string.Join(", ", qb2)}]");
 
Console.WriteLine("\nSubprogram(2) examples:");
x = 0.25;
y = EvalBern2(pb2, x);
m = EvalMono2(pm, x);
Console.WriteLine($"p({x:F2}) = {y:G14} (mono {m:G14})");
x = 7.5;
y = EvalBern2(pb2, x);
m = EvalMono2(pm, x);
Console.WriteLine($"p({x:F2}) = {y:G14} (mono {m:G14})");
 
x = 0.25;
y = EvalBern2(qb2, x);
m = EvalMono2(qm, x);
Console.WriteLine($"q({x:F2}) = {y:G14} (mono {m:G14})");
x = 7.5;
y = EvalBern2(qb2, x);
m = EvalMono2(qm, x);
Console.WriteLine($"q({x:F2}) = {y:G14} (mono {m:G14})");
 
Console.WriteLine("\nSubprogram(3) examples:");
var pb3 = ToBern3(new double[] { pm[0], pm[1], pm[2], 0 });
var qb3 = ToBern3(new double[] { qm[0], qm[1], qm[2], 0 });
var rb3 = ToBern3(rm);
Console.WriteLine($"mono [{string.Join(", ", pm)}] --> bern [{string.Join(", ", pb3)}]");
Console.WriteLine($"mono [{string.Join(", ", qm)}] --> bern [{string.Join(", ", qb3)}]");
Console.WriteLine($"mono [{string.Join(", ", rm)}] --> bern [{string.Join(", ", rb3)}]");
 
Console.WriteLine("\nSubprogram(4) examples:");
x = 0.25;
y = EvalBern3(pb3, x);
m = EvalMono3(new double[] { pm[0], pm[1], pm[2], 0 }, x);
Console.WriteLine($"p({x:F2}) = {y:G14} (mono {m:G14})");
x = 7.5;
y = EvalBern3(pb3, x);
m = EvalMono3(new double[] { pm[0], pm[1], pm[2], 0 }, x);
Console.WriteLine($"p({x:F2}) = {y:G14} (mono {m:G14})");
 
x = 0.25;
y = EvalBern3(qb3, x);
m = EvalMono3(new double[] { qm[0], qm[1], qm[2], 0 }, x);
Console.WriteLine($"q({x:F2}) = {y:G14} (mono {m:G14})");
x = 7.5;
y = EvalBern3(qb3, x);
m = EvalMono3(new double[] { qm[0], qm[1], qm[2], 0 }, x);
Console.WriteLine($"q({x:F2}) = {y:G14} (mono {m:G14})");
 
x = 0.25;
y = EvalBern3(rb3, x);
m = EvalMono3(rm, x);
Console.WriteLine($"r({x:F2}) = {y:G14} (mono {m:G14})");
x = 7.5;
y = EvalBern3(rb3, x);
m = EvalMono3(rm, x);
Console.WriteLine($"r({x:F2}) = {y:G14} (mono {m:G14})");
 
Console.WriteLine("\nSubprogram(5) examples:");
var pc = Bern2To3(pb2);
var qc = Bern2To3(qb2);
Console.WriteLine($"bern [{string.Join(", ", pb2)}] --> bern3 [{string.Join(", ", pc)}]");
Console.WriteLine($"bern [{string.Join(", ", qb2)}] --> bern3 [{string.Join(", ", qc)}]");
}
}
</syntaxhighlight>
{{out}}
<pre>
Subprogram(1) examples:
mono [1, 0, 0] --> bern [1, 1, 1]
mono [1, 2, 3] --> bern [1, 2, 6]
 
Subprogram(2) examples:
p(0.25) = 1 (mono 1)
p(7.50) = 1 (mono 1)
q(0.25) = 1.6875 (mono 1.6875)
q(7.50) = 184.75 (mono 184.75)
 
Subprogram(3) examples:
mono [1, 0, 0] --> bern [1, 1, 1, 1]
mono [1, 2, 3] --> bern [1, 1.66666666666667, 3.33333333333333, 6]
mono [1, 2, 3, 4] --> bern [1, 1.66666666666667, 3.33333333333333, 10]
 
Subprogram(4) examples:
p(0.25) = 1 (mono 1)
p(7.50) = 1 (mono 1)
q(0.25) = 1.6875 (mono 1.6875)
q(7.50) = 184.75 (mono 184.75)
r(0.25) = 1.75 (mono 1.75)
r(7.50) = 1872.25 (mono 1872.25)
 
Subprogram(5) examples:
bern [1, 1, 1] --> bern3 [1, 1, 1, 1]
bern [1, 2, 6] --> bern3 [1, 1.66666666666667, 3.33333333333333, 6]
 
</pre>
 
=={{header|C++}}==
<syntaxhighlight lang="c++">
#include <cstdint>
#include <iostream>
#include <vector>
 
std::string to_string(const std::vector<double>& list) {
std::string result = "[";
for ( uint64_t i = 0; i < list.size() - 1; ++i ) {
result += std::to_string(list[i]) + ", ";
}
result += std::to_string(list.back()) + "]";
return result;
}
 
// Subprogram (1)
std::vector<double> monomial_to_bernstein_degree2(const std::vector<double>& monomial) {
return std::vector<double>{ monomial[0],
monomial[0] + ( monomial[1] / 2.0 ),
monomial[0] + monomial[1] + monomial[2] };
}
 
// Subprogram (2)
double evaluate_bernstein_degree2(const std::vector<double>& bernstein, const double& t) {
// de Casteljau’s algorithm
const double s = 1 - t;
const double b01 = ( s * bernstein[0] ) + ( t * bernstein[1] );
const double b12 = ( s * bernstein[1] ) + ( t * bernstein[2] );
return ( s * b01 ) + ( t * b12 );
}
 
// Subprogram (3)
std::vector<double> monomial_to_bernstein_degree3(const std::vector<double>& monomial) {
return std::vector<double>{ monomial[0],
monomial[0] + ( monomial[1] / 3.0 ),
monomial[0] + ( 2.0 * monomial[1] / 3.0 ) + ( monomial[2] / 3.0 ),
monomial[0] + monomial[1] + monomial[2] + monomial[3] };
}
 
// Subprogram (4)
double evaluate_bernstein_degree3(const std::vector<double>& bernstein, const double& t) {
// de Casteljau’s algorithm
const double s = 1 - t;
const double b01 = ( s * bernstein[0] ) + ( t * bernstein[1] );
const double b12 = ( s * bernstein[1] ) + ( t * bernstein[2] );
const double b23 = ( s * bernstein[2] ) + ( t * bernstein[3] );
const double b012 = ( s * b01 ) + ( t * b12 );
const double b123 = ( s * b12 ) + ( t * b23 );
return ( s * b012 ) + ( t * b123 );
}
 
// Subprogram (5)
std::vector<double> bernstein_degree2_to_degree3(const std::vector<double>& bernstein) {
return std::vector<double>{ bernstein[0],
( bernstein[0] / 3.0 ) + ( 2.0 * bernstein[1] / 3.0 ),
( 2.0 * bernstein[1] / 3.0 ) + ( bernstein[2] / 3.0 ),
bernstein[2] };
}
 
double evaluate_monomial_degree2(const std::vector<double>& monomial, const double& t) {
// Horner’s rule
return monomial[0] + ( t * ( monomial[1] + ( t * monomial[2] ) ) );
}
 
double evaluate_monomial_degree3(const std::vector<double>& monomial, const double& t) {
// Horner’s rule
return monomial[0] + ( t * ( monomial[1] + ( t * ( monomial[2] + ( t * monomial[3] ) ) ) ) );
}
 
int main() {
/**
* For the following polynomials, use Subprogram (1) to find coefficients in the degree-2 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
*/
std::vector<double> pMonomial2 = { 1.0, 0.0, 0.0 };
std::vector<double> qMonomial2 = { 1.0, 2.0, 3.0 };
std::vector<double> pBernstein2 = monomial_to_bernstein_degree2(pMonomial2);
std::vector<double> qBernstein2 = monomial_to_bernstein_degree2(qMonomial2);
std::cout << "Subprogram (1) examples:" << std::endl;
std::cout << " monomial " + to_string(pMonomial2) + " --> bernstein " + to_string(pBernstein2) << std::endl;
std::cout << " monomial " + to_string(qMonomial2) + " --> bernstein " + to_string(qBernstein2) << std::endl;
 
/**
* Use Subprogram (2) to evaluate p(x) and q(x) at x = 0.25, 7.50. Display the results.
* Optionally also display results from evaluating in the original monomial basis.
*/
std::cout << "Subprogram (2) examples:" << std::endl;
for ( const double& x : { 0.25, 7.50 } ) {
std::cout << " p(" << x << ") = " << evaluate_bernstein_degree2(pBernstein2, x)
<< " ( mono: " << evaluate_monomial_degree2(pMonomial2, x) << " )" << std::endl;
}
for ( const double& x : { 0.25, 7.50 } ) {
std::cout << " q(" << x << ") = " << evaluate_bernstein_degree2(qBernstein2, x)
<< " ( mono: " << evaluate_monomial_degree2(qMonomial2, x) << " )" << std::endl;
}
 
/**
* For the following polynomials, use Subprogram (3) to find coefficients in the degree-3 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
* r(x) = 1 + 2x + 3x² + 4x³
*
* Display the results.
*/
std::vector<double> pMonomial3 = { 1.0, 0.0, 0.0, 0.0 };
std::vector<double> qMonomial3 = { 1.0, 2.0, 3.0, 0.0 };
std::vector<double> rMonomial3 = { 1.0, 2.0, 3.0, 4.0 };
std::vector<double> pBernstein3 = monomial_to_bernstein_degree3(pMonomial3);
std::vector<double> qBernstein3 = monomial_to_bernstein_degree3(qMonomial3);
std::vector<double> rBernstein3 = monomial_to_bernstein_degree3(rMonomial3);
std::cout << "Subprogram (3) examples:" << std::endl;
std::cout << " monomial " + to_string(pMonomial3) + " --> bernstein " + to_string(pBernstein3) << std::endl;
std::cout << " monomial " + to_string(qMonomial3) + " --> bernstein " + to_string(qBernstein3) << std::endl;
std::cout << " monomial " + to_string(rMonomial3) + " --> bernstein " + to_string(rBernstein3) << std::endl;
 
/**
* Use Subprogram (4) to evaluate p(x), q(x), and r(x) at x = 0.25, 7.50. Display the results.
* Optionally also display results from evaluating in the original monomial basis.
*/
std::cout << "Subprogram (4) examples:" << std::endl;
for ( const double& x : { 0.25, 7.50 } ) {
std::cout << " p(" << x << ") = " << evaluate_bernstein_degree3(pBernstein3, x)
<< " ( mono: " << evaluate_monomial_degree3(pMonomial3, x) << " )" << std::endl;
}
for ( const double& x : { 0.25, 7.50 } ) {
std::cout << " q(" << x << ") = " << evaluate_bernstein_degree3(qBernstein3, x)
<< " ( mono: " << evaluate_monomial_degree3(qMonomial3, x) << " )" << std::endl;
}
for ( const double& x : { 0.25, 7.50 } ) {
std::cout << " r(" << x << ") = " << evaluate_bernstein_degree3(rBernstein3, x)
<< " ( mono: " << evaluate_monomial_degree3(rMonomial3, x) << " )" << std::endl;
}
 
/**
* For the following polynomials, using the result of Subprogram (1) applied to the polynomial,
* use Subprogram (5) to get coefficients for the degree-3 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
*
* Display the results.
*/
std::cout << "Subprogram (5) examples:" << std::endl;
std::vector<double> pBernstein3a = bernstein_degree2_to_degree3(pBernstein2);
std::vector<double> qBernstein3a = bernstein_degree2_to_degree3(qBernstein2);
std::cout << " bernstein " + to_string(pBernstein2) + " --> bernstein " + to_string(pBernstein3a) << std::endl;
std::cout << " bernstein " + to_string(qBernstein2) + " --> bernstein " + to_string(qBernstein3a) << std::endl;
}
</syntaxhighlight>
{{ out }}
<pre>
Subprogram (1) examples:
monomial [1.000000, 0.000000, 0.000000] --> bernstein [1.000000, 1.000000, 1.000000]
monomial [1.000000, 2.000000, 3.000000] --> bernstein [1.000000, 2.000000, 6.000000]
Subprogram (2) examples:
p(0.25) = 1 ( mono: 1 )
p(7.5) = 1 ( mono: 1 )
q(0.25) = 1.6875 ( mono: 1.6875 )
q(7.5) = 184.75 ( mono: 184.75 )
Subprogram (3) examples:
monomial [1.000000, 0.000000, 0.000000, 0.000000] --> bernstein [1.000000, 1.000000, 1.000000, 1.000000]
monomial [1.000000, 2.000000, 3.000000, 0.000000] --> bernstein [1.000000, 1.666667, 3.333333, 6.000000]
monomial [1.000000, 2.000000, 3.000000, 4.000000] --> bernstein [1.000000, 1.666667, 3.333333, 10.000000]
Subprogram (4) examples:
p(0.25) = 1 ( mono: 1 )
p(7.5) = 1 ( mono: 1 )
q(0.25) = 1.6875 ( mono: 1.6875 )
q(7.5) = 184.75 ( mono: 184.75 )
r(0.25) = 1.75 ( mono: 1.75 )
r(7.5) = 1872.25 ( mono: 1872.25 )
Subprogram (5) examples:
bernstein [1.000000, 1.000000, 1.000000] --> bernstein [1.000000, 1.000000, 1.000000, 1.000000]
bernstein [1.000000, 2.000000, 6.000000] --> bernstein [1.000000, 1.666667, 3.333333, 6.000000]
</pre>
 
Line 1,420 ⟶ 1,743:
mono [1 1 1] --> bern [1 1 1 1]
mono [1 2 6] --> bern [1 1.6666666666667 3.3333333333333 6]
</pre>
 
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.util.Arrays;
import java.util.List;
 
public final class BernsteinBasisPolynomials {
 
public static void main(String[] args) {
/**
* For the following polynomials, use Subprogram (1) to find coefficients in the degree-2 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
*/
QuadraticCoefficients pMonomial2 = new QuadraticCoefficients(1.0, 0.0, 0.0);
QuadraticCoefficients qMonomial2 = new QuadraticCoefficients(1.0, 2.0, 3.0);
QuadraticCoefficients pBernstein2 = monomialToBernsteinDegree2(pMonomial2);
QuadraticCoefficients qBernstein2 = monomialToBernsteinDegree2(qMonomial2);
System.out.println("Subprogram (1) examples:");
System.out.println(" monomial " + pMonomial2 + " --> bernstein " + pBernstein2);
System.out.println(" monomial " + qMonomial2 + " --> bernstein " + qBernstein2);
/**
* Use Subprogram (2) to evaluate p(x) and q(x) at x = 0.25, 7.50. Display the results.
* Optionally also display results from evaluating in the original monomial basis.
*/
System.out.println("Subprogram (2) examples:");
for ( double x : List.of( 0.25, 7.50 ) ) {
System.out.println(" p(" + x + ") = " + evaluateBernsteinDegree2(pBernstein2, x) +
" ( mono: " + evaluateMonomialDegree2(pMonomial2, x) + " )");
}
for ( double x : List.of( 0.25, 7.50 ) ) {
System.out.println(" q(" + x + ") = " + evaluateBernsteinDegree2(qBernstein2, x) +
" ( mono: " + evaluateMonomialDegree2(qMonomial2, x) + " )");
}
/**
* For the following polynomials, use Subprogram (3) to find coefficients in the degree-3 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
* r(x) = 1 + 2x + 3x² + 4x³
*
* Display the results.
*/
CubicCoefficients pMonomial3 = new CubicCoefficients(1.0, 0.0, 0.0, 0.0);
CubicCoefficients qMonomial3 = new CubicCoefficients(1.0, 2.0, 3.0, 0.0);
CubicCoefficients rMonomial3 = new CubicCoefficients(1.0, 2.0, 3.0, 4.0);
CubicCoefficients pBernstein3 = monomialToBernsteinDegree3(pMonomial3);
CubicCoefficients qBernstein3 = monomialToBernsteinDegree3(qMonomial3);
CubicCoefficients rBernstein3 = monomialToBernsteinDegree3(rMonomial3);
System.out.println("Subprogram (3) examples:");
System.out.println(" monomial " + pMonomial3 + " --> bernstein " + pBernstein3);
System.out.println(" monomial " + qMonomial3 + " --> bernstein " + qBernstein3);
System.out.println(" monomial " + rMonomial3 + " --> bernstein " + rBernstein3);
/**
* Use Subprogram (4) to evaluate p(x), q(x), and r(x) at x = 0.25, 7.50. Display the results.
* Optionally also display results from evaluating in the original monomial basis.
*/
System.out.println("Subprogram (4) examples:");
for ( double x : List.of( 0.25, 7.50 ) ) {
System.out.println(" p(" + x + ") = " + evaluateBernsteinDegree3(pBernstein3, x) +
" ( mono: " + evaluateMonomialDegree3(pMonomial3, x) + " )");
}
for ( double x : List.of( 0.25, 7.50 ) ) {
System.out.println(" q(" + x + ") = " + evaluateBernsteinDegree3(qBernstein3, x) +
" ( mono: " + evaluateMonomialDegree3(qMonomial3, x) + " )");
}
for ( double x : List.of( 0.25, 7.50 ) ) {
System.out.println(" r(" + x + ") = " + evaluateBernsteinDegree3(rBernstein3, x) +
" ( mono: " + evaluateMonomialDegree3(rMonomial3, x) + " )");
}
/**
* For the following polynomials, using the result of Subprogram (1) applied to the polynomial,
* use Subprogram (5) to get coefficients for the degree-3 Bernstein basis:
*
* p(x) = 1
* q(x) = 1 + 2x + 3x²
*
* Display the results.
*/
System.out.println("Subprogram (5) examples:");
CubicCoefficients pBernstein3a = bernsteinDegree2ToDegree3(pBernstein2);
CubicCoefficients qBernstein3a = bernsteinDegree2ToDegree3(qBernstein2);
System.out.println(" bernstein " + pBernstein2 + " --> bernstein " + pBernstein3a);
System.out.println(" bernstein " + qBernstein2 + " --> bernstein " + qBernstein3a);
}
private static record QuadraticCoefficients(double q0, double q1, double q2) {
@Override
public String toString() {
return Arrays.asList(q0, q1, q2).toString();
}
}
private static record CubicCoefficients(double c0, double c1, double c2, double c3) {
@Override
public String toString() {
return Arrays.asList(c0, c1, c2, c3).toString();
}
}
// Subprogram (1)
private static QuadraticCoefficients monomialToBernsteinDegree2(QuadraticCoefficients monomial) {
return new QuadraticCoefficients( monomial.q0,
monomial.q0 + ( monomial.q1 / 2.0 ),
monomial.q0 + monomial.q1 + monomial.q2 );
}
// Subprogram (2)
private static double evaluateBernsteinDegree2(QuadraticCoefficients bernstein, double t) {
// de Casteljau’s algorithm
final double s = 1 - t;
final double b01 = ( s * bernstein.q0 ) + ( t * bernstein.q1 );
final double b12 = ( s * bernstein.q1 ) + ( t * bernstein.q2 );
return ( s * b01 ) + ( t * b12 );
}
// Subprogram (3)
private static CubicCoefficients monomialToBernsteinDegree3(CubicCoefficients monomial) {
return new CubicCoefficients( monomial.c0,
monomial.c0 + ( monomial.c1 / 3.0 ),
monomial.c0 + ( 2.0 * monomial.c1 / 3.0 ) + ( monomial.c2 / 3.0 ),
monomial.c0 + monomial.c1 + monomial.c2 + monomial.c3 );
}
// Subprogram (4)
private static double evaluateBernsteinDegree3(CubicCoefficients bernstein, double t) {
// de Casteljau’s algorithm
final double s = 1 - t;
final double b01 = ( s * bernstein.c0 ) + ( t * bernstein.c1 );
final double b12 = ( s * bernstein.c1 ) + ( t * bernstein.c2 );
final double b23 = ( s * bernstein.c2 ) + ( t * bernstein.c3 );
final double b012 = ( s * b01 ) + ( t * b12 );
final double b123 = ( s * b12 ) + ( t * b23 );
return ( s * b012 ) + ( t * b123 );
}
// Subprogram (5)
private static CubicCoefficients bernsteinDegree2ToDegree3(QuadraticCoefficients bernstein) {
return new CubicCoefficients( bernstein.q0,
( bernstein.q0 / 3.0 ) + ( 2.0 * bernstein.q1 / 3.0 ),
( 2.0 * bernstein.q1 / 3.0 ) + ( bernstein.q2 / 3.0 ),
bernstein.q2 );
}
private static double evaluateMonomialDegree2(QuadraticCoefficients monomial, double t) {
// Horner’s rule
return monomial.q0 + ( t * ( monomial.q1 + ( t * monomial.q2 ) ) );
}
private static double evaluateMonomialDegree3(CubicCoefficients monomial, double t) {
// Horner’s rule
return monomial.c0 + ( t * ( monomial.c1 + ( t * ( monomial.c2 + ( t * monomial.c3 ) ) ) ) );
}
 
}
</syntaxhighlight>
{{ out }}
<pre>
Subprogram (1) examples:
monomial [1.0, 0.0, 0.0] --> bernstein [1.0, 1.0, 1.0]
monomial [1.0, 2.0, 3.0] --> bernstein [1.0, 2.0, 6.0]
Subprogram (2) examples:
p(0.25) = 1.0 ( mono: 1.0 )
p(7.5) = 1.0 ( mono: 1.0 )
q(0.25) = 1.6875 ( mono: 1.6875 )
q(7.5) = 184.75 ( mono: 184.75 )
Subprogram (3) examples:
monomial [1.0, 0.0, 0.0, 0.0] --> bernstein [1.0, 1.0, 1.0, 1.0]
monomial [1.0, 2.0, 3.0, 0.0] --> bernstein [1.0, 1.6666666666666665, 3.333333333333333, 6.0]
monomial [1.0, 2.0, 3.0, 4.0] --> bernstein [1.0, 1.6666666666666665, 3.333333333333333, 10.0]
Subprogram (4) examples:
p(0.25) = 1.0 ( mono: 1.0 )
p(7.5) = 1.0 ( mono: 1.0 )
q(0.25) = 1.6874999999999998 ( mono: 1.6875 )
q(7.5) = 184.75000000000034 ( mono: 184.75 )
r(0.25) = 1.7499999999999998 ( mono: 1.75 )
r(7.5) = 1872.25 ( mono: 1872.25 )
Subprogram (5) examples:
bernstein [1.0, 1.0, 1.0] --> bernstein [1.0, 1.0, 1.0, 1.0]
bernstein [1.0, 2.0, 6.0] --> bernstein [1.0, 1.6666666666666665, 3.333333333333333, 6.0]
</pre>
 
Line 1,827 ⟶ 2,340:
(1.0, 1.0, 1.0) --> (1.0, 1.0, 1.0, 1.0)
(1.0, 2.0, 6.0) --> (1.0, 1.6667, 3.3332, 6.0)</pre>
 
=={{header|Nim}}==
{{trans|Python}}
<syntaxhighlight lang="Nim">type
Coeffs2 = (float, float, float)
Coeffs3 = (float, float, float, float)
 
# Subprogram (1).
func monomialToBernsteinDegree2(monomialCoefficients: Coeffs2): Coeffs2 =
let (a0, a1, a2) = monomialCoefficients
result = (a0, a0 + ((1/2) * a1), a0 + a1 + a2)
 
# Subprogram (2).
func evaluateBernsteinDegree2(bernsteinCoefficients: Coeffs2; t: float): float =
let (b0, b1, b2) = bernsteinCoefficients
# de Casteljau’s algorithm.
let s = 1 - t
let b01 = (s * b0) + (t * b1)
let b12 = (s * b1) + (t * b2)
result = (s * b01) + (t * b12)
 
# Subprogram (3).
func monomialToBernsteinDegree3(monomialCoefficients: Coeffs3): Coeffs3 =
let (a0, a1, a2, a3) = monomialCoefficients
result = (a0, a0 + ((1/3) * a1), a0 + ((2/3) * a1) + ((1/3) * a2), a0 + a1 + a2 + a3)
 
# Subprogram (4).
func evaluateBernsteinDegree3(bernsteinCoefficients: Coeffs3; t: float): float =
let (b0, b1, b2, b3) = bernsteinCoefficients
# de Casteljau’s algorithm.
let s = 1 - t
let b01 = (s * b0) + (t * b1)
let b12 = (s * b1) + (t * b2)
let b23 = (s * b2) + (t * b3)
let b012 = (s * b01) + (t * b12)
let b123 = (s * b12) + (t * b23)
result = (s * b012) + (t * b123)
 
# Subprogram (5).
func bernsteinDegree2ToDegree3(bernsteinCoefficients: Coeffs2): Coeffs3 =
let (b0, b1, b2) = bernstein_coefficients
result = (b0, ((1/3) * b0) + ((2/3) * b1), ((2/3) * b1) + ((1/3) * b2), b2)
 
func evaluateMonomialDegree2(monomialCoefficients: Coeffs2; t: float): float =
let (a0, a1, a2) = monomialCoefficients
# Horner’s rule.
result = a0 + (t * (a1 + (t * a2)))
 
func evaluateMonomialDegree3(monomialCoefficients: Coeffs3; t: float): float =
let (a0, a1, a2, a3) = monomialCoefficients
# Horner’s rule.
result = a0 + (t * (a1 + (t * (a2 + (t * a3)))))
 
#
# For the following polynomials, use Subprogram (1) to find
# coefficients in the degree-2 Bernstein basis:
#
# p(x) = 1
# q(x) = 1 + 2x + 3x²
#
# Display the results.
#
let pmono2 = (1.0, 0.0, 0.0)
let qmono2 = (1.0, 2.0, 3.0)
let pbern2 = monomialToBernsteinDegree2(pmono2)
let qbern2 = monomialToBernsteinDegree2(qmono2)
echo "Subprogram (1) examples:"
echo " mono ", pmono2, " --> bern ", pbern2
echo " mono ", qmono2, " --> bern ", qbern2
 
#
# Use Subprogram (2) to evaluate p(x) and q(x) at x = 0.25, 7.50.
# Display the results. Optionally also display results from evaluating
# in the original monomial basis.
#
echo "Subprogram (2) examples:"
for x in [0.25, 7.50]:
echo " p(", x, ") = ", evaluateBernsteinDegree2(pbern2, x),
" (mono: ", evaluateMonomialDegree2(pmono2, x), ')'
for x in [0.25, 7.50]:
echo " q(", x, ") = ", evaluateBernsteinDegree2(qbern2, x),
" (mono: ", evaluateMonomialDegree2(qmono2, x), ')'
 
#
# For the following polynomials, use Subprogram (3) to find
# coefficients in the degree-3 Bernstein basis:
#
# p(x) = 1
# q(x) = 1 + 2x + 3x²
# r(x) = 1 + 2x + 3x² + 4x³
#
# Display the results.
#
let pmono3 = (1.0, 0.0, 0.0, 0.0)
let qmono3 = (1.0, 2.0, 3.0, 0.0)
let rmono3 = (1.0, 2.0, 3.0, 4.0)
let pbern3 = monomialToBernsteinDegree3(pmono3)
let qbern3 = monomialToBernsteinDegree3(qmono3)
let rbern3 = monomialToBernsteinDegree3(rmono3)
echo "Subprogram (3) examples:"
echo " mono ", pmono3, " --> bern ", pbern3
echo " mono ", qmono3, " --> bern ", qbern3
echo " mono ", rmono3, " --> bern ", rbern3
 
#
# Use Subprogram (4) to evaluate p(x), q(x), and r(x) at x = 0.25,
# 7.50. Display the results. Optionally also display results from
# evaluating in the original monomial basis.
#
echo "Subprogram (4) examples:"
for x in [0.25, 7.50]:
echo " p(", x, ") = ", evaluateBernsteinDegree3(pbern3, x),
" (mono: ", evaluateMonomialDegree3(pmono3, x), ')'
for x in [0.25, 7.50]:
echo " q(", x, ") = ", evaluateBernsteinDegree3(qbern3, x),
" (mono: ", evaluateMonomialDegree3(qmono3, x), ')'
for x in [0.25, 7.50]:
echo " r(", x, ") = ", evaluateBernsteinDegree3(rbern3, x),
" (mono: ", evaluateMonomialDegree3(rmono3, x), ')'
 
#
# For the following polynomials, using the result of Subprogram (1)
# applied to the polynomial, use Subprogram (5) to get coefficients
# for the degree-3 Bernstein basis:
#
# p(x) = 1
# q(x) = 1 + 2x + 3x²
#
# Display the results.
#
echo "Subprogram (5) examples:"
let pbern3a = bernsteinDegree2ToDegree3(pbern2)
let qbern3a = bernsteinDegree2ToDegree3(qbern2)
echo " bern ", pbern2, "--> bern ", pbern3a
echo " bern ", qbern2, "--> bern ", qbern3a
</syntaxhighlight>
 
{{out}}
<pre>Subprogram (1) examples:
mono (1.0, 0.0, 0.0) --> bern (1.0, 1.0, 1.0)
mono (1.0, 2.0, 3.0) --> bern (1.0, 2.0, 6.0)
Subprogram (2) examples:
p(0.25) = 1.0 (mono: 1.0)
p(7.5) = 1.0 (mono: 1.0)
q(0.25) = 1.6875 (mono: 1.6875)
q(7.5) = 184.75 (mono: 184.75)
Subprogram (3) examples:
mono (1.0, 0.0, 0.0, 0.0) --> bern (1.0, 1.0, 1.0, 1.0)
mono (1.0, 2.0, 3.0, 0.0) --> bern (1.0, 1.666666666666667, 3.333333333333333, 6.0)
mono (1.0, 2.0, 3.0, 4.0) --> bern (1.0, 1.666666666666667, 3.333333333333333, 10.0)
Subprogram (4) examples:
p(0.25) = 1.0 (mono: 1.0)
p(7.5) = 1.0 (mono: 1.0)
q(0.25) = 1.6875 (mono: 1.6875)
q(7.5) = 184.7500000000003 (mono: 184.75)
r(0.25) = 1.75 (mono: 1.75)
r(7.5) = 1872.25 (mono: 1872.25)
Subprogram (5) examples:
bern (1.0, 1.0, 1.0)--> bern (1.0, 1.0, 1.0, 1.0)
bern (1.0, 2.0, 6.0)--> bern (1.0, 1.666666666666667, 3.333333333333333, 6.0)
</pre>
 
=={{header|ObjectIcon}}==
Line 2,338 ⟶ 3,012:
bern (1, 1.0, 1) --> bern (1, 1.0, 1.0, 1)
bern (1, 2.0, 6) --> bern (1, 1.6666666666666665, 3.333333333333333, 6)</pre>
 
=={{header|Raku}}==
{{trans|Wren}}
<syntaxhighlight lang="raku" line># 20230601 Raku programming solution
 
sub toBern2 { [ @_[0], @_[0]+@_[1]/2, @_[0..2].sum ] }
 
sub toBern3 (@a) { @a[0], @a[0]+@a[1]/3, @a[0]+@a[1]*2/3+@a[2]/3, @a[0..3].sum }
 
sub evalBern-N (@b, $t) { # uses de Casteljau's algorithm
my ($s, @bern) = 1 - $t, @b.Slip;
while ( @bern.elems > 2 ) {
@bern = @bern.rotor(2 => -1).map: { $s * .[0] + $t * .[1] };
}
return $s * @bern[0] + $t * @bern[1]
}
 
sub evaluations (@m, @b, $x) {
my $m = ([o] map { $_ + $x * * }, @m)(0); # Horner's rule
return "p({$x.fmt: '%.2f'}) = { evalBern-N @b, $x } (mono $m)";
}
 
sub bern2to3 (@a) { @a[0], @a[0]/3+@a[1]*2/3, @a[1]*2/3+@a[2]/3, @a[2] }
 
my (@pm,@qm,@rm) := ([1, 0, 0], [1, 2, 3], [1, 2, 3, 4]);
 
say "Subprogram(1) examples:";
 
my (@pb2,@qb2) := (@pm,@qm).map: { toBern2 $_ };
say "mono [{.[0]}] --> bern [{.[1]}]" for (@pm,@pb2,@qm,@qb2).rotor: 2;
 
say "\nSubprogram(2) examples:";
 
for (@pm,@pb2,@qm,@qb2).rotor(2) X (0.25,7.5) -> [[@m,@b], $x] {
say evaluations @m, @b, $x
}
 
say "\nSubprogram(3) examples:";
 
.push(0) for (@pm,@qm);
my (@pb3,@qb3,@rb3) := (@pm,@qm,@rm).map: { toBern3 $_ };
say "mono [{.[0]}] --> bern [{.[1]}]" for (@pm,@pb3,@qm,@qb3,@rm,@rb3).rotor: 2;
 
say "\nSubprogram(4) examples:";
 
for (@pm,@pb3,@qm,@qb3,@rm,@rb3).rotor(2) X (0.25,7.5) -> [[@m,@b], $x] {
say evaluations @m, @b, $x
}
 
say "\nSubprogram(5) examples:";
 
my (@pc,@qc) := (@pb2,@qb2).map: { bern2to3 $_ };
say "mono [{.[1]}] --> bern [{.[0]}]" for (@pc,@pb2,@qc,@qb2).rotor: 2;
</syntaxhighlight>
{{out}}
<pre>Subprogram(1) examples:
mono [1 0 0] --> bern [1 1 1]
mono [1 2 3] --> bern [1 2 6]
 
Subprogram(2) examples:
p(0.25) = 1 (mono 1)
p(7.50) = 1 (mono 1)
p(0.25) = 1.6875 (mono 1.6875)
p(7.50) = 184.75 (mono 184.75)
 
Subprogram(3) examples:
mono [1 0 0 0] --> bern [1 1 1 1]
mono [1 2 3 0] --> bern [1 1.666667 3.333333 6]
mono [1 2 3 4] --> bern [1 1.666667 3.333333 10]
 
Subprogram(4) examples:
p(0.25) = 1 (mono 1)
p(7.50) = 1 (mono 1)
p(0.25) = 1.6875 (mono 1.6875)
p(7.50) = 184.75 (mono 184.75)
p(0.25) = 1.75 (mono 1.75)
p(7.50) = 1872.25 (mono 1872.25)
 
Subprogram(5) examples:
mono [1 1 1] --> bern [1 1 1 1]
mono [1 2 6] --> bern [1 1.666667 3.333333 6]</pre>
 
=={{header|Wren}}==
Line 2,344 ⟶ 3,099:
{{libheader|Wren-math}}
Note that the library method, ''Math.evalPoly'', evaluates polynomials of any degree using Horner's rule but requires the coefficients to be presented in highest to lowest degree order.
<syntaxhighlight lang="ecmascriptwren">import "./fmt" for Fmt
import "./math" for Math
 
Aspen138
337

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/343374...358581"