SpecialPolynomials.jl

AbstractSpecialPolynomial{B,T,X}

An abstract type to distinguish the different polynomial types in this package.

The concrete types specify different bases for the space of polynomials of degree n or less.

This package includes:

• Several classical orthogonal polynomials.
• Newton and Lagrange interpolating polynomials
• Bernstein polynomials

As many of the methods for the base Polynomials class are directly coded if possible, but quite a few depend on conversion to the base Polynomial type (which uses the standard polynomial basis).

AbstractOrthogonalPolynomial{T,X}

This is an alias for polynomials with an orthogonal basis (AbstractOrthogonalBasis) specified.

These polynomials have several properties, including an accompanying inner product satisfying ⟨yᵢ, yⱼ⟩ = cᵢδᵢⱼ.

In addition to methods inherited from the underlying AbstractPolynomial type, orthogonal polynomial types may have methods weight_function, generating_function, leading_term, norm2, jacobi_matrix, and gauss_nodes_weights, though none are exported.

Subtypes of AbstractCOPBasis <: AbstractOrthogonalBasis utilize the fact that the basis polynomials satisfy

(ax² + bx + c)yᵢ''(x) + (dx+e)yᵢ'(x) + λᵢyᵢ(x) = 0 (or a discrete analogue)

where the structural relations are functions of a,b,c,d,e. These allow default definitions for polynomial evaluation, addition, multiplication, differentiation, integration, and conversion to and from the Polynomial type (the FallingFactorial type in the discrete c case),

A key structural relation is the three-term recursion, yᵢ₊₁ = (Aᵢx + Bᵢ)yᵢ - Cᵢyᵢ₋₁. For systems specified by a weight function, the values of Aᵢ, Bᵢ, and Cᵢ can be generated, yielding formulas for polynomial evaluation, addition, and conversion to the Polynomial type through evaluation.

AbstractCCOPBasis <:  AbstractCOPBasis

Following Koepf and Schmersau, a family y(x)=p_n(x)=k_x⋅x^n + ... for n ∈ {0, 1,…}, k_n ≠ 0 of polynomials is a family of classical continuous orthogonal polynomials if each is a solution of the differential equation

(a⋅x²+b⋅x+c) ⋅ yᵢ'' + (d⋅x + e) ⋅ yᵢ' + λᵢ⋅ yᵢ = 0.

A family is characterized, up to choice of leading term, by the 5 coefficients: a,b,c,d,e. Let σ = (a⋅x²+b⋅x+c), τ = (d⋅x + e).

From these 5 coefficients several structural equations are represented. For example the three-point recursion.

P₍ᵢ₊₁₎ = (Aᵢ⋅x + Bᵢ) * Pᵢ - Cᵢ * P₍ᵢ₋₁₎,

where Aᵢ,Bᵢ,Cᵢ can be represented in formulas involving just a,b,c,d,e and i.

Rearranging gives the structural equation:

x⋅p_n = [an, bn, cn] ⋅ [p_{n+1}, p_n, p_{n-1}] (Eqn (7))

The other structural equations are (equation references are from Koepf and Schmersau):

σ⋅p'_n = [αn, βn, γn] ⋅ [p_{n+1}, p_n, p_{n-1}] (Eqn (9), n ≥ 1)

p_n = [ân, b̂n, ĉn] ⋅ [p'_{n+1}, p'_n, p'_{n-1}] (Eqn (19))

x⋅p'_n = [αᴵn, βᴵn, γᴵn] ⋅ [p'_{n+1}, p'_n, p'_{n-1}] (Eqn (14))

Using (7), Clenshaw polynomial evaluation using the three point recursion is defined.

Using (19), expressions for derivatives are found.

Using (19), expressions for integration are found (p7).

Using their theorems 2,4, and 5, connection coefficients, C(n,m) satisfying P_n(x) = ∑ C(n,m) Q_m(x) (n ≥ 0, 0 ≤ m ≤ n) are found. These allow fallback definitions for convert(Polynomial,p), convert(P, p::Polynomial), convert(P{α…}, p::P(β…)) and through composition polynomial multiplication, p*q.

