### Discrete Chebyshev Polynomial Usage Example

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/discrete-orthogonal.md

Shows how to instantiate and evaluate Discrete Chebyshev polynomials.

```julia
dc = DiscreteChebyshev{10}([1, 2, 3])

# Evaluate
dc(0)
dc(5)
dc(10)

# Basis elements
basis(DiscreteChebyshev{10}, n)
```

--------------------------------

### Meixner Polynomial Usage Example

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/discrete-orthogonal.md

Demonstrates creating and evaluating Meixner polynomials with specified parameters.

```julia
# Meixner with β=2, c=0.4
m = MeixnerPolynomial{2, 0.4}([1, 2, 3])

# Evaluate
m(0)
m(5)
m(10)

# Basis element
basis(MeixnerPolynomial{2, 0.4}, 3)
```

--------------------------------

### Parametric Type Delegation Example

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/types.md

Demonstrates how shifted and variant polynomial types are related through the 'ϟ' delegation macro, using LegendreBasis as an example. Shows that variants like MonicLegendreBasis and OrthonormalLegendreBasis delegate to the base LegendreBasis.

```julia
LegendreBasis
├── MonicLegendreBasis (ϟ(MonicLegendreBasis) = LegendreBasis)
├── OrthonormalLegendreBasis (ϟ(OrthonormalLegendreBasis) = LegendreBasis)
└── ShiftedLegendreBasis (ϟ(ShiftedLegendreBasis) = LegendreBasis)
```

--------------------------------

### Example: Interpolate Points and Evaluate

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/interpolating-polynomials.md

Demonstrates fitting a Newton polynomial to discrete points and then evaluating it at a new point.

```julia
# Interpolate points
xs = [1.0, 2.0, 3.0, 4.0]
ys = [1.0, 4.0, 9.0, 16.0]
p = fit(Newton, xs, ys)

# Evaluate
p(2.5)
```

--------------------------------

### Install SpecialPolynomials Package

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Install the SpecialPolynomials package using Julia's package manager. Requires Julia 1.10 or later.

```julia
using Pkg
Pkg.add("SpecialPolynomials")
```

--------------------------------

### Example: Interpolate Function and Evaluate

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/interpolating-polynomials.md

Demonstrates fitting a Newton polynomial to a function sampled over a range of nodes and evaluating it.

```julia
# Interpolate function
f(x) = sin(x)
xs = range(0, π, length=10)
p = fit(Newton, xs, f)
p(π/2)
```

--------------------------------

### Bezier Curve Example

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/docs/src/index.md

Example demonstrating the use of Bernstein polynomials to define a Bezier curve.

```APIDOC
## Example of a Bezier curve

Parameterized by `r(t) = ∑₀ᴺ bᵢBᵢ(t)`:

```julia
using Plots, Polynomials, SpecialPolynomials

bs = [[220, 260], [220, 40], [35, 100],  [120, 140]]

p = Bernstein(bs)
ts = range(0, stop=1, length=50)
ps = p.(ts)
xs, ys = [[pᵢ[1] for pᵢ ∈ ps], [pᵢ[2] for pᵢ ∈ ps]]
p = plot(xs, ys, legend=false)
scatter!(p, [b[1] for b in bs], [b[2] for b in bs])
```
```

--------------------------------

### DualBernstein Polynomial Example and Conversion

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/bernstein-polynomials.md

Demonstrates constructing a DualBernstein polynomial, evaluating it, and converting it to a standard Bernstein polynomial.

```julia
p = DualBernstein{3}([1, 2, 3, 4])
p(0.5)  # evaluate

# Convert between Bernstein and dual Bernstein
b = Bernstein{3}([1, 2, 3, 4])
db = convert(DualBernstein{3}, b)
```

--------------------------------

### Plot Legendre Polynomial Basis Functions

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/docs/src/examples.md

Generates plots for the first few Legendre polynomial basis functions using Plots.jl. This example requires Plots.jl to be installed and loaded.

```example
using Plots, Polynomials, SpecialPolynomials  # hide
n = 6
ψ₁, ψ₂, ψ₃, ψ₄, ψ₅, ψ₆ = basis.(Legendre, 0:n - 1)
kw = (xlabel="r", ylabel="ψₙ(r)", legend=false)
sps = [
  plot( ψ₁, label="ψ₁(r)"; kw...),
  plot( ψ₂, label="ψ₂(r)"; kw...),
  plot( ψ₃, label="ψ₃(r)"; kw...),
  plot( ψ₄, label="ψ₄(r)"; kw...),
  plot( ψ₅, label="ψ₅(r)"; kw...),
  plot( ψ₆, label="ψ₆(r)"; kw...),
]
p = plot(sps..., layout=(3, 2))
savefig(p, "LegendrePolynomials.svg"); nothing  # hide
```

--------------------------------

### Bezier Curve Example with Bernstein Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/docs/src/index.md

Demonstrates how to define and plot a Bezier curve using Bernstein polynomials. Requires Plots, Polynomials, and SpecialPolynomials packages.

