Issue 312a (#328)

* add some references; close issue #311

* close issue #312

* merge in master; better fix for issue #320

* update doctests for v1.6 change to print of Vectors

Author: jverzani

README.md

@@ -206,15 +206,20 @@ Polynomial objects also have other methods:
 
 * [MultivariatePolynomials.jl](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl) for multivariate polynomials and moments of commutative or non-commutative variables
 
-* [PolynomialRings](https://github.com/tkluck/PolynomialRings.jl) A library for arithmetic and algebra with multi-variable polynomials.
+* [PolynomialRings.jl](https://github.com/tkluck/PolynomialRings.jl) A library for arithmetic and algebra with multi-variable polynomials.
 
 * [AbstractAlgebra.jl](https://github.com/wbhart/AbstractAlgebra.jl), [Nemo.jl](https://github.com/wbhart/Nemo.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series, [Hecke.jl](https://github.com/thofma/Hecke.jl) for algebraic number theory.
 
-* [CommutativeAlgebra](https://github.com/KlausC/CommutativeRings.jl) the start of a computer algebra system specialized to discrete calculations with support for polynomials.
+* [CommutativeAlgebra.jl](https://github.com/KlausC/CommutativeRings.jl) the start of a computer algebra system specialized to discrete calculations with support for polynomials.
 
 * [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) for a fast complex polynomial root finder. For larger degree problems, also [FastPolynomialRoots](https://github.com/andreasnoack/FastPolynomialRoots.jl) and [AMRVW](https://github.com/jverzani/AMRVW.jl).
 
 
+* [SpecialPolynomials.jl](https://github.com/jverzani/SpecialPolynomials.jl) A package providing various polynomial types beyond the standard basis polynomials in `Polynomials.jl`. Includes interpolating polynomials, Bernstein polynomials, and classical orthogonal polynomials.
+
+* [ClassicalOrthogonalPolynomials.jl](https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl) A Julia package for classical orthogonal polynomials and expansions. Includes `chebyshevt`, `chebyshevu`, `legendrep`, `jacobip`, `ultrasphericalc`, `hermiteh`, and `laguerrel`. The same repository includes `FastGaussQuadrature.jl`, `FastTransforms.jl`, and the `ApproxFun` packages.
+
+
 ## Legacy code
 
 As of v0.7, the internals of this package were greatly generalized and new types and method names were introduced. For compatability purposes, legacy code can be run after issuing `using Polynomials.PolyCompat`.

docs/src/extending.md

@@ -25,7 +25,6 @@ As always, if the default implementation does not work or there are more efficie
 | `domain` | x | Should return an  [`AbstractInterval`](https://invenia.github.io/Intervals.jl/stable/#Intervals-1) |
 | `vander` | | Required for [`fit`](@ref) |
 | `companion` | | Required for [`roots`](@ref) |
-| `fromroots` | | By default, will form polynomials using `prod(variable(::P) - r)` for reach root `r`|
 | `*(::P, ::P)` | | Multiplication of polynomials |
 | `divrem` | | Required for [`gcd`](@ref)|
 | `one`| | Convenience to find constant in new basis |
@@ -141,7 +140,7 @@ Now `p` is treated as the vector `p.coeffs`, as regards broadcasting, so some th
 
 ```jldoctest AliasPolynomial
 julia> p .+ 2
-4-element Array{Int64,1}:
+4-element Vector{Int64}:
  5
  4
  5

docs/src/index.md

@@ -144,17 +144,17 @@ the returned roots may be real or complex.
 
 ```jldoctest
 julia> roots(Polynomial([1, 0, -1]))
-2-element Array{Float64,1}:
+2-element Vector{Float64}:
  -1.0
   1.0
 
 julia> roots(Polynomial([1, 0, 1]))
-2-element Array{Complex{Float64},1}:
+2-element Vector{ComplexF64}:
  0.0 - 1.0im
  0.0 + 1.0im
 
 julia> roots(Polynomial([0, 0, 1]))
-2-element Array{Float64,1}:
+2-element Vector{Float64}:
  0.0
  0.0
 ```
@@ -227,7 +227,7 @@ The `pairs` iterator, iterates over the indices and coefficients, attempting to
 
 ```jldoctest
 julia> v = [1,2,0,4]
-4-element Array{Int64,1}:
+4-element Vector{Int64}:
  1
  2
  0
@@ -236,7 +236,7 @@ julia> v = [1,2,0,4]
 julia> p,ip,sp,lp = Polynomial(v), ImmutablePolynomial(v), SparsePolynomial(v), LaurentPolynomial(v, -1);
 
 julia> collect(pairs(p))
-4-element Array{Pair{Int64,Int64},1}:
+4-element Vector{Pair{Int64, Int64}}:
  0 => 1
  1 => 2
  2 => 0
@@ -246,13 +246,13 @@ julia> collect(pairs(ip)) == collect(pairs(p))
 true
 
 julia> collect(pairs(sp)) # unordered dictionary with only non-zero terms
-3-element Array{Pair{Int64,Int64},1}:
+3-element Vector{Pair{Int64, Int64}}:
  0 => 1
  3 => 4
  1 => 2
 
 julia> collect(pairs(lp))
-4-element Array{Pair{Int64,Int64},1}:
+4-element Vector{Pair{Int64, Int64}}:
  -1 => 1
   0 => 2
   1 => 0
@@ -267,13 +267,126 @@ julia> p = Polynomial([1,2,0,4], :u)
 Polynomial(1 + 2*u + 4*u^3)
 
 julia> collect(Polynomials.monomials(p))
-4-element Array{Any,1}:
+4-element Vector{Any}:
  Polynomial(1)
  Polynomial(2*u)
  Polynomial(0)
  Polynomial(4*u^3)
 ```
 
+## Relationship between the `T` and `P{T,X}`
+
+The addition of a polynomial and a scalar, such as
+
+```jldoctest natural_inclusion
+julia> using Polynomials
+
+julia> p = Polynomial([1,2,3], :x)
+Polynomial(1 + 2*x + 3*x^2)
+
+julia> p + 3
+Polynomial(4 + 2*x + 3*x^2)
+```
+
+seems natural, but in `Julia`, as `3` is of type `Int` and `p` of type `Polynomial{Int,:x}` some addition must be defined. The basic idea  is that `3` is promoted to the *constant* polynomial `3` with indeterminate `:x` as `3*one(p)` and then addition of `p + 3*one(p)` is performed.
+
+This identification of a scalar with a constant polynomial can go both ways. If `q` is a *constant* polynomial of type `Polynomial{Int, :y}` then we should expect that `p+q` would be defined, as `p` plus the constant term of `q`. Indeed this is the case
+
+```jldoctest natural_inclusion
+julia> q = Polynomial(3, :y)
+Polynomial(3)
+
+julia> p + q
+Polynomial(4 + 2*x + 3*x^2)
+```
+
+If `q` is non-constant, such as `variable(Polynomial, :y)`, then there would be an error due to the mismatched symbols. (The mathematical result would need a multivariate polynomial, not a univariate polynomial, as this package provides.)
+
+The same conversion is done for polynomial multiplication: constant polynomials are treated as numbers; non-constant polynomials must have their symbols match. 
+
+There is an oddity though the following two computations look the same, they are technically different:
+
+```jldoctest natural_inclusion
+julia> one(Polynomial, :x) + one(Polynomial, :y)
+Polynomial(2.0)
+
+julia> one(Polynomial, :y) + one(Polynomial, :x)
+Polynomial(2.0)
+```
+
+Both are constant polynomials over `Int`, but the first has the indeterminate `:y`, the second `:x`. 
+
+This technical difference causes no issues with polynomial addition or multiplication, as there constant polynomials are treated as numbers, but can be an issue when constant polynomials are used as array elements.
+
+For arrays, the promotion of numbers to polynomials, allows natural constructions like:
+
+```jldoctest natural_inclusion
+julia> p = Polynomial([1,2],:x)
+Polynomial(1 + 2*x)
+
+julia> q = Polynomial([1,2],:y)  # non-constant polynomials with different indeterminates
+Polynomial(1 + 2*y)
+
+julia> [1 p]
+1×2 Matrix{Polynomial{Int64, :x}}:
+ Polynomial(1)  Polynomial(1 + 2*x)
+
+julia> [1 one(q)]
+1×2 Matrix{Polynomial{Int64, :y}}:
+ Polynomial(1)  Polynomial(1)
+```
+
+However, as there would be an ambiguous outcome of the following
+
+```jldoctest natural_inclusion
+julia> [one(p) one(q)]
+ERROR: ArgumentError: Polynomials have different indeterminates
+[...]
+```
+
+an error thrown.
+
+In general, arrays with mixtures of non-constant polynomials with *different* indeterminates will error. By default, an error will occur when constant polynomials with different indeterminates are used as components. However, for *typed* arrays, conversion will allow such constructs to be used.
+
+Using `one(q)` for a constant polynomial with indeterminate `:y` we have:
+
+```jldoctest natural_inclusion
+julia> P = typeof(p)
+Polynomial{Int64, :x}
+
+julia> P[one(p) one(q)]
+1×2 Matrix{Polynomial{Int64, :x}}:
+ Polynomial(1)  Polynomial(1)
+```
+
+Of course, by not being explicit, there are sill gotchas. For example, we can construct this matrix without a specific types:
+
+```jldoctest natural_inclusion
+julia> [one(p), one(q)+one(p)]
+2-element Vector{Polynomial{Int64, :x}}:
+ Polynomial(1)
+ Polynomial(2)
+```
+
+but not this one:
+
+```jldoctest natural_inclusion
+julia> [one(p), one(p) + one(q)]
+ERROR: ArgumentError: Polynomials have different indeterminates
+[...]
+```
+
+Also, mixing types can result in unspecific symbols, as this example shows:
+
+```jldoctest natural_inclusion
+julia> [1 p; p 1] + [1 2one(q); 3 4] # array{P{T,:x}} + array{P{T,:y}}
+2×2 Matrix{Polynomial{Int64, X} where X}:
+ Polynomial(2)        Polynomial(3 + 2*x)
+ Polynomial(4 + 2*x)  Polynomial(5)
+```
+
+Though were a non-constant polynomial with indeterminate `y` replacing
+`2one(q)` above, that addition would throw an error.
 
 ## Related Packages
 
@@ -285,14 +398,17 @@ julia> collect(Polynomials.monomials(p))
 
 * [MultivariatePolynomials.jl](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl) for multivariate polynomials and moments of commutative or non-commutative variables
 
-* [PolynomialRings](https://github.com/tkluck/PolynomialRings.jl) A library for arithmetic and algebra with multi-variable polynomials.
+* [PolynomialRings.jl](https://github.com/tkluck/PolynomialRings.jl) A library for arithmetic and algebra with multi-variable polynomials.
 
 * [AbstractAlgebra.jl](https://github.com/wbhart/AbstractAlgebra.jl), [Nemo.jl](https://github.com/wbhart/Nemo.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series, [Hecke.jl](https://github.com/thofma/Hecke.jl) for algebraic number theory.
 
-* [CommutativeAlgebra](https://github.com/KlausC/CommutativeRings.jl) the start of a computer algebra system specialized to discrete calculations with support for polynomials.
+* [CommutativeAlgebra.jl](https://github.com/KlausC/CommutativeRings.jl) the start of a computer algebra system specialized to discrete calculations with support for polynomials.
 
 * [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) for a fast complex polynomial root finder. For larger degree problems, also [FastPolynomialRoots](https://github.com/andreasnoack/FastPolynomialRoots.jl) and [AMRVW](https://github.com/jverzani/AMRVW.jl).
 
+* [SpecialPolynomials.jl](https://github.com/jverzani/SpecialPolynomials.jl) A package providing various polynomial types beyond the standard basis polynomials in `Polynomials.jl`. Includes interpolating polynomials, Bernstein polynomials, and classical orthogonal polynomials.
+
+* [ClassicalOrthogonalPolynomials.jl](https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl) A Julia package for classical orthogonal polynomials and expansions. Includes `chebyshevt`, `chebyshevu`, `legendrep`, `jacobip`, `ultrasphericalc`, `hermiteh`, and `laguerrel`. The same repository includes `FastGaussQuadrature.jl`, `FastTransforms.jl`, and the `ApproxFun` packages.
 
 
 

src/abstract.jl

@@ -20,6 +20,7 @@ abstract type AbstractPolynomial{T,X} end
 # (and for LaurentPolynomial the offset)
 _convert(p::P, as) where {T,X,P <: AbstractPolynomial{T,X}} = ⟒(P)(as, X)
 
+
 """
     Polynomials.@register(name)
 
@@ -41,7 +42,10 @@ macro register(name)
     poly = esc(name)
     quote
         Base.convert(::Type{P}, p::P) where {P<:$poly} = p
-        Base.convert(P::Type{<:$poly}, p::$poly{T}) where {T} = constructorof(P){eltype(P), indeterminate(P,p)}(coeffs(p))
+        function Base.convert(P::Type{<:$poly}, p::$poly{T,X}) where {T,X}
+            isconstant(p) && return constructorof(P){eltype(P),indeterminate(P)}(constantterm(p))
+            constructorof(P){eltype(P), indeterminate(P,p)}(coeffs(p))
+        end
         Base.promote(p::P, q::Q) where {X, T, P <:$poly{T,X}, Q <: $poly{T,X}} = p,q
         Base.promote_rule(::Type{<:$poly{T,X}}, ::Type{<:$poly{S,X}}) where {T,S,X} =  $poly{promote_type(T, S),X}
         Base.promote_rule(::Type{<:$poly{T,X}}, ::Type{S}) where {T,S<:Number,X} =

src/common.jl

@@ -413,17 +413,31 @@ end
 Returns the complex conjugate of the polynomial
 """
 LinearAlgebra.conj(p::P) where {P <: AbstractPolynomial} = map(conj, p)
-LinearAlgebra.adjoint(p::P) where {P <: AbstractPolynomial} = map(adjoint, p) 
+LinearAlgebra.adjoint(p::P) where {P <: AbstractPolynomial} = map(adjoint, p)
+LinearAlgebra.adjoint(A::VecOrMat{<:AbstractPolynomial}) = adjoint.(permutedims(A)) # default has indeterminate indeterminate
 LinearAlgebra.transpose(p::AbstractPolynomial) = p
 LinearAlgebra.transpose!(p::AbstractPolynomial) = p
 
 #=
 Conversions =#
 Base.convert(::Type{P}, p::P) where {P <: AbstractPolynomial} = p
 Base.convert(P::Type{<:AbstractPolynomial}, x) = P(x)
+function Base.convert(::Type{S}, p::P) where {S <: Number,T, P<:Polynomials.AbstractPolynomial{T}}
+    Polynomials.isconstant(p) && return convert(S, Polynomials.constantterm(p))
+    throw(ArgumentError("Can't convert a nonconstant polynomial to type $S"))
+end
+
+# Methods to ensure that matrices of polynomials behave as desired
 Base.promote_rule(::Type{<:AbstractPolynomial{T}},
-    ::Type{<:AbstractPolynomial{S}},
-) where {T,S} = Polynomial{promote_type(T, S)}
+                  ::Type{<:AbstractPolynomial{S}},
+                  ) where {T,S} = Polynomial{promote_type(T, S)}
+Base.promote_rule(::Type{P},::Type{Q}) where {T,X, P<:AbstractPolynomial{T,X},
+                                              S,   Q<:AbstractPolynomial{S,X}} =
+                                                   Polynomial{promote_type(T, S),X}
+Base.promote_rule(::Type{P},::Type{Q}) where {T,X, P<:AbstractPolynomial{T,X},
+                                              S,Y, Q<:AbstractPolynomial{S,Y}} =
+                                                  assert_same_variable(X,Y)
+
 
 #=
 Inspection =#
@@ -453,6 +467,7 @@ function _eltype(P::Type{<:AbstractPolynomial}, p::AbstractPolynomial)
 end
 Base.iszero(p::AbstractPolynomial) = all(iszero, p)
 
+
 # See discussions in https://github.com/JuliaMath/Polynomials.jl/issues/258
 """
     all(pred, poly::AbstractPolynomial)
@@ -543,9 +558,6 @@ return `p(0)`, the constant term in the standard basis
 """
 constantterm(p::AbstractPolynomial{T}) where {T} = p(zero(T))
 
-
-
-
 """
     degree(::AbstractPolynomial)
 
@@ -720,8 +732,12 @@ function indeterminate(::Type{P}) where {P <: AbstractPolynomial}
     X == nothing ? :x : X
 end
 indeterminate(p::P) where {P <: AbstractPolynomial} = _indeterminate(P)
-function indeterminate(PP::Type{P}, p::AbstractPolynomial) where {P <: AbstractPolynomial}
-    X = _indeterminate(PP) == nothing ? indeterminate(p) :  _indeterminate(PP)
+function indeterminate(PP::Type{P}, p::AbstractPolynomial{T,Y}) where {P <: AbstractPolynomial, T,Y}
+    X = _indeterminate(PP)
+    X == nothing && return Y
+    assert_same_variable(X,Y)
+    return X
+    #X = _indeterminate(PP) == nothing ? indeterminate(p) :  _indeterminate(PP)
 end
 function indeterminate(PP::Type{P}, x::Symbol) where {P <: AbstractPolynomial}
     X = _indeterminate(PP) == nothing ? x :  _indeterminate(PP)
@@ -776,7 +792,7 @@ julia> p = 100 + 24x - 3x^2
 Polynomial(100 + 24*x - 3*x^2)
 
 julia> roots((x - 3) * (x + 2))
-2-element Array{Float64,1}:
+2-element Vector{Float64}:
  -2.0
   3.0
 

src/pade.jl

@@ -94,11 +94,14 @@ julia> using Polynomials, Polynomials.PolyCompat, SpecialFunctions
 
 
 
+
 julia> p = Polynomial(@.(1 // BigInt(gamma(1:17))));
 
 
+
 julia> pade = Pade(p, 8, 8);
 
+
 julia> pade(1.0) ≈ exp(1.0)
 true
 

src/polynomials/ChebyshevT.jl

@@ -35,6 +35,9 @@ julia> Polynomials.evalpoly(5.0, p, false) # bypasses the domain check done in p
 
 The latter shows how to evaluate a `ChebyshevT` polynomial outside of its domain, which is `[-1,1]`. (For newer versions of `Julia`, `evalpoly` is an exported function from Base with methods extended in this package, so the module qualification is unnecessary.
 
+!!! Note:
+    The Chebyshev polynomials are also implemented in `ApproxFun`, `ClassicalOrthogonalPolynomials.jl`, `FastTransforms.jl`, and `SpecialPolynomials.jl`.
+
 """
 struct ChebyshevT{T <: Number, X} <: AbstractPolynomial{T, X}
     coeffs::Vector{T}
@@ -102,7 +105,7 @@ julia> c(0)
 1.5
 
 julia> c.(-1:0.5:1)
-5-element Array{Float64,1}:
+5-element Vector{Float64}:
  2.0
  1.25
  1.5

src/polynomials/LaurentPolynomial.jl

@@ -135,12 +135,14 @@ function Base.convert(P::Type{<:Polynomial}, q::LaurentPolynomial)
     P([q[i] for i  in 0:n], indeterminate(q))
 end
 
-# save variable if specified in P
-Base.convert(::Type{P}, p::LaurentPolynomial) where {T, X, P<:LaurentPolynomial{T,X}} = P(p.coeffs, p.m[])
+# need to add p.m[], so abstract.jl method isn't sufficent
 function Base.convert(::Type{P}, p::LaurentPolynomial) where {P<:LaurentPolynomial}
-    T = eltype(P)
-    X = indeterminate(P,p)
-    ⟒(P){T, X}(convert(Vector{T},p.coeffs), p.m[])
+    S′ = _eltype(P)
+    Y′ = _indeterminate(P)
+    S = S′ == nothing ? eltype(p) : S′
+    Y = Y′ == nothing ? indeterminate(p) : Y′
+    isconstant(p) && return LaurentPolynomial{S,Y}(constantterm(p))
+    LaurentPolynomial{S,Y}(p.coeffs, p.m[])
 end
 
 function Base.convert(::Type{P}, q::StandardBasisPolynomial{S}) where {T, P <:LaurentPolynomial{T},S}