If non-monic versions are desired, then the leading term can be specified through kn() (which by default is defined by the method k1k0(P,i), the ratio of kᵢ₊₁/kᵢ). The @register_monic macro is useful for creating monic versions through method delegation from the common non-monic systems. Similarly, the @register_shifted macro is useful to provide shifted versions (cf. ShiftedLegendre).

Registering a system, defining an abcde method, and optionally defining k1k0 is usually sufficient to define a new system, though the general equations may need specializations when algebraic cancellation is required.

The defaults for evaluation and multiplication are generally an order of magnitude slower than a directly defined function. For some families this is done (e.g. Chebyshev,ChebyshevU, Hermite, Laguerre), but not all.

Koekoek and Swarttouw present an encyclopedia of formula characterizing families of orthogonal polynomials.

 AbstractCDOPBasis{T,X} <: AbstractCOPBasis{T,X}

Following Koepf and Schmersau, a family y(x)=p_n(x)=k_x⋅x^n + ... for n ∈ {0, 1,…}, k_n ≠ 0 of polynomials is a family of classical discrete orthogonal polynomials if it is a solution of a difference equation

(a⋅x²+b⋅x+c) ⋅ Δ∇y + (d⋅x + e) ⋅ ∇' + λᵢ⋅ y = 0,

where Δy(x) = y(x+1) - y(x) and ∇y(x) = y(x) - y(x-1).

A family is characterized by the 5 coefficients: a,b,c,d,e. Let σ = (a⋅x²+b⋅x+c), τ = (d⋅x + e).

As in the classical-continuous-orthogonal-polynomial case AbstractCCOPBasis, from these 5 values the coefficients in the there-point recursion, and other structural equations can be represented. These allow polynomial multiplication, integration, differentiation, conversion, etc. to be defined generically.

Koekoek and Swarttouw present an encyclopedia of formula characterizing families of orthogonal polynomials.

For example, on p29 they give formulas for Hahn polynomials through:

n(n+α+β+1)y(x) = B(x)y(x+1) -[B(x)+D(x)]y(x) + D(x)y(x-1), with explicit values for B and D. Reexpressing gives: BΔy(x) - D∇y(x) -λ y(x) = 0. From the rexpressed Eqn (4) for Koepf & Schemersau we have the identification: σ+τ = B; σ=D, so τ=B-D. From this a,b,c,d,e can be gleaned.

The above, is termed the eigenvalue equation (e.g. Goertz and Offner), as it can be reexpressed as