```julia
using Plots, Polynomials, SpecialPolynomials;  # hide
bs = [[220, 260], [220, 40], [35, 100],  [120, 140]]

p = Bernstein(bs)
ts = range(0, stop=1, length=50)
ps = p.(ts)
xs, ys = [[pᵢ[1] for pᵢ ∈ ps], [pᵢ[2] for pᵢ ∈ ps]]
p = plot(xs, ys, legend=false)
scatter!(p, [b[1] for b in bs], [b[2] for b in bs])
savefig("bezier.svg"); nothing  # hide
```

--------------------------------

### Bernstein Constructor Example

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/bernstein-polynomials.md

Demonstrates how to construct Bernstein polynomials using both type-specified and type-unspecified constructors, including specifying a custom variable.

```julia
using Polynomials, SpecialPolynomials

# Degree 3 Bernstein polynomial
p = Bernstein{3}([1.0, 2.0, 3.0, 4.0])

# With custom variable
q = Bernstein{2}([1, 2, 3], :s)
```

```julia
# Automatically determines degree 2 (from 3 coefficients)
p = Bernstein([1.0, 2.0, 3.0])

# Equivalent to
p2 = Bernstein{2}([1.0, 2.0, 3.0])
```

--------------------------------

### Constructing Jacobi Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/core-types.md

Provides an example of constructing a Jacobi polynomial with specific parameters and coefficients. ChebyshevU is noted as a special case.

```julia
# Jacobi with parameters
j = Jacobi{1/2, 1/2}([1, 2, 3])  # ChebyshevU is a special case
```

--------------------------------

### Falling Factorial Polynomial Usage Example

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/discrete-orthogonal.md

Illustrates the creation and evaluation of polynomials using the falling factorial basis.

```julia
ff = FallingFactorial([1, 2, 3, 4])

# Evaluate at integers
ff(0)   # = 1
ff(1)   # = 1 + 2⋅1
ff(2)   # = 1 + 2⋅1 + 3⋅2⋅1
ff(3)   # = 1 + 2⋅1 + 3⋅2⋅1 + 4⋅3⋅2⋅1
```

--------------------------------

### Compute Gamma Function Values

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/utility-functions.md

Examples of computing the gamma function for integer and fractional inputs, demonstrating its relationship with factorials and constants like sqrt(pi).

```julia
gamma(5)     # 24 = 4!
gamma(1/2)   # √π ≈ 1.7725
gamma(3)     # 2! = 2
```

--------------------------------

### Quadratic Bernstein Polynomial Example

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/bernstein-polynomials.md

Shows the construction and evaluation of a quadratic (degree 2) Bernstein polynomial, representing a parabolic arc.

```julia
p = Bernstein{2}([p0, p1, p2])
# Represents: p0(1-x)² + 2p1(1-x)x + p2 x²

# Parabolic arc with control points p0, p1, p2
p = Bernstein{2}([0.0, 0.5, 1.0])
p(0.0)    # 0.0
p(0.5)    # 0.5
p(1.0)    # 1.0
```

--------------------------------

### Cubic Bernstein Polynomial Example

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/bernstein-polynomials.md

Demonstrates the construction of a cubic (degree 3) Bernstein polynomial, representing a cubic Bézier curve with 4 control points.

```julia
p = Bernstein{3}([0.0, 0.33, 0.67, 1.0])
# Cubic Bézier curve with 4 control points
```

--------------------------------

### Visualize Interpolation with Different Bases

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/docs/src/examples.md

Plots an example function `f(x)` and its interpolations using Chebyshev, ChebyshevU, Legendre, and Newton bases with degree 3. This visualization helps to observe the Runge phenomenon and the effectiveness of different node choices.

```julia
using Plots, Polynomials, SpecialPolynomials;  # hide
f(x) = exp(-x)*sinpi(x)
plot(f, -1, 1, legend=false, color=:black, linewidth=3)
p=fit(Val(:interpolating), Chebyshev, f, 3); plot!(p, color=:blue)
p=fit(Val(:interpolating), ChebyshevU, f, 3); plot!(p, color=:red)
fit(Val(:interpolating), Legendre, f, 3); plot!(p, color=:green)
xs = [-0.5, 0.0, 0.5]
p=fit(Newton, xs, f);
ts = range(-1, 1, length=100); plot!(ts, p.(ts), color=:brown)
savefig("fitting.svg"); nothing  # hide
```

--------------------------------

### Linear Bernstein Polynomial Example

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/bernstein-polynomials.md

Illustrates the construction and evaluation of a linear (degree 1) Bernstein polynomial, representing linear interpolation.

```julia
p = Bernstein{1}([a, b])
# Represents: a(1-x) + bx = a + (b-a)x
# Linear interpolation between a and b on [0,1]

p = Bernstein{1}([0.0, 1.0])
p(0.5)  # 0.5 (linear interpolation)
```

--------------------------------

### Get Individual Basis Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Create specific basis polynomials (e.g., Legendre P₀, P₁, P₂) using the `basis` function. These can then be used to construct linear combinations manually.

