Merge pull request #196 from jverzani/miles-lucas-refactor-delta

Miles lucas refactor delta

Author: jverzani

.travis.yml

@@ -9,6 +9,7 @@ julia:
   - 1.1
   - 1.2
   - 1.3
+  - 1.4
   - nightly
 
 matrix:

Project.toml

@@ -10,8 +10,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 
 [compat]
-Intervals = "0.5.0"
-RecipesBase = "0.7, 0.8"
+Intervals = "0.5.0, 1.0"
+RecipesBase = "0.7, 0.8, 1.0"
 julia = "1"
 
 

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
@@ -130,8 +130,8 @@ julia> roots(Polynomial([1, 0, -1]))
 
 julia> roots(Polynomial([1, 0, 1]))
 2-element Array{Complex{Float64},1}:
- 0.0+1.0im
- 0.0-1.0im
+ 0.0 - 1.0im
+ 0.0 + 1.0im
 
 julia> roots(Polynomial([0, 0, 1]))
 2-element Array{Float64,1}:
@@ -141,16 +141,16 @@ 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);
 
-julia> fit(xs, ys)
-Polynomial(1.0000000000000016 + 0.059334723072240664*x + 0.39589720602859824*x^2 - 0.2845598112184312*x^3 + 0.03867830809692903*x^4)
+julia> fit(xs, ys) |> x -> round(x, digits=4)
+Polynomial(1.0 + 0.0593*x + 0.3959*x^2 - 0.2846*x^3 + 0.0387*x^4)
 
-julia> fit(ChebyshevT, xs, ys, deg=2)
-ChebyshevT([0.541280671210034, -0.8990834124779993, -0.4237852336242923])
+julia> fit(ChebyshevT, xs, ys, deg=2) |> x -> round(x, digits=4)
+ChebyshevT(0.5413⋅T_0(x) - 0.8991⋅T_1(x) - 0.4238⋅T_2(x))
 ```
 
 Visual example:
@@ -166,15 +166,16 @@ Polynomial objects also have other methods:
 
 * `coeffs`: returns the entire coefficient vector
 
-* `degree`: returns the polynomial degree, `length` is 1 plus the degree
+* `degree`: returns the polynomial degree, `length` is number of stored coefficients
 
-* `variable`: returns the polynomial symbol as a degree 1 polynomial
+* `variable`: returns the polynomial symbol as polynomial in the underlying type
 
 * `norm`: find the `p`-norm of a polynomial
 
-* `conj`: finds the conjugate of a polynomial over a complex fiel
+* `conj`: finds the conjugate of a polynomial over a complex field
+
+* `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 standard  basis  `1,  x,  x^2, ..., x^n`.  New polynomial  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.

docs/src/index.md

@@ -104,8 +104,8 @@ ERROR: Polynomials must have same variable
 ### Integrals and Derivatives
 
 Integrate the polynomial `p` term by term, optionally adding constant
-term `C`. The order of the resulting polynomial is one higher than the
-order of `p`.
+term `C`. The degree of the resulting polynomial is one higher than the
+degree of `p`.
 
 ```jldoctest
 julia> integrate(Polynomial([1, 0, -1]))
@@ -115,8 +115,8 @@ julia> integrate(Polynomial([1, 0, -1]), 2)
 Polynomial(2.0 + 1.0*x - 0.3333333333333333*x^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`.
 
 ```jldoctest
 julia> derivative(Polynomial([1, 3, -1]))
@@ -137,8 +137,8 @@ julia> roots(Polynomial([1, 0, -1]))
 
 julia> roots(Polynomial([1, 0, 1]))
 2-element Array{Complex{Float64},1}:
- -0.0 + 1.0im
   0.0 - 1.0im
+  0.0 + 1.0im
 
 julia> roots(Polynomial([0, 0, 1]))
 2-element Array{Float64,1}:
@@ -148,13 +148,13 @@ 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 polynomial interpolation or a (weighted) least-squares approximation.
 
 ```@example
 using Plots, Polynomials
 xs = range(0, 10, length=10)
 ys = exp.(-xs)
