WIP

Author: jverzani

README.md

@@ -96,8 +96,8 @@ ERROR: Polynomials must have same variable.
 #### Integrals and Derivatives
 
 Integrate the polynomial `p` term by term, optionally adding constant
-term `k`. The order of the resulting polynomial is one higher than the
-order of `p`.
+term `k`. The degree of the resulting polynomial is one higher than the
+degree of `p`.
 
 ```julia
 julia> integrate(Polynomial([1, 0, -1]))
@@ -107,8 +107,8 @@ julia> integrate(Polynomial([1, 0, -1]), 2)
 Polynomial(2.0 + x - 0.3333333333333333x^3)
 ```
 
-Differentiate the polynomial `p` term by term. The order of the
-resulting polynomial is one lower than the order of `p`.
+Differentiate the polynomial `p` term by term. The degree of the
+resulting polynomial is one lower than the degree of `p`.
 
 ```julia
 julia> derivative(Polynomial([1, 3, -1]))
@@ -119,7 +119,7 @@ Polynomial(3 - 2x)
 
 
 Return the roots (zeros) of `p`, with multiplicity. The number of
-roots returned is equal to the order of `p`. By design, this is not type-stable,
+roots returned is equal to the degree of `p`. By design, this is not type-stable,
 the returned roots may be real or complex.
 
 ```julia
@@ -141,7 +141,7 @@ julia> roots(Polynomial([0, 0, 1]))
 
 #### Fitting arbitrary data
 
-Fit a polynomial (of order `deg`) to `x` and `y` using a least-squares approximation.
+Fit a polynomial (of degree `deg`) to `x` and `y` using a least-squares approximation.
 
 ```julia
 julia> xs = 0:4; ys = @. exp(-xs) + sin(xs);
@@ -174,7 +174,7 @@ Polynomial objects also have other methods:
 
 * `conj`: finds the conjugate of a polynomial over a complex fiel
 
-* `truncate`: set to 0 all small terms in a polynomial; 
+* `truncate`: set to 0 all small terms in a polynomial;
 * `chop` chops off any small leading values that may arise due to floating point operations.
 
 * `gcd`: greatest common divisor of two polynomials.

docs/src/extending.md

@@ -1,12 +1,15 @@
 # Extending Polynomials
 
-The [`AbstractPolynomial`](@ref) type was made to be extended via a rich interface. 
+The [`AbstractPolynomial`](@ref) type was made to be extended via a rich interface.
 
 ```@docs
 AbstractPolynomial
 ```
 
-To implement a new polynomial type, `P`, the following methods should be implemented. 
+A polynomial's  coefficients  are  relative to some *basis*. The `Polynomial` type relates coefficients  `[a0, a1,  ..., an]`, say,  to the  polynomial  `a0 +  a1*x + a2*x^  + ... +  an*x^n`,  through the natural  basis  `1,  x,  x^2, ..., x^n`.  New p olynomial  types typically represent the polynomial through a different  basis. For example,  `CheyshevT` uses a basis  `T_0=1, T_1=x,  T_2=2x^2-1,  ...,  T_n  =  2xT_{n-1} - T_{n-2}`.  For this type  the  coefficients  `[a0,a1,...,an]` are associated with  the polynomial  `a0*T0  + a1*T_1 +  ...  +  an*T_n`.
+
+To implement a new polynomial type, `P`, the following methods should
+be implemented.
 
 !!! note
     Promotion rules will always coerce towards the [`Polynomial`](@ref) type, so not all methods have to be implemented if you provide a conversion function.
@@ -26,5 +29,6 @@ As always, if the default implementation does not work or there are more efficie
 | `-(::P, ::P)` | | Subtraction of polynomials |
 | `*(::P, ::P)` | | Multiplication of polynomials |
 | `divrem` | | Required for [`gcd`](@ref)|
+| `variable`| | Convenience to find monomial `x` in new  basis|
 
-Check out both the [`Polynomial`](@ref) and [`ChebyshevT`](@ref) for examples of this interface being extended!
+Check out both the [`Polynomial`](@ref) and [`ChebyshevT`](@ref) for examples of this interface being extended. [`Bernstein`](@ref) is an example where the basis depends on the  degree of  the polynomials being represented.

docs/src/index.md

@@ -163,14 +163,17 @@ savefig("polyfit.svg"); nothing # hide
 
 ![](polyfit.svg)
 
+### Iteration
+
+A polynomial, e.g. `a_0 + a_1 x + a_2 x^2 + ... + a_n x^n` can be seen as a collection of coefficients, `[a_0, a_1, ..., a_n]`, reelative to some polynomial basis. The most  familiar basis being  `1`, `x`, `x^2`, ... If the basis is implicit, then a polynomial is just a vector. Vectors or 1-based, but polynomial types are 0-based, for purposes of indexing (e.g. `getindex`, `setindex!`, `eachindex`). Iteration over a polynomial steps through the basis vectors, e.g. `a_0`, `a_1 x`, ...
 
 ## Related Packages
 
 * [MultiPoly.jl](https://github.com/daviddelaat/MultiPoly.jl) for sparse multivariate polynomials
 
 * [MultivariatePolynomials.jl](https://github.com/blegat/MultivariatePolynomials.jl) for multivariate polynomials and moments of commutative or non-commutative variables
 
-* [Nemo.jl](https://github.com/wbhart/Nemo.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series
+* [AbstractAlgeebra.jl](https://github.com/wbhart/AbstractAlgebra.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series.
 
 * [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) for a fast complex polynomial root finder
 

docs/src/polynomials/bernstein.md

@@ -0,0 +1,45 @@
+# Bernstein Polynomials
+
+```@meta
+DocTestSetup = quote
+  using Polynomials
+end
+```
+
+
+The [Bernstein polynomials](https://en.wikipedia.org/wiki/Bernstein_polynomial) are a family of polynomials defined on the interval `[0,1]`. For each `n` there are `n+1` polynomials, given by: `b(n,nu) = choose(n,nu) x^nu * (1-x)^(n-nu)` for `nu` in `0:n`. Together, these form a basis that can represent any polynomial of degree `n` or less through a linear combination.
+
+
+## Bernstein(N,T)
+
+```@docs
+Bernstein
+Bernstein{N,T}()
+```
+
+The `Bernstein{N,T}` type holds coefficients representing the polynomial `a_0 b(n,0) + a_1 b(n,1) + ... + a_n b(n,n)`.
+
+For example, using `n=3`, the monomial `x` is represented by `0 b(3,0) + 1/3 b(3,1) + 2/3 b(3,2) + 1  b(3,3)`, which can be constructed through `Bernstein{3, T}([0, 1/3, 2/3, 1])`
+
+
+### Conversion
+
+[`Bernsteing`](@ref) can be converted to [`Polynomial`](@ref) and vice-versa.
+
+```jldoctest
+julia> b = Bernstein([1, 0, 3, 4])
+Bernstein(1⋅β(3, 0)(x) + 3⋅β(3, 2)(x) + 4⋅β(3, 3)(x))
+
+julia> p = convert(Polynomial{Int}, b)
+Polynomial(1 - 3*x + 12*x^2 - 6*x^3)
+
+julia> convert(Bernstein{3, Float64},  p)
+Bernstein(1.0⋅β(3, 0)(x) + 3.0⋅β(3, 2)(x) + 4.0⋅β(3, 3)(x))
+```
+
+Bernstein polynomials may be  converted  into a basis with a larger  degree:
+
+```jldoctest
+julia> convert(Bernstein{4, Float64}, b)
+Bernstein(1.0⋅β(4, 0)(x) + 0.25⋅β(4, 1)(x) + 1.5⋅β(4, 2)(x) + 3.25⋅β(4, 3)(x) + 4.0⋅β(4, 4)(x))
+```

docs/src/polynomials/chebyshev.md

@@ -6,24 +6,34 @@ DocTestSetup = quote
 end
 ```
 