```julia
# Legendre basis elements
P0 = basis(Legendre, 0)
P1 = basis(Legendre, 1)
P2 = basis(Legendre, 2)

# Evaluate
P0(0.5)  # 1.0
P1(0.5)  # 0.5
P2(0.5)  # -0.125

# Construct linear combination manually
p = 1.5*P0 + 2.0*P1 + 0.5*P2
p(0.5)
```

--------------------------------

### Least Squares Fit for Chebyshev Polynomials (Wavy Function)

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/docs/src/examples.md

Demonstrates a least squares fit of a Chebyshev polynomial to a 'wavy' function. The example shows that a lower degree fit (50) can be inaccurate, indicated by the assertion failing.

```jldoctest
julia> f(x) = sin(6x) + sin(60*exp(x))
f (generic function with 1 method)

julia> p50 = fit(Val(:lsq), Chebyshev{Float64}, f, 50);

julia> maximum(norm(p50(x)-f(x) for x in range(-1,1,length=500))) <= sqrt(eps()) # cf. graph below
false
```

--------------------------------

### Bernstein Polynomial Arithmetic

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/bernstein-polynomials.md

Shows examples of basic arithmetic operations on Bernstein polynomials, including addition, subtraction, multiplication, and scalar multiplication.

```julia
p = Bernstein{2}([1, 2, 3])
q = Bernstein{2}([1, 1, 1])

p + q    # addition (promotes to same degree)
p - q    # subtraction
p * q    # multiplication (result has degree 2*(N₁+N₂))
2 * p    # scalar multiplication
```

--------------------------------

### Retrieving Polynomial Domain Information

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Illustrates how to get the domain for Legendre, Chebyshev, and Bernstein polynomial types using the `domain` function.

```julia
using Polynomials: domain

p = Legendre([1, 2, 3])
domain(Legendre)    # Returns interval [-1, 1]

p = Chebyshev([1, 2, 3])
domain(Chebyshev)   # Returns interval (-1, 1) [open]

p = Bernstein([1, 2, 3])
domain(Bernstein)   # Returns interval [0, 1]
```

--------------------------------

### Chebyshev Polynomial Initialization

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/README.md

Demonstrates how to initialize a Chebyshev polynomial using its coefficients. Coefficients are stored lowest-degree-first.

```julia
p = Chebyshev([1, 2, 3])  # 1⋅T₀(x) + 2⋅T₁(x) + 3⋅T₂(x)
```

--------------------------------

### Constructing Chebyshev Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/core-types.md

Demonstrates how to construct a Chebyshev polynomial from a list of coefficients and a variable symbol. It also shows how to evaluate the polynomial.

```julia
using Polynomials, SpecialPolynomials

# Chebyshev polynomial: 1⋅T₀(x) + 2⋅T₁(x) + 3⋅T₂(x)
p = Chebyshev([1, 2, 3], :x)

# Evaluate
p(0.5)  # type-specific evaluation via Clenshaw algorithm
```

--------------------------------

### Get Basis Element

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Retrieve a specific basis element for a given polynomial type and degree.

```julia
basis(Legendre, n)
```

--------------------------------

### Get Orthonormal Legendre Basis

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Retrieve the orthonormal variant of the Legendre polynomial basis for a given degree.

```julia
basis(OrthonormalLegendre, n)
```

--------------------------------

### Get Monic Legendre Basis

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Retrieve the monic variant of the Legendre polynomial basis for a given degree.

```julia
basis(MonicLegendre, n)
```

--------------------------------

### Instantiate and Use Laguerre Polynomial Basis Elements

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/orthogonal-polynomials.md

Demonstrates how to instantiate a Laguerre polynomial with given coefficients and how to access its basis elements for different degrees. Use this to create and evaluate Laguerre polynomials.

```julia
l = Laguerre([1, 2, 3])

# Basis elements
basis(Laguerre, 0)(x)  # 1
basis(Laguerre, 1)(x)  # 1 - x
basis(Laguerre, 2)(x)  # (x² - 4x + 2)/2
```

--------------------------------

### Computing Hypergeometric Functions

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Demonstrates the computation of hypergeometric functions using `pFq` and their application in generating Legendre polynomials.

```julia
import SpecialPolynomials: pFq, Pochhammer

# Compute hypergeometric function
pFq((1/3, 2/3), (5/6,), 27/32)  # ≈ 8/5

# Use with polynomials
x = variable(Polynomial)
n = 5

# Legendre via hypergeometric
result = pFq((-n, n+1), (1,), (1-x)/2)

# Should match Legendre basis
expected = basis(Legendre, n)

# Verify at test point
result(0.3) ≈ expected(0.3)  # true
```

--------------------------------

### Convert Polynomial Type

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Convert a polynomial from one type to another, for example, from Legendre to a standard Polynomial type.

```julia
convert(Polynomial, p)
```

--------------------------------

### Lagrange Polynomial Construction via fit()

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/interpolating-polynomials.md

Constructs a Lagrange polynomial using the `fit` function, either by directly providing points or by fitting a function to sampled points.

```julia
fit(::Type{Lagrange}, xs::Vector, ys::Vector)
```

```julia
fit(::Type{Lagrange}, xs::Vector, ws::Vector, f)
```

