Bernstein basis polynomials: Difference between revisions

Add C# implementation
m (Minor correction to code.)
(Add C# implementation)
 
(2 intermediate revisions by 2 users not shown)
Line 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 2,763 ⟶ 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
 
337

edits