+
+The [Chebyshev polynomials](https://en.wikipedia.org/wiki/Chebyshev_polynomials) are two sequences of polynomials, `T_n` and `U_n`. The Chebyshev polynomials of the first kind, `T_n`, can be defined by the recurrence relation `T_0(x)=1`, `T_1(x)=x`, and `T_{n+1}(x) = 2xT_n{x}-T_{n-1}(x)`. The Chebyshev polynomioals of the second kind, `U_n(x)`, can be defined by `U_0(x)=1`, `U_1(x)=2x`, and `U_{n+1}(x) = 2xU_n(x) - U_{n-1}(x)`. Both `T_n` and `U_n` have degree `n`, and any polynomial  of  degree `n` may be written as a linear combination of the polynomials `T_0`, `T_1`, ..., `T_n` (similarly with `U`).
+
+
 ## First Kind
 
 ```@docs
 ChebyshevT
 ChebyshevT()
 ```
 
+The `ChebyshevT` type holds coefficients representing the polynomial `a_0 T_0 + a_1 T_1 + ... + a_n T_n`.
+
+For example, the basis polynomial `T_4` can be represented with `ChebyshevT([0,0,0,0,1])`.
+
+
 ### Conversion
 
 [`ChebyshevT`](@ref) can be converted to [`Polynomial`](@ref) and vice-versa.
 
 ```jldoctest
 julia> c = ChebyshevT([1, 0, 3, 4])
-ChebyshevT([1, 0, 3, 4])
+ChebyshevT(1⋅T_0(x) + 3⋅T_2(x) + 4⋅T_3(x))
+
 
 julia> p = convert(Polynomial, c)
 Polynomial(-2 - 12*x + 6*x^2 + 16*x^3)
 
 julia> convert(ChebyshevT, p)
-ChebyshevT([1.0, 0.0, 3.0, 4.0])
+ChebyshevT(1.0⋅T_0(x) + 3.0⋅T_2(x) + 4.0⋅T_3(x))
 ```

docs/src/reference.md

@@ -93,17 +93,19 @@ will plot the polynomial within the range `[a, b]`.
 ### Example: The Polynomials.jl logo
 ```@example
 using Plots, Polynomials
-xs = range(-1, 1, length=100)
+# T1, T2, T3, and T4:
 chebs = [
   ChebyshevT([0, 1]),
   ChebyshevT([0, 0, 1]),
   ChebyshevT([0, 0, 0, 1]),
   ChebyshevT([0, 0, 0, 0, 1]),
 ]
 colors = ["#4063D8", "#389826", "#CB3C33", "#9558B2"]
-plot() # hide
-for (cheb, col) in zip(chebs, colors)
-  plot!(xs, cheb.(xs), c=col, lw=5, label="")
+itr = zip(chebs, colors)
+(cheb,col), state = iterate(itr)
+p = plot(cheb, c=col,  lw=5, legend=false, label="")
+for (cheb, col) in Base.Iterators.rest(itr, state)
+  plot!(cheb, c=col, lw=5)
 end
 savefig("chebs.svg"); nothing # hide
 ```

src/Polynomials.jl

@@ -12,6 +12,7 @@ include("contrib.jl")
 include("polynomials/Polynomial.jl")
 include("polynomials/ChebyshevT.jl")
 include("polynomials/ChebyshevU.jl")
+include("polynomials/Bernstein.jl")
 
 include("polynomials/Poly.jl") # Deprecated -> Will be removed
 include("pade.jl")