```julia
# Direct interpolation of points
xs = [0.0, 0.5, 1.0]
ys = [0.0, 0.25, 1.0]  # y = x²
p = fit(Lagrange, xs, ys)

# Interpolation of a function
f(x) = sin(π*x)
xs = range(-1, 1, length=11)
p = fit(Lagrange, xs, f.(xs))
```

--------------------------------

### Getting the nth Basis Element

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/core-types.md

Shows how to retrieve a specific basis element for a given polynomial type and index.

```julia
# Get the nth basis element
b₅ = basis(Legendre, 5)
```

--------------------------------

### Import Necessary Packages

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Import SpecialPolynomials for polynomial operations and Polynomials for the Polynomial type. Plots is optional for visualization.

```julia
using SpecialPolynomials
using Polynomials  # Also recommended for Polynomial type

# Optional: for visualization
using Plots
```

--------------------------------

### gamma(z)

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/utility-functions.md

Computes the gamma function, an alias for Γ. It is imported from SpecialFunctions.jl and re-exported. Commonly used values and examples are provided.

```APIDOC
## gamma(z)

### Description
Computes the gamma function, an alias for `Γ`.

### Signature
`gamma(z) -> Number`

### Details
Imported from `SpecialFunctions.jl` and re-exported.

Computes the gamma function defined by the integral:
```
Γ(z) = ∫₀^∞ tᶻ⁻¹ e⁻ᵗ dt
```

**Commonly used values:**
- `gamma(n)` = `(n-1)!` for positive integer `n`
- `gamma(1/2)` = `√π`
- `gamma(n+1)` = `n * gamma(n)`

### Example
```julia
gamma(5)     # 24 = 4!
gamma(1/2)   # √π ≈ 1.7725
gamma(3)     # 2! = 2
```
```

--------------------------------

### Calculate Gauss-Kronrod Nodes and Weights for Legendre Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/docs/src/examples.md

Demonstrates how to obtain the nodes and weights for Gauss-Kronrod quadrature using the `gauss_nodes_weights` function for Legendre polynomials. This is useful for numerical integration.

```julia
julia> xs, ws = SpecialPolynomials.gauss_nodes_weights(Legendre, 4)
([-0.86113631159405, -0.33998104358485537, 0.33998104358485715, 0.8611363115940526], [0.34785484513745607, 0.6521451548625448, 0.652145154862546, 0.3478548451374536])

julia> basis(Legendre, 4).(xs)
4-element Vector{Float64}:
 -1.199040866595169e-14
  1.6930901125533637e-15
 -1.6653345369377348e-15
  1.1102230246251565e-16

julia> f(x) = x^7 - x^6; F(x) = x^8/8 - x^7/7;

julia> sum(f(x)*w for (x,w) in zip(xs, ws)) - (F(1) - F(-1))
3.3306690738754696e-15
```

--------------------------------

### Legendre Polynomial Initialization and Evaluation

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/orthogonal-polynomials.md

Initializes a Legendre polynomial from a list of coefficients and evaluates it at a given point. The coefficients represent the polynomial in terms of Legendre basis polynomials.

```julia
p = Legendre([1, 2, 3])  # 1⋅P₀ + 2⋅P₁ + 3⋅P₂
p(0.0)
```

--------------------------------

### Construction from Basis Element

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/bernstein-polynomials.md

Shows how to construct a single Bernstein basis function using the `basis` function.

```julia
p = basis(Bernstein{3}, 1)  # the β₃,₁(x) basis function
```

--------------------------------

### Approximate sin(πx) using Bernstein Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/bernstein-polynomials.md

Demonstrates how to approximate the function sin(πx) on the interval [0,1] using Bernstein polynomials. This involves creating a Bernstein polynomial from sample points and plotting the approximation against the original function.

```julia
using SpecialPolynomials, Polynomials, Plots

# Function to approximate
f(x) = sin(π*x)

# Sample points for interpolation
n = 5
xs = range(0, 1, length=n+1)
ys = f.(xs)

# Create Bernstein polynomial
p = Bernstein(ys)

# Evaluate on finer grid
test_xs = range(0, 1, length=200)
p_vals = p.(test_xs)
f_vals = f.(test_xs)

# Plot comparison
plot(test_xs, f_vals, label="f(x)", linewidth=2)
plot!(test_xs, p_vals, label="Bernstein approximation", linewidth=2)
scatter!(xs, ys, label="interpolation points", markersize=6)
```

--------------------------------

### Lagrange Interpolation

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/core-types.md

Demonstrates how to perform Lagrange interpolation by providing a set of x and y coordinates. The interpolated value at a given point can then be calculated.

```julia
# Lagrange interpolation
xs = [1.0, 2.0, 3.0]
ys = [1.0, 4.0, 9.0]
q = Lagrange(xs, ys)
q(1.5)  # interpolated value
```

--------------------------------

### Get Weight Function for Chebyshev Basis

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/README.md

Retrieves the weight function associated with the Chebyshev polynomial basis. This function defines the orthogonality relations.

```julia
weight_function(ChebyshevBasis)  # → x -> 1/√(1-x²)
```

--------------------------------

### DualBernstein Polynomial Constructor

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/core-types.md

Shows the constructor for DualBernstein polynomials, used in approximation theory.

```julia
DualBernstein{N}(coeffs::Vector, var::Symbol=:x)
```

--------------------------------

### Parameter Table Format

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/README.md

This is the standardized format for all function/method documentation, detailing parameter names, types, requirements, defaults, and descriptions.

```markdown
| Parameter | Type | Required | Default | Description |
|-----------|------|----------|---------|-------------|
```

--------------------------------

### Legendre Basis Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/orthogonal-polynomials.md

Demonstrates how to access and evaluate Legendre basis polynomials for different degrees. These are fundamental building blocks for Legendre polynomials.

```julia
basis(Legendre, 0)(x)  # 1
basis(Legendre, 1)(x)  # x
basis(Legendre, 2)(x)  # (3x² - 1)/2
basis(Legendre, 3)(x)  # (5x³ - 3x)/2
```

--------------------------------

### Get Chebyshev Basis Vectors

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/orthogonal-polynomials.md

Retrieves the first few basis vectors for Chebyshev polynomials. These represent the polynomials T₀(x), T₁(x), and T₂(x).

```julia
# Get the nth basis vector
basis(Chebyshev, 0)(x)  # returns 1
basis(Chebyshev, 1)(x)  # returns x
basis(Chebyshev, 2)(x)  # returns 2x² - 1
```

--------------------------------

### Constructing DualBernstein Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/bernstein-polynomials.md

Shows how to construct DualBernstein polynomials with a specified degree N or inferred degree.

```julia
DualBernstein{N}(coeffs::Vector, var::Symbol=:x)
```

```julia
DualBernstein(coeffs::Vector, var::Symbol=:x)  # infers N
```

--------------------------------

### Construct Basis Element

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/core-types.md

Constructs the n-th basis element of a given polynomial family type P. For example, `basis(Legendre, 3)` returns the third Legendre polynomial basis element.

```julia
basis(Legendre, 3)  # returns Legendre(1⋅P₃(x))

```

```julia
basis(Chebyshev, 5) # returns Chebyshev(1⋅T₅(x))

```

--------------------------------

### Check Hahn Polynomial Orthogonality

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/discrete-orthogonal.md

Verify the orthogonality of Hahn polynomials. This example computes the orthogonality integral for two different Hahn polynomials, demonstrating that the integral is approximately zero when the polynomial indices are different.

```julia
# Verify orthogonality of Hahn polynomials
α, β, N = 1, 1, 10
k_type = HahnPolynomial{α, β, N}

p = basis(k_type, 2)
q = basis(k_type, 5)

# Compute orthogonality integral
total = 0.0
for i in 0:N
    # Hahn weight (from hypergeometric distribution)
    w_i = 1.0  # simplified for example
    total += w_i * p(i) * q(i)
end

# total ≈ 0 for i ≠ j
```

--------------------------------

### Accessing Individual Bernstein Basis Functions

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/bernstein-polynomials.md

Illustrates how to construct and evaluate individual Bernstein basis functions (βₙ,ₖ(x)) for a given degree N.

```julia
# Construct individual basis functions
b = basis(Bernstein{3}, 0)  # β₃,₀(x) = (1-x)³
b(0.0)     # 1
b(1.0)     # 0

b = basis(Bernstein{3}, 3)  # β₃,₃(x) = x³
b(0.0)     # 0
b(1.0)     # 1
```

--------------------------------

### ∫(f, a, b)

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/utility-functions.md

A shorthand for numerical quadrature using QuadGK.jl, providing a convenient mathematical notation for integration.

```APIDOC
## `∫(f, a, b) -> Float64`

### Description
Shorthand for numerical quadrature using `quadgk`. This function wraps `QuadGK.quadgk` for convenient mathematical notation.

### Parameters
#### Path Parameters
- **f** (Function) - Required - Integrand
- **a** (Number) - Required - Lower limit of integration
- **b** (Number) - Required - Upper limit of integration

### Request Example
```julia
import SpecialPolynomials: ∫

# Integrate x² from 0 to 1
∫(x -> x^2, 0, 1)  # ≈ 0.333...

# Integrate with polynomial
p = Legendre([1, 2, 3])
∫(p, -1, 1)  # orthogonality integral
```

### Response
#### Success Response
- **Float64** - Integral value as `Float64`

**Alias:** The `_quadgk(f, a, b)` function provides the same functionality without special character notation.
```

--------------------------------

### Work with Different Polynomial Bases

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Create and evaluate polynomials from various bases like Chebyshev, Hermite, Jacobi, Legendre, and Bernstein. Note that Jacobi polynomials require parameters α and β.

```julia
# Chebyshev polynomial
c = Chebyshev([1, 2, 3])
c(0.5)

# Hermite polynomial
h = Hermite([1, 2, 3])
h(2.0)

# Jacobi polynomial with parameters α=1/2, β=1/2
j = Jacobi{1/2, 1/2}([1, 2, 3])
j(0.0)

# Legendre
l = Legendre([1, 2, 3])
l(-0.5)

# Bernstein (on [0,1])
b = Bernstein{3}([0, 0.25, 0.75, 1])
b(0.5)
```