Δ(D(x)⋅ω(x)⋅∇yᵢ(x) = λᵢ⋅ω(x)⋅yᵢ(x)

Jacobi{α,  β, T}

Implements the Jacobi polynomials. These have weight function w(x) = (1-x)^α ⋅ (1+x)^β over the domain [-1,1]. Many orthogonal polynomial types are special cases. The parameters are specified to the constructors:

julia> using Polynomials, SpecialPolynomials

julia> p = Jacobi{-1/2, -1/2}([0,0,1])
Jacobi{-0.5,-0.5}(1⋅Jᵅᵝ₂(x))

julia> convert(Polynomial, p)
Polynomial(-0.375 + 0.75*x^2)

julia> monic(p) = (q=convert(Polynomial,p); q/q[end])
monic (generic function with 1 method)

julia> monic(p) ≈  monic(basis(Chebyshev, 2))
true

Special cases include Jacobi{α-1/2,α-1/2} = Gegenbauer{α}, Jacobi{-1/2,-1/2} = Chebyshev, Jacobi{1/2,1/2} = ChebyshevU, and Jacobi{0,0} = Legendre.

Gegenbauer{α, T <: Number}

The Gegenbauer polynomials have weight function (1-x^2)^(α-1/2) over the domain [-1,1]. The parameter α is specified in the constructor. These are also called the ultra-spherical polynomials. The Legendre polynomials are the specialization Gegenbauer{1/2}.

julia> using Polynomials, SpecialPolynomials

julia> p =  Gegenbauer{1/2}([1,2,3])
Gegenbauer{0.5}(1⋅Cᵅ₀(x) + 2⋅Cᵅ₁(x) + 3⋅Cᵅ₂(x))

julia> convert(Polynomial, p)
Polynomial(-0.5 + 2.0*x + 4.5*x^2)

Chebyshev{<:Number}(coeffs::AbstractVector, var=:x)

Chebyshev polynomial of the first kind. These have domain $(-1,1)$ and weight function $1/\sqrt{1-x^2}$.

A polynomial is constructed from its coefficients a, lowest order first, optionally in terms of the given variable x. x can be a character, symbol, or string.

Examples

julia> using Polynomials, SpecialPolynomials

julia> Chebyshev([1, 0, 3, 4])
Chebyshev(1⋅T₀(x) + 3⋅T₂(x) + 4⋅T₃(x))

julia> Chebyshev([1, 2, 3, 0], :s)
Chebyshev(1⋅T₀(s) + 2⋅T₁(s) + 3⋅T₂(s))

julia> one(Chebyshev)
Chebyshev(1.0⋅T₀(x))

ChebyshevU{T}

Implements the Chebyshev polynomials of the second kind. These have weight function w(x) = sqrt(1-x^2) over the domain [-1,1].

julia> using Polynomials, SpecialPolynomials

julia> p = ChebyshevU([1,2,3])
ChebyshevU(1⋅U₀(x) + 2⋅U₁(x) + 3⋅U₂(x))

julia> convert(Polynomial, p)
Polynomial(-2.0 + 4.0*x + 12.0*x^2)

julia> derivative(p)
ChebyshevU(4.0⋅U₀(x) + 12.0⋅U₁(x))

julia> roots(p)
2-element Vector{ComplexF64}:
 -0.6076252185107651 + 0.0im
 0.27429188517743175 + 0.0im

Legendre{T}

Implements the Legendre polynomials. These have weight function w(x) = 1 over the domain [-1,1].

julia> using Polynomials, SpecialPolynomials

julia> p = Legendre([1,2,3])
Legendre(1⋅P₀(x) + 2⋅P₁(x) + 3⋅P₂(x))

julia> convert(Polynomial, p)
Polynomial(-0.5 + 2.0*x + 4.5*x^2)

julia> p2m, p2m1 = basis.(Legendre, (8,9)) # evaluation P_{2m+k}(-1) =  (-1)^k
(Legendre(1.0⋅P₈(x)), Legendre(1.0⋅P₉(x)))

julia> p2m(-1) == 1
false

julia> p2m1(-1) == -1
false

julia> n = 5  # verify  Rodrigues' formula
5

julia> x = Polynomial(:x)
Polynomial(1.0*x)

julia> derivative((x^2-1)^n, n) - 2^n *  factorial(n) * basis(Legendre, n)
LaurentPolynomial(0.0)

julia> p4, p5  =  basis.(Legendre, (4,5)) # verify  orthogonality  of  P₄,P₅
(Legendre(1.0⋅P₄(x)), Legendre(1.0⋅P₅(x)))

julia> SpecialPolynomials.innerproduct(Legendre, p4,  p5)
-1.111881270890998e-16

ShiftedJacobi

Shifted Jacobi polynomial constructor. P̃(x) = P(2x-1). Shifted Jacobi polynomials are orthogonal on $[0,1]$.

ShiftedLegendre

Shifted Legendre polynomial constructor. Shifted Legendre polynoials are orthogonal on $[0,1]$.

Laguerre{α, T <: Number}

The Laguerre polynomials have weight function x^α * exp(-x) over the domain [0, oo). The parameter α is specified through the constructor.

julia> using Polynomials, SpecialPolynomials

julia> p = Laguerre{1/2}([1,2,3])
Laguerre{0.5}(1⋅Lᵅ₀(x) + 2⋅Lᵅ₁(x) + 3⋅Lᵅ₂(x))

julia> convert(Polynomial, p)
Polynomial(9.625 - 9.5*x + 1.5*x^2)

The Laguerre polynomials are the case α=0.

julia> using Polynomials, SpecialPolynomials

julia> p = Laguerre{0}([1,2,3])
Laguerre{0}(1⋅L₀(x) + 2⋅L₁(x) + 3⋅L₂(x))

julia> convert(Polynomial, p)
Polynomial(6.0 - 8.0*x + 1.5*x^2)

julia> phi(u, i) = derivative(u) -  u # verify Rodrigues' formula for small n; n! L_n = (d/dx-1)^n x^n
phi (generic function with 1 method)

julia> x = Polynomial(:x)
Polynomial(1.0*x)

julia> n = 7
7

julia> factorial(n) * basis(Laguerre{0}, n) - foldl(phi, 1:n, init=x^n)
LaurentPolynomial(-5.4569682106375694e-12 + 1.4551915228366852e-11*x - 7.275957614183426e-12*x²)

Hermite

The Hermite polynomials have two versions the physicists (Hermite or H) and the probablalists (ChebyshevHermite or Hₑ). They are related through Hᵢ(x) = 2^(i/2) Hₑᵢ(√2 x).

The Hermite polynomials have weight function w(x)=exp(-x^2/2) and domain the real line.

Examples

julia> using Polynomials,  SpecialPolynomials

julia> x = variable(Polynomial{Rational{Int}})
Polynomial(x)

julia> [basis(Hermite, i)(x) for i in 0:5]
6-element Vector{Polynomial{Float64, :x}}:
 Polynomial(1.0)
 Polynomial(2.0*x)
 Polynomial(-2.0 + 4.0*x^2)
 Polynomial(-12.0*x + 8.0*x^3)
 Polynomial(12.0 - 48.0*x^2 + 16.0*x^4)
 Polynomial(120.0*x - 160.0*x^3 + 32.0*x^5)

julia> [basis(ChebyshevHermite, i)(x) for i in 0:5]
6-element Vector{Polynomial{Float64, :x}}:
 Polynomial(1.0)
 Polynomial(1.0*x)
 Polynomial(-1.0 + 1.0*x^2)
 Polynomial(-3.0*x + 1.0*x^3)
 Polynomial(3.0 - 6.0*x^2 + 1.0*x^4)
 Polynomial(15.0*x - 10.0*x^3 + 1.0*x^5)

ChebyshevHermite

Type for the Probabalist's Hermite polynomials.

Bessel{α}

Implements the Bessel polynomials, introduced by Krall and Frink (with b=2). The case a=2 corresponds to the Bessel polynomials of Wikipedia. The Bessel polynomials are not orthogonal over a domain of the real line, rather over an arbitrary curve in the complex plane enclosing the origin. The weight function is ρ(x)=(2πi)^(-1)∑Γ(α)/Γ(α+n-1)(-β/x)^n, where β=2.

julia> using Polynomials, SpecialPolynomials

julia> 𝐐 = Rational{Int}
Rational{Int64}

julia> x = variable(Polynomial{𝐐})
Polynomial(x)

julia> [basis(Bessel{3//2, 𝐐}, i)(x) for i in 0:5]
6-element Vector{Polynomial{Rational{Int64}, :x}}:
 Polynomial(1//1)
 Polynomial(1//1 + 3//4*x)
 Polynomial(1//1 + 5//2*x + 35//16*x^2)
 Polynomial(1//1 + 21//4*x + 189//16*x^2 + 693//64*x^3)
 Polynomial(1//1 + 9//1*x + 297//8*x^2 + 1287//16*x^3 + 19305//256*x^4)
 Polynomial(1//1 + 55//4*x + 715//8*x^2 + 10725//32*x^3 + 182325//256*x^4 + 692835//1024*x^5)

Charlier{μ}

References: Koekoek and Swarttouw §1.12

 Krawchouk{p,𝐍}

Also spelled Krawtchouk, Kravhcuk,….

References: Koekoek and Swarttouw §1.10; see also Coleman for a different parameterization.

Meixner{γ,μ}

References: Koekoek and Swarttouw §1.9

Hahn{α,β,𝐍}

References: Koekoek and Swarttouw §1.5

DiscreteChebyshev

This uses p22 of Koepf and Schmersau to define a two-parameter family of non orthogonal polynomials.

julia> using Polynomials, SpecialPolynomials

julia> import SpecialPolynomials: Δₓ, ∇ₓ

julia> α,β = 1/2, 1
(0.5, 1)

julia> P  = DiscreteChebyshev{α,β}
Polynomials.MutableDensePolynomial{SpecialPolynomials.DiscreteChebyshevBasis{0.5, 1}}

julia> i = 5
5

julia> yᵢ = basis(P, i)
DiscreteChebyshev{0.5,1}(1.0⋅K⁽ᵅᵝ⁾₅(x))

julia> x = variable(P)
DiscreteChebyshev{0.5,1}(- 2.0⋅K⁽ᵅᵝ⁾₀(x) + 2.0⋅K⁽ᵅᵝ⁾₁(x))

julia> a,b,c,d,e = SpecialPolynomials.abcde(P)
(a = 0, b = 0, c = 1, d = 0.5, e = 1)

julia> λᵢ  = -(a*i*(i-1)  + d*i)
-2.5

julia> Δₓ(∇ₓ(yᵢ)) +  (α*x + β) * Δₓ(yᵢ) ≈ -λᵢ*yᵢ # p22: "are not orthogonal, but satisfy the difference equation..."
true

FallingFactorial{T}

Construct a polynomial with respect to the basis x⁰̲, x¹̲, x²̲, … where xⁱ̲ = x ⋅ (x-1) ⋅ (x-2) ⋯ (x-i+1) is the falling Pochhammer symbol. See Falling factorial for several facts about this polynomial basis.

In Koepf and Schmersau connection coefficients between the falling factorial polynomial system and classical discrete orthogonal polynomials are given.

Examples

julia> using Polynomials, SpecialPolynomials

julia> p = basis(FallingFactorial, 3)
Polynomials.MutableDensePolynomial(1.0⋅x³̲)

julia> x = variable(Polynomial)
Polynomial(1.0*x)

julia> p(x) ≈ x*(x-1)*(x-2)
true

weight_function(p)
weight_function(::Type{P})

For an orthogonal polynomial type, a function w with ∫ B_n(t) B_m(t) w(t) dt = 0 when n and m are not equal.

generating_function(p)
generating_function(::Type{P})

The generating function is a function defined by: (t,x) -> sum(t^n Pn(x) for n in 0:oo).

abcde

A named tuple returning the constants a,b,c,d,e for a CCOP type with (a⋅x²+b⋅x+c)*P₍ᵢ₊₂₎'' + (d⋅x + e) * P₍ᵢ₊₁₎ + λᵢ Pᵢ = 0.

An(::Type{P},n)

Orthogonal polynomials defined by a weight function satisfy a three point recursion formula of the form:

P_{n+1} = (A_n x + B_n) P_{n} - C_n P_{n-1}

If the polynomials are monic, this is usually parameterized as:

π_{n+1} = (x - α̃_n) π_n - β̃_n π_{n-1}

These functions are used through recursion when evaluating the polynomials, converting to Polynomial format, for constructing the Vandermonde matrix, for construction the Jacobi matrix, and elsewhere.

Bn(::Type{B},n)
Bn(p::P, n)

cf. An()

Cn(::Type{B},n)
Cn(p::P, n)

cf. An()

jacobi_matrix(::Type{P}, n)

The Jacobi Matrix is a symmetric tri-diagonal matrix. The diagonal entries are the alphaᵢ values, the off diagonal entries, the square root of the betaᵢ values. This matrix has the properties that

• the eigenvalues are the roots of the corresponding basis vector. As these roots are important in quadrature, and finding eigenvalues of a symmetric tri-diagonal matrix yields less error than finding the eigenvalues of the companion matrix, this can be used for higher degree basis polynomials.
• the normalized eigenvectors have initial term proportional to the weights in a quadrature formula

gauss_nodes_weights(::Type{P}, n)

Returns a tuple of nodes and weights for Gauss quadrature for the given orthogonal type.

When loaded, the values are computed through the FastGaussQuadrature package.

For some types, a method from A. Glaser, X. Liu, and V. Rokhlin. "A fast algorithm for the calculation of the roots of special functions." SIAM J. Sci. Comput., 29 (2007), 1420-1438. is used.

For others the Jacobi matrix, Jn, for which the Golub-Welsch] algorithm The nodes are computed from the eigenvalues of Jn, the weights a scaling of the first component of the normalized eigen vectors (β_0 * [v[1] for v in vs])

lagrange_barycentric_nodes_weights(::Type{<:SpecialPolynomial}, n::Int)

Return a collection of n+1 nodes and simplified weights the given family of polynomials. There are explicit formula for Chebyshev and Chebyshev, for other classical continuous orthogonal polynomials there is a relationship between the Gauss nodes and weights and the barycentric ones that is used.

AbstractInterpolatingPolynomial{T,X}

Abstract type for interpolating polynomials.

These are polynomial representations of p(x) satisfying p(x_i) = y_i for a specified set of x values and y values.

For a collection of points (x_0,y_0), ..., (x_n, y_n) there is a unique polynomial of degree n or less satisfying p(x_i)=y_i. This fact allows the specification of p(x) using a vector of coefficients relative to some set of basis vectors.

The two main types, Lagrange and Newton, store the nodes within the instance. In particular, the type does not contain all the information needed to describe the family. So methods like convert(::Type, p) will not work. Use fit(Type, xs, p), as appropriate, instead.

Lagrange(xs, [ws], coeffs, [var])

Lagrange interpolation of points (xᵢ, fᵢ) for i ∈ 0..n.

• xs, coeffs: the interpolating coordinates.
• ws: weights used in the barycentric representation. (From SpecialPolynomials.lagrange_barycentric_weights or SpecialPolynomials.lagrange_barycentric_nodes_weights.)
• var: the polynomial indeterminate

Extended help

The Lagrange interpolation of points (xᵢ, fᵢ) for i ∈ 0:n is the polynomial p(x) = ∑ᵢ lⱼ(x) fⱼ.

The basis vectors lⱼ(x) are 1 on xⱼ and 0 on xᵢ when i ≠ j. That is, lⱼ(x) = Π_{i ≠ j}(x-xᵢ)/Π_{i ≠j}(xⱼ-xᵢ). These can be rewritten in terms of weights, wⱼ, depending on the xᵢ only, yielding lⱼ = l(x) wⱼ/(x - xⱼ) with l(x) = Π(x-xᵢ). Going further, yields the barycentric formula:

p(x) = (∑ wⱼ / (x - xⱼ) ⋅ fⱼ) / (∑ wⱼ / (x - xⱼ) ).

This representation has several properties, as detailed in Berrut and Trefethen Barycentric Lagrange Interpolation.

Examples

julia> using Polynomials, SpecialPolynomials

julia> p =  Lagrange([1,2,3], [1,2,3])
Lagrange(1⋅ℓ_0(x) + 2⋅ℓ_1(x) + 3⋅ℓ_2(x))

julia> p.([1,2,3]) # the coefficients
3-element Vector{Int64}:
 1
 2
 3

julia> convert(Polynomial,  p)
Polynomial(1.0*x)

The instances hold the nodes and weights, which are necessary for representation, so the type alone can not be used for functions such as variable or convert(Lagrange, ...). For the former we can use an instance, for the latter we can use fit:

julia> p =  Lagrange([1,2,3], [1,2,3])
Lagrange(1⋅ℓ_0(x) + 2⋅ℓ_1(x) + 3⋅ℓ_2(x))

julia> variable(p)
Lagrange(1⋅ℓ_0(x) + 2⋅ℓ_1(x) + 3⋅ℓ_2(x))

julia> q = Polynomial([0,0,1])
Polynomial(x^2)

julia> qq = fit(Lagrange, p.xs, p.ws, q)
Lagrange(1⋅ℓ_0(x) + 4⋅ℓ_1(x) + 9⋅ℓ_2(x))

julia> convert(Polynomial, qq)
Polynomial(1.0*x^2)

For a given set of nodes, SpecialPolynomials.lagrange_barycentric_weights can compute the weights. For all but modest values of n, interpolating polynomials suffer from the Runge phenomenon unless the nodes are well chosen. (They should have asymptotic density of 1/√(1-x^2) over [-1,1].) For P=Chebyshvev and P=ChebyshevU, the function SpecialPolynomials.lagrange_barycentric_nodes_weights(P, n) will return a good choice of n+1 points over [-1,1] along with precomputed weights.

julia> xs, ws = SpecialPolynomials.lagrange_barycentric_nodes_weights(SpecialPolynomials.ChebyshevBasis, 64);


julia> f(x) = exp(-x)*sinpi(x)
f (generic function with 1 method)

julia> p = fit(Lagrange, xs, f.(xs));


julia> degree(p)
64

julia> maximum(abs.(f(x) - p(x) for x in range(-1, 1, length=20))) <= 1e-14
true

Newton{N,S,T,X}

A Newton interpolating polynomial uses a basis 1, (x-x_0), (x-x_0)(x-x_1), ..., (x-x0)(x-x1)⋅⋅⋅(x-x_{n-1}) and coefficients (in forward form) f[x_0], f[x_0,x_1], ...,f[x_0,...,x_n]. The Newton class stores the nodes (after sorting) and the Newton tableau used to generate the coefficients on fitting.

The easiest way to construct an instance is with fit, as in:

julia> using Polynomials, SpecialPolynomials

julia> xs = [1,2,3,4]; f(x)= x^3 - 2x + 1;

julia> p = fit(Newton, xs, f)
Newton(5.0⋅p_1(x) + 6.0⋅p_2(x) + 1.0⋅p_3(x))

julia> p.(xs) == f.(xs)  # p interpolates
true

julia> convert(Polynomial, p)
Polynomial(1.0 - 2.0*x + 1.0*x^3)

Bernstein{N, T, X}

A Bernstein polynomial is a polynomial expressed in terms of Bernstein basic polynomials. For each degree, 𝐍, this is a set of 𝐍+1 degree 𝐍 polynomials of the form: β_{𝐍,ν} = (ν choose 𝐍) x^ν (1-x)^{𝐍-ν}, 0 ≤ x ≤ 1.

The Bernstein{𝐍,T} type represents a polynomial of degree 𝐍 or less with a linear combination of the basis vectors using coefficients of type T.

julia> using Polynomials, SpecialPolynomials

julia> p = basis(Bernstein{3},  2)
Bernstein(1⋅β₃,₂(x))

julia> convert(Polynomial, p)
Polynomial(3.0*x^2 - 3.0*x^3)

DualBernstein{N, α, β, T, X}

Given the inner product <f,g> = ∫₀¹ (1-x)ᵅ xᵝ f(x) g(x) dx and the Bernstein polynomials Bᵢⁿ(x) the dual Bernstein polynomials satisfy <Bᵢⁿ, Dⱼⁿ(x;α,β)> = δᵢⱼ

Example

julia> n,α,β = 5, 0, 0;

julia> D = DualBernstein{n,α,β}
DualBernstein{5, 0, 0}

julia> bi = basis(D, 3)
DualBernstein(1.0⋅βᵅᵝ₅,₃(x))

julia> bi(0.2)
7.910400000000036

The Bernstein-Bezier form of f minimizes the value of the least-square error for the α-β norm.

julia> f(x) = sinpi(x);

julia> n, α, β = 5, 1/2, 1/2
(5, 0.5, 0.5)

julia> B, D = Bernstein{n}, DualBernstein{n,α,β};

julia> Iₖ = [ip(f, basis(D,k), α, β) for k in 0:n];

julia> pₙ = B(Iₖ);

julia> λ = x -> f(x) - pₙ(x);

julia> SpecialPolynomials.innerproduct(ShiftedJacobi{α, β}, λ, λ)
0.00017514530881540565

Reference

The implementation follows that of Chudy and Woźny.