-f = fit(xs, ys)
+f = fit(xs, ys)  # fit(xs, ys, k)  for fitting a kth degreee polynomial
 
 scatter(xs, ys, label="Data");
 plot!(f, extrema(xs)..., label="Fit");
@@ -163,14 +163,55 @@ savefig("polyfit.svg"); nothing # hide
 
 ![](polyfit.svg)
 
+### Other bases
+
+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]`, relative to some polynomial basis. The most  familiar basis being  the standard one: `1`, `x`, `x^2`, ...  Alternative bases are possible.  The `ChebyshevT` polynomials are  implemented, as an example. Instead of `Polynomial`  or  `Polynomial{T}`, `ChebyshevT` or  `ChebyshevT{T}` constructors are used:
+
+```jldoctest
+julia> p1 = ChebyshevT([1.0, 2.0, 3.0])
+ChebyshevT(1.0⋅T_0(x) + 2.0⋅T_1(x) + 3.0⋅T_2(x))
+
+julia> p2 = ChebyshevT{Float64}([0, 1, 2])
+ChebyshevT(1.0⋅T_1(x) + 2.0⋅T_2(x))
+
+julia> p1 + p2
+ChebyshevT(1.0⋅T_0(x) + 3.0⋅T_1(x) + 5.0⋅T_2(x))
+
+julia> p1 * p2
+ChebyshevT(4.0⋅T_0(x) + 4.5⋅T_1(x) + 3.0⋅T_2(x) + 3.5⋅T_3(x) + 3.0⋅T_4(x))
+
+julia> derivative(p1)
+ChebyshevT(2.0⋅T_0(x) + 12.0⋅T_1(x))
+
+julia> integrate(p2)
+ChebyshevT(0.25⋅T_0(x) - 1.0⋅T_1(x) + 0.25⋅T_2(x) + 0.3333333333333333⋅T_3(x))
+
+julia> convert(Polynomial, p1)
+Polynomial(-2.0 + 2.0*x + 6.0*x^2)
+
+julia> convert(ChebyshevT, Polynomial([1.0, 2,  3]))
+ChebyshevT(2.5⋅T_0(x) + 2.0⋅T_1(x) + 1.5⋅T_2(x))
+```
+
+
+### Iteration
+
+If its basis is implicit, then a polynomial may be  seen as just a vector of  coefficients. Vectors or 1-based, but, for convenience, 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`, ...
+
+```jldoctest
+julia> as = [1,2,3,4,5]; p = Polynomial(as);
+
+julia> as[3], p[2], collect(p)[3]
+(3, 3, Polynomial(3*x^2))
+```
 
 ## 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/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 uniquely 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

@@ -16,14 +16,14 @@ end
 
 ```@docs
 coeffs
-order
 degree
 length
 size
 domain
 mapdomain
 chop
 chop!
+round
 truncate
 truncate!
 ```
@@ -91,19 +91,22 @@ plot(::AbstractPolynomial, a, b; kwds...)
 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

@@ -8,16 +8,15 @@ include("show.jl")
 include("plots.jl")
 include("contrib.jl")
 
+# Interface for all AbstractPolynomials
+include("common.jl")
+
+
 # Polynomials
 include("polynomials/Polynomial.jl")
 include("polynomials/ChebyshevT.jl")
-include("polynomials/ChebyshevU.jl")
 
-include("polynomials/Poly.jl") # Deprecated -> Will be removed
-include("pade.jl")
-include("compat.jl") # Where we keep deprecations
-
-# Interface for all AbstractPolynomials
-include("common.jl")
+# to be deprecated, then removed
+include("compat.jl")
 
 end # module

src/abstract.jl

@@ -35,7 +35,7 @@ macro register(name)
     poly = esc(name)
     quote
         Base.convert(::Type{P}, p::P) where {P<:$poly} = p
-        Base.convert(P::Type{<:$poly}, p::$poly) where {T} = P(p.coeffs, p.var)
+        Base.convert(P::Type{<:$poly}, p::$poly) where {T} = P(coeffs(p), p.var)
         Base.promote_rule(::Type{$poly{T}}, ::Type{$poly{S}}) where {T,S} =
             $poly{promote_type(T, S)}
         Base.promote_rule(::Type{$poly{T}}, ::Type{S}) where {T,S<:Number} =