--------------------------------

### Find Roots of Legendre Polynomial

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/docs/src/examples.md

Demonstrates how to find the roots of a Legendre polynomial and verify the residual. Use this to find the roots of a given polynomial.

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

julia> rts = roots(p); rts ≈ [-1, -1/5, 0]
true

julia> maximum(abs∘p, rts) <= 10eps() # small residual
true
```

--------------------------------

### Basic Lagrange Interpolation

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/interpolating-polynomials.md

Use `fit(Lagrange, xs, ys)` to create a Lagrange interpolating polynomial from a set of points. Evaluate the polynomial at intermediate points using standard function call syntax.

```julia
using SpecialPolynomials, Polynomials

# Define some points
xs = [0.0, 1.0, 2.0, 3.0]
ys = [1.0, 2.0, 5.0, 10.0]

# Create interpolant
p = fit(Lagrange, xs, ys)

# Evaluate at intermediate points
p(0.5)   # ≈ interpolated value
p(1.5)   # ≈ interpolated value
p(2.5)   # ≈ interpolated value
```

--------------------------------

### Jacobi Polynomial Constructor

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/core-types.md

Shows the constructor for Jacobi polynomials, which take parameters α and β that define the weight function.

```julia
Jacobi{α,β}(coeffs::Vector, var::Symbol=:x)
```

--------------------------------

### Gegenbauer Polynomial Constructor

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/core-types.md

Illustrates the constructor for Gegenbauer polynomials, which take a parameter λ that specifies the weight function.

```julia
Gegenbauer{λ}(coeffs::Vector, var::Symbol=:x)
```

--------------------------------

### Construct Legendre Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/docs/src/examples.md

Constructs the first few Legendre polynomials using their coefficients relative to the Legendre basis.

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> p0 = Legendre([1])
Legendre(1⋅P₀(x))

julia> p1 = Legendre([0,1])
Legendre(1⋅P₁(x))

julia> p2 = Legendre([0,0,1])
Legendre(1⋅P₂(x))

julia> p3 = Legendre([0,0,0,1])
Legendre(1⋅P₃(x))
```

--------------------------------

### Create a Legendre Polynomial

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Create a Legendre polynomial by providing its coefficients. The order of coefficients corresponds to increasing degree of Legendre basis polynomials.

```julia
p = Legendre([1, 2, 3])
```

--------------------------------

### Direct Lagrange Polynomial Construction with Weights

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/interpolating-polynomials.md

Constructs a Lagrange polynomial using precomputed barycentric weights. Ensure weights are computed using `lagrange_barycentric_weights` or `lagrange_barycentric_nodes_weights`.

```julia
Lagrange(xs::Vector, ws::Vector, coeffs::Vector, var::Symbol=:x)
```

```julia
using SpecialPolynomials, Polynomials

xs = [1.0, 2.0, 3.0]
ws = lagrange_barycentric_weights(xs)  # compute weights
ys = [1.0, 4.0, 9.0]  # function values

p = Lagrange(xs, ws, ys)
p(1.5)  # evaluate at intermediate point
```

--------------------------------

### Hermite Polynomial Initialization and Evaluation

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/orthogonal-polynomials.md

Initializes a Physicists' Hermite polynomial from coefficients and shows how to access its basis polynomials. Hermite polynomials are used with a Gaussian weight function.

```julia
# Physicists' Hermite
h = Hermite([1, 2, 3])  # 1⋅H₀ + 2⋅H₁ + 3⋅H₂
basis(Hermite, 0)(x)    # 1
basis(Hermite, 1)(x)    # 2x
basis(Hermite, 2)(x)    # 4x² - 2
basis(Hermite, 3)(x)    # 8x³ - 12x
```

--------------------------------

### Newton Polynomial Constructor

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/core-types.md

Illustrates the constructor for Newton interpolating polynomials, which use divided differences over specified nodes.

```julia
Newton(xs::Vector{S}, coeffs::Vector{T}, var::Symbol=:x)
```

--------------------------------

### Lagrange Polynomial Construction

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/interpolating-polynomials.md

Constructs a Lagrange polynomial using direct construction with weights, automatic weight computation, or the `fit()` function.

```APIDOC
## Lagrange Polynomial Construction

### Direct construction with weights

```julia
Lagrange(xs::Vector, ws::Vector, coeffs::Vector, var::Symbol=:x)
```

Constructs a Lagrange polynomial interpolating `(xs[i], coeffs[i])` using precomputed barycentric weights `ws`.

**Parameters:**

| Parameter | Type | Required | Description |
|-----------|------|----------|-------------|
| `xs` | `Vector{<:Number}` | Yes | Interpolation nodes |
| `ws` | `Vector{<:Number}` | Yes | Barycentric weights from `lagrange_barycentric_weights` or `lagrange_barycentric_nodes_weights` |
| `coeffs` | `Vector{<:Number}` | Yes | Values at nodes (lowest index first) |
| `var` | `Symbol` | No | Indeterminate variable, default `:x` |

**Returns:** `Lagrange` polynomial instance

**Example:**

```julia
using SpecialPolynomials, Polynomials

xs = [1.0, 2.0, 3.0]
ws = lagrange_barycentric_weights(xs)  # compute weights
ys = [1.0, 4.0, 9.0]  # function values

p = Lagrange(xs, ws, ys)
p(1.5)  # evaluate at intermediate point
```

### Construction with automatic weight computation

```julia
Lagrange(xs::Vector, coeffs::Vector, var::Symbol=:x)
```

Automatically computes barycentric weights from nodes.

**Example:**

```julia
xs = [1.0, 2.0, 3.0]
ys = [1.0, 4.0, 9.0]

p = Lagrange(xs, ys)  # weights computed automatically
p(1.5)
```

### Construction via `fit()`

```julia
fit(::Type{Lagrange}, xs::Vector, ys::Vector)
fit(::Type{Lagrange}, xs::Vector, ws::Vector, f)
```

**Parameters:**

| Parameter | Type | Required | Description |
|-----------|------|----------|-------------|
| `xs` | `Vector` | Yes | Interpolation nodes |
| `ys` | `Vector` | Yes | Function values at nodes OR |
| `ws` | `Vector` | Yes | Barycentric weights (alternative to node-only) |
| `f` | `Function` | Yes | Function to sample at nodes |

**Returns:** `Lagrange` polynomial

**Examples:**

```julia
# Direct interpolation of points
xs = [0.0, 0.5, 1.0]
ys = [0.0, 0.25, 1.0]  # y = x²
p = fit(Lagrange, xs, ys)

# Interpolation of a function
f(x) = sin(π*x)
xs = range(-1, 1, length=11)
p = fit(Lagrange, xs, f.(xs))
```
```

--------------------------------

### Interpolating polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/core-types.md

Constructors for Lagrange and Newton interpolating polynomials.

```APIDOC
### Interpolating polynomials

#### Lagrange{N, S, R, T, X}

```julia
Lagrange(xs::Vector{S}, ws::Vector{R}, coeffs::Vector{T}, var::Symbol=:x)
Lagrange(xs::Vector{S}, coeffs::Vector{T}, var::Symbol=:x)  # computes weights
```

Lagrange interpolation over nodes `xs` with barycentric weights `ws`.

#### Newton{N, S, T, X}

```julia
Newton(xs::Vector{S}, coeffs::Vector{T}, var::Symbol=:x)
```

Newton divided-difference interpolation over nodes `xs`.
```

--------------------------------

### Differentiate and Integrate Orthogonal Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/orthogonal-polynomials.md

Illustrates how to compute the derivative and integral of an orthogonal polynomial. The results are returned as standard Polynomial types.

```julia
p = Legendre([1, 2, 3])

dp = derivative(p)      # returns Polynomial
ip = integrate(p)       # returns Polynomial
```

--------------------------------

### Verifying Orthogonality Numerically

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/index.md

Illustrates how to numerically verify the orthogonality property of polynomial families using basis functions. The conceptual integral relation is provided as a comment.

```julia
# Not exported, but conceptually:
# ∫_{-1}^{1} Pₙ(x) Pₘ(x) dx = cₙ δₙₘ

# Verify numerically
p = basis(Legendre, 3)
q = basis(Legendre, 5)
```

--------------------------------

### High-Degree Interpolation with Barycentric Lagrange

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/getting-started.md

Demonstrates stable high-degree interpolation using Barycentric Lagrange method with Chebyshev nodes and weights.

```julia
# Good nodes + barycentric formula = stable interpolation
xs, ws = SpecialPolynomials.lagrange_barycentric_nodes_weights(
    SpecialPolynomials.ChebyshevBasis, 100
)
p = fit(Lagrange, xs, ws, f)

# Now high-degree interpolation is stable
p(0.5)
```

--------------------------------

### Convert Discrete Polynomial to Standard Polynomial and Evaluate

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/discrete-orthogonal.md

Demonstrates converting a Hahn polynomial to a standard polynomial representation and then evaluating it at real numbers, enabling interpolation.

```julia
h = HahnPolynomial{1, 1, 10}([1, 2, 3])

# Convert to standard polynomial
poly = convert(Polynomial, h)

# Evaluate at real numbers (extends polynomial)
poly(2.5)  # interpolation outside original points
```

--------------------------------

### Compute Inner Product of Legendre Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/docs/src/examples.md

Demonstrates computing the inner product of two Legendre polynomials using numerical integration (quadgk) and the unexported `innerproduct` function. The inner product of orthogonal polynomials with different degrees should be zero.

```jldoctest
julia> using QuadGK


julia> P = Legendre
Polynomials.MutableDensePolynomial{SpecialPolynomials.LegendreBasis}

julia> p4,p5 = basis.(P, [4,5])
2-element Vector{Polynomials.MutableDensePolynomial{SpecialPolynomials.LegendreBasis, Float64, :x}}:
 Legendre(1.0⋅P₄(x))
 Legendre(1.0⋅P₅(x))

julia> wf, dom = SpecialPolynomials.weight_function(P), Polynomials.domain(P);


julia> quadgk(x -> p4(x) * p5(x) *  wf(x), first(dom), last(dom))
(0.0, 0.0)

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

--------------------------------

### Probabilists' Hermite Polynomial Initialization

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/orthogonal-polynomials.md

Initializes a Probabilists' Hermite polynomial. This variant is related to the Physicists' Hermite polynomials by a scaling factor and argument transformation.

```julia
# Probabilists' Hermite: related by Hᵢ(x) = 2^(i/2) Heᵢ(√2x)
h_prob = ChebyshevHermite([1, 2, 3])
```

--------------------------------

### Compare Root Finding Stability for High-Degree Legendre Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/docs/src/examples.md

Illustrates the difference in stability between using the companion matrix and the `jacobi_matrix` for finding roots of a high-degree Legendre polynomial. The `jacobi_matrix` method is more accurate.

```julia
julia> using LinearAlgebra

julia> p50 = basis(Legendre{Float64}, 50)
Legendre(1.0⋅P₅₀(x))

julia> as = eigvals(Polynomials.companion(p50));

julia> maximum(abs ∘ p50, as) < sqrt(eps())
false

julia> bs = eigvals(SpecialPolynomials.jacobi_matrix(Legendre, 50 ));

julia> maximum(abs ∘ p50, bs) < sqrt(eps())
true

julia> maximum(abs, roots(p50) - bs) < sqrt(eps())
true
```

--------------------------------

### Convert Newton Polynomial to Standard Polynomial

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/interpolating-polynomials.md

Illustrates converting a Newton polynomial representation to a standard polynomial form.

```julia
p = Newton([1, 2, 3], [1, 2, 1])
poly = convert(Polynomial, p)  # to standard polynomial
```

--------------------------------

### Single-parameter Polynomial Type Constructors

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/core-types.md

Provides generic construction signatures for single-parameter polynomial types like Chebyshev, Legendre, Hermite, Laguerre, and Bessel polynomials.

```julia
P(coeffs::Vector{T}, var::Symbol=:x)
P(coeffs::Tuple, var::Symbol=:x)
```

--------------------------------

### Relate `fromroots` and `roots` for Jacobi Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/docs/src/examples.md

Shows the relationship between `fromroots` and `roots` for Jacobi polynomials, requiring a monic polynomial. This is useful for reconstructing polynomials from their roots.

```julia
julia> using Polynomials, SpecialPolynomials; const SP=SpecialPolynomials
SpecialPolynomials

julia> P = Jacobi{1/2,-1/2}
Polynomials.MutableDensePolynomial{SpecialPolynomials.JacobiBasis{0.5, -0.5}}


julia> p = P([1,1,2,3])
Jacobi{0.5,-0.5}(1⋅Jᵅᵝ₀(x) + 1⋅Jᵅᵝ₁(x) + 2⋅Jᵅᵝ₂(x) + 3⋅Jᵅᵝ₃(x))

julia> q = SP.monic(p) # monic is not exported
Jacobi{0.5,-0.5}(0.13333333333333333⋅Jᵅᵝ₀(x) + 0.13333333333333333⋅Jᵅᵝ₁(x) + 0.26666666666666666⋅Jᵅᵝ₂(x) + 0.4⋅Jᵅᵝ₃(x))

julia> fromroots(P, roots(q)) - q |> u -> truncate(u, atol=sqrt(eps()))
Jacobi{0.5,-0.5}(0.0 + 0.0im)
```

--------------------------------

### Lagrange Polynomial Construction with Automatic Weight Computation

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/interpolating-polynomials.md

Constructs a Lagrange polynomial where barycentric weights are computed automatically from the provided nodes. This simplifies the construction process.

```julia
Lagrange(xs::Vector, coeffs::Vector, var::Symbol=:x)
```

```julia
xs = [1.0, 2.0, 3.0]
ys = [1.0, 4.0, 9.0]

p = Lagrange(xs, ys)  # weights computed automatically
p(1.5)
```

--------------------------------

### Access Weight Functions for Orthogonal Polynomial Bases

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/orthogonal-polynomials.md

Demonstrates how to access the weight functions associated with different orthogonal polynomial basis types. These functions define the inner product for the polynomials.

```julia
# Access via the basis type
weight_function(LegendreBasis)     # x -> 1
weight_function(ChebyshevBasis)    # x -> 1/√(1-x²)
weight_function(HermiteBasis)      # x -> exp(-x²)
```

--------------------------------

### Convert Orthogonal Polynomials to Standard Polynomials

Source: https://github.com/jverzani/specialpolynomials.jl/blob/master/_autodocs/api-reference/orthogonal-polynomials.md

Demonstrates converting an orthogonal polynomial to a standard Polynomial type and back. Intermediate type conversions may be necessary.

```julia
p = Chebyshev([1, 2, 3])

# Convert to standard polynomial
poly = convert(Polynomial, p)

# Convert back (through intermediate type if needed)
p2 = convert(Chebyshev, convert(Polynomial, p))
```