Merge pull request #264 from jverzani/issue_258a

integrate in map, any, all, as possible

Author: jverzani

src/abstract.jl

@@ -17,6 +17,10 @@ abstract type AbstractPolynomial{T} end
 # works for most cases
 ⟒(P::Type{<:AbstractPolynomial}) = constructorof(P)
 
+# convert `as` into polynomial of type P based on instance, inheriting variable
+# (and for LaurentPolynomial the offset)
+_convert(p::P, as) where {P <: AbstractPolynomial} = ⟒(P)(as, p.var)
+
 """
     Polynomials.@register(name)
 

src/common.jl

@@ -120,7 +120,7 @@ Returns the roots of the given polynomial. This is calculated via the eigenvalue
 
 !!! note
 
-        The [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) package provides an alternative that is a bit faster and a bit more accurate; the [FastPolynomialRoots](https://github.com/andreasnoack/FastPolynomialRoots.jl) provides an interface to FORTRAN code implementing an algorithm that can handle very large polynomials (it is  `O(n^2)` not `O(n^3)`. the [AMRVW.jl](https://github.com/jverzani/AMRVW.jl) package implements the algorithm in Julia, allowing the use of other  number types.
+        The [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) package provides an alternative that is a bit faster and a bit more accurate; the [FastPolynomialRoots](https://github.com/andreasnoack/FastPolynomialRoots.jl) provides an interface to FORTRAN code implementing an algorithm that can handle very large polynomials (it is  `O(n^2)` not `O(n^3)`. The [AMRVW.jl](https://github.com/jverzani/AMRVW.jl) package implements the algorithm in Julia, allowing the use of other  number types.
 
 """
 function roots(q::AbstractPolynomial{T}; kwargs...) where {T <: Number}
@@ -210,8 +210,8 @@ In-place version of [`chop`](@ref)
 function chop!(p::AbstractPolynomial{T};
     rtol::Real = Base.rtoldefault(real(T)),
                atol::Real = 0,) where {T}
-    isempty(coeffs(p)) && return p
-    tol = norm(coeffs(p)) * rtol + atol
+    isempty(values(p)) && return p
+    tol = norm(p) * rtol + atol
     for i = lastindex(p):-1:0
         val = p[i]
         if abs(val) > tol #!isapprox(val, zero(T); rtol = rtol, atol = atol)
@@ -253,15 +253,19 @@ Linear Algebra =#
 
 Calculates the p-norm of the polynomial's coefficients
 """
-LinearAlgebra.norm(q::AbstractPolynomial, p::Real = 2) = norm(coeffs(q), p)
+function LinearAlgebra.norm(q::AbstractPolynomial, p::Real = 2)
+    vs = values(q)
+    isempty(vs) && return zero(q[0])
+    return norm(vs, p)
+end
 
 """
     conj(::AbstractPolynomial)
 
 Returns the complex conjugate of the polynomial
 """
-LinearAlgebra.conj(p::P) where {P <: AbstractPolynomial} = ⟒(P)(conj(coeffs(p)), p.var)
-LinearAlgebra.adjoint(p::P) where {P <: AbstractPolynomial} = ⟒(P)(adjoint.(coeffs(p)), p.var)
+LinearAlgebra.conj(p::P) where {P <: AbstractPolynomial} = map(conj, p)
+LinearAlgebra.adjoint(p::P) where {P <: AbstractPolynomial} = map(adjoint, p) 
 LinearAlgebra.transpose(p::AbstractPolynomial) = p
 LinearAlgebra.transpose!(p::AbstractPolynomial) = p
 
@@ -293,12 +297,7 @@ Base.eltype(p::AbstractPolynomial{T}) where {T} = T
 Base.eltype(::Type{<:AbstractPolynomial}) = Float64
 Base.eltype(::Type{<:AbstractPolynomial{T}}) where {T} = T
 #Base.eltype(::Type{P}) where {P <: AbstractPolynomial} = P # changed  in v1.1.0
-function Base.iszero(p::AbstractPolynomial)
-    if length(p) == 0
-        return true
-    end
-    return all(iszero.(coeffs(p))) && p[0] == 0
-end
+Base.iszero(p::AbstractPolynomial) = all(iszero, p)
 
 # See discussions in https://github.com/JuliaMath/Polynomials.jl/issues/258
 """
@@ -308,21 +307,25 @@ Test whether all coefficients of an `AbstractPolynomial` satisfy predicate `pred
 
 You can implement `isreal`, etc., to a `Polynomial` by using `all`.
 """
-Base.all(pred, poly::AbstractPolynomial) = all(pred, poly[:])
+Base.all(pred, p::AbstractPolynomial) = all(pred, values(p))
 """
     any(pred, poly::AbstractPolynomial)
 
 Test whether any coefficient of an `AbstractPolynomial` satisfies predicate `pred`.
 """
-Base.any(pred, poly::AbstractPolynomial) = any(pred, poly[:])
+Base.any(pred, p::AbstractPolynomial) = any(pred, values(p))
+
+
+
+
 """
-    map(fn, p::AbstractPolynomial)
+    map(fn, p::AbstractPolynomial, args...)
 
 Transform coefficients of `p` by applying a function (or other callables) `fn` to each of them.
 
 You can implement `real`, etc., to a `Polynomial` by using `map`.
 """
-Base.map(fn, p::P) where {P<:AbstractPolynomial} = ⟒(P)(map(fn, coeffs(p)), p.var)
+Base.map(fn, p::P, args...)  where {P<:AbstractPolynomial} = _convert(p, map(fn, coeffs(p), args...))
 
 """
     isreal(p::AbstractPolynomial)
@@ -371,7 +374,7 @@ coeffs(p::AbstractPolynomial) = p.coeffs
 Return the degree of the polynomial, i.e. the highest exponent in the polynomial that
 has a nonzero coefficient. The degree of the zero polynomial is defined to be -1.
 """
-degree(p::AbstractPolynomial) = iszero(p) ? -1 : length(p) - 1
+degree(p::AbstractPolynomial) = iszero(p) ? -1 : lastindex(p) 
 
 
 """
@@ -384,7 +387,7 @@ isconstant(p::AbstractPolynomial) = degree(p) <= 0
 
 
 
-hasnan(p::AbstractPolynomial) = any(isnan.(coeffs(p)))
+hasnan(p::AbstractPolynomial) = any(isnan, p)
 
 """
     domain(::Type{<:AbstractPolynomial})
@@ -427,6 +430,8 @@ Base.firstindex(p::AbstractPolynomial) = 0
 Base.lastindex(p::AbstractPolynomial) = length(p) - 1
 Base.eachindex(p::AbstractPolynomial) = 0:length(p) - 1
 Base.broadcastable(p::AbstractPolynomial) = Ref(p)
+# like coeffs, though possibly only non-zero values (e.g. SparsePolynomial)
+Base.values(p::AbstractPolynomial) = coeffs(p) 
 
 # iteration
 # iteration occurs over the basis polynomials
@@ -472,7 +477,7 @@ Base.setindex!(p::AbstractPolynomial, values, ::Colon) =
 
 #=
 identity =#
-Base.copy(p::P) where {P <: AbstractPolynomial} = P(copy(coeffs(p)), p.var)
+Base.copy(p::P) where {P <: AbstractPolynomial} = _convert(p, copy(coeffs(p)))
 Base.hash(p::AbstractPolynomial, h::UInt) = hash(p.var, hash(coeffs(p), h))
 
 #=
@@ -540,20 +545,18 @@ basis(p::P, k::Int, _var::SymbolLike=:x; var=_var) where {P<:AbstractPolynomial}
 
 #=
 arithmetic =#
-Base.:-(p::P) where {P <: AbstractPolynomial} = P(-coeffs(p), p.var)
+Base.:-(p::P) where {P <: AbstractPolynomial} = _convert(p, -coeffs(p)) 
 Base.:+(c::Number, p::AbstractPolynomial) = +(p, c)
 Base.:-(p::AbstractPolynomial, c::Number) = +(p, -c)
 Base.:-(c::Number, p::AbstractPolynomial) = +(-p, c)
 Base.:*(c::Number, p::AbstractPolynomial) = *(p, c)
 
 function Base.:*(p::P, c::S) where {P <: AbstractPolynomial,S}
-    T = promote_type(P, S)
-    return T(coeffs(p) .* c, p.var)
+    _convert(p, coeffs(p) .* c)
 end
 
 function Base.:/(p::P, c::S) where {T,P <: AbstractPolynomial{T},S}
-    R = promote_type(P, eltype(one(T) / one(S)))
-    return R(coeffs(p) ./ c, p.var)
+    _convert(p, coeffs(p) ./ c)
 end
 
 Base.:-(p1::AbstractPolynomial, p2::AbstractPolynomial) = +(p1, -p2)
@@ -644,7 +647,6 @@ function Base.isapprox(p1::AbstractPolynomial{T},
 
     p1, p2 = promote(p1, p2)
     check_same_variable(p1, p2)  || error("p1 and p2 must have same var")
-
     # copy over from abstractarray.jl
     Δ  = norm(p1-p2)
     if isfinite(Δ)
@@ -661,7 +663,7 @@ function Base.isapprox(p1::P,
                        n::S;
                        rtol::Real = (Base.rtoldefault(T, S, 0)),
                        atol::Real = 0,) where {T,S, P<:AbstractPolynomial{T}}
-    return isapprox(p1, ⟒(P){T}(n,p1.var))
+    return isapprox(p1, _convert(p1, [n])) 
 end
 
 Base.isapprox(n::S,

src/polynomials/ImmutablePolynomial.jl

@@ -92,6 +92,7 @@ function ImmutablePolynomial(coeffs::NTuple{M,T}, var::SymbolLike=:x) where {M,
     ImmutablePolynomial{T,N}(cs, var)
 end
 
+
 # Convenience; pass tuple to Polynomial
 # Not documented, not sure this is a good idea as P(...)::P is not true...
 Polynomial(coeffs::NTuple{N,T}, var::SymbolLike = :x) where{N,T} =
@@ -110,7 +111,6 @@ Base.collect(p::P) where {P <: ImmutablePolynomial} = [pᵢ for pᵢ ∈ p]
 Base.copy(p::P) where {P <: ImmutablePolynomial} = P(coeffs(p), p.var)
 
 # catch q == 0 case
-LinearAlgebra.norm(q::ImmutablePolynomial{T}, p::Real = 2) where {T} = degree(q) == -1 ? zero(T) : norm(coeffs(q), p)
 
 #  zero, one, variable
 function Base.zero(P::Type{<:ImmutablePolynomial},var::SymbolLike=:x)

src/polynomials/LaurentPolynomial.jl

@@ -93,12 +93,12 @@ end
 @register LaurentPolynomial
 
 ## constructors
-function LaurentPolynomial{T}(coeffs::AbstractVector{S}, m::Int, var::Symbol=:x) where {
+function LaurentPolynomial{T}(coeffs::AbstractVector{S}, m::Int, var::SymbolLike=:x) where {
     T <: Number, S <: Number}
     LaurentPolynomial{T}(T.(coeffs), m, var)
 end
 
-function LaurentPolynomial{T}(coeffs::AbstractVector{S}, var::Symbol=:x) where {
+function LaurentPolynomial{T}(coeffs::AbstractVector{S}, var::SymbolLike=:x) where {
     T <: Number, S <: Number}
     LaurentPolynomial{T}(T.(coeffs), 0, var)
 end
@@ -114,7 +114,8 @@ end
 
 
 # Add interface for OffsetArray
-function  LaurentPolynomial{T}(coeffs::OffsetArray{S, 1, Array{S,1}}, var::SymbolLike=:x) where {T, S}
+function  LaurentPolynomial{T}(coeffs::OffsetArray{S, 1, Array{S,1}}, var::SymbolLike=:x) where {
+    T<:Number, S<:Number}
     m,n = axes(coeffs, 1)
     LaurentPolynomial{T}(T.(coeffs.parent), m:n, Symbol(var))
 end
@@ -177,21 +178,19 @@ function Base.range(p::LaurentPolynomial)
 end
 
 function Base.inv(p::LaurentPolynomial)
-    m,n =  degreerange(p)
+    m,n =  (extrema∘degreerange)(p)
     m != n && throw(ArgumentError("Only monomials can be inverted"))
     LaurentPolynomial([1/p for p in p.coeffs], -m, p.var)
 end
 
 ##
 ## changes to common.jl mostly as the range in the type is different
 ##
-Base.copy(p::P) where {P <: LaurentPolynomial} = P(copy(coeffs(p)), degreerange(p), p.var)
 Base.:(==)(p1::LaurentPolynomial, p2::LaurentPolynomial) =
     check_same_variable(p1, p2) && (degreerange(p1) == degreerange(p2)) && (coeffs(p1) == coeffs(p2))
 Base.hash(p::LaurentPolynomial, h::UInt) = hash(p.var, hash(degreerange(p), hash(coeffs(p), h)))
 
-degree(p::LaurentPolynomial) = p.n[]
-isconstant(p::LaurentPolynomial) = degreerange(p) == 0:0
+isconstant(p::LaurentPolynomial) = iszero(lastindex(p)) && iszero(firstindex(p)) 
 basis(P::Type{<:LaurentPolynomial{T}}, n::Int, var::SymbolLike=:x) where{T} = LaurentPolynomial(ones(T,1), n:n, var)
 basis(P::Type{LaurentPolynomial}, n::Int, var::SymbolLike=:x) = LaurentPolynomial(ones(Float64, 1), n:n, var)
 
@@ -234,6 +233,8 @@ Base.lastindex(p::LaurentPolynomial) = p.n[]
 Base.eachindex(p::LaurentPolynomial) = degreerange(p)
 degreerange(p::LaurentPolynomial) = firstindex(p):lastindex(p)
 
+_convert(p::P, as) where {P <: LaurentPolynomial} = ⟒(P)(as, firstindex(p), p.var)
+
 ## chop/truncation
 # trim  from *both* ends
 function chop!(p::P;
@@ -326,11 +327,7 @@ julia> conj(p)(z) ≈ (conj ∘ p ∘ conj)(z)
 true
 ```
 """
-function LinearAlgebra.conj(p::P) where {P <: LaurentPolynomial}
-    ps = coeffs(p)
-    m = firstindex(p)
-    ⟒(P)(conj(ps), m, p.var)
-end
+LinearAlgebra.conj(p::P) where {P <: LaurentPolynomial} = map(conj, p)
 
 
 """
@@ -443,23 +440,8 @@ end
 
 
 # scalar operattoinis
-Base.:-(p::P) where {P <: LaurentPolynomial} = P(-coeffs(p), firstindex(p), p.var)
-
-function Base.:+(p::LaurentPolynomial{T}, c::S) where {T, S <: Number}
-    q = LaurentPolynomial([c], 0, p.var)
-    p + q
-end
-
-function Base.:*(p::P, c::S) where {T,P <: LaurentPolynomial,  S <: Number}
-    as = c * copy(coeffs(p))
-    return ⟒(P)(as, firstindex(p), p.var)
-end
-
-
-function Base.:/(p::P, c::S) where {T,P <: LaurentPolynomial{T},S <: Number}
-    R = promote_type(P, eltype(one(T) / one(S)))
-    return R(coeffs(p) ./ c, firstindex(p), p.var)
-end
+# standard-basis defn. assumes basis 1, x, x², ...
+Base.:+(p::LaurentPolynomial{T}, c::S) where {T, S <: Number} = sum(promote(p,c))
 
 ##
 ## Poly + and  *

src/polynomials/SparsePolynomial.jl

@@ -182,6 +182,13 @@ function (p::SparsePolynomial{T})(x::S) where {T,S}
 
 end
 
+#  map: over values -- not keys
+function Base.map(fn, p::P, args...) where {P <: SparsePolynomial}
+    ks, vs = keys(p.coeffs), values(p.coeffs)
+    vs′ = map(fn, vs, args...)
+    _convert(p, Dict(Pair.(ks, vs′)))
+end
+
 
 
 function Base.:+(p1::SparsePolynomial{T}, p2::SparsePolynomial{S}) where {T, S}

src/polynomials/standard-basis.jl

@@ -24,6 +24,8 @@ isconstant(p::StandardBasisPolynomial) = degree(p) <= 0
 
 Base.convert(P::Type{<:StandardBasisPolynomial}, q::StandardBasisPolynomial) = isa(q, P) ? q : P([q[i] for i in 0:degree(q)], q.var)
 
+Base.values(p::StandardBasisPolynomial) = values(p.coeffs)
+
 variable(::Type{P}, var::SymbolLike = :x) where {P <: StandardBasisPolynomial} = P([0, 1], var)
 
 function fromroots(P::Type{<:StandardBasisPolynomial}, r::AbstractVector{T}; var::SymbolLike = :x) where {T <: Number}
@@ -45,7 +47,7 @@ function Base.:+(p::P, c::S) where {T, P <: StandardBasisPolynomial{T}, S<:Numbe
     R = promote_type(T,S)
     as = R[c  for c in coeffs(p)]
     as[1] += c
-    ⟒(P)(as, p.var)
+    _convert(p, as)
 end
 
 
@@ -66,7 +68,7 @@ function derivative(p::P, order::Integer = 1) where {T, P <: StandardBasisPolyno
     @inbounds for i in order:n - 1
         a2[i - order + 1] = reduce(*, (i - order + 1):i, init = p[i])
     end
-    return ⟒(P)(a2, p.var)
+    return _convert(p, a2)
 end
 
 
@@ -82,7 +84,7 @@ function integrate(p::P, k::S) where {T, P <: StandardBasisPolynomial{T}, S<:Num
     @inbounds for i in 1:n
         a2[i + 1] = p[i - 1] / i
     end
-    return ⟒(P)(a2, p.var)
+    return _convert(p, a2)
 end
 
 
@@ -119,7 +121,7 @@ function Base.divrem(num::P, den::Q) where {T, P <: StandardBasisPolynomial{T},
     end
     resize!(r_coeff, min(length(r_coeff), m))
 
-    return ⟒(P)(q_coeff, var), ⟒(P)(r_coeff, var)
+    return _convert(num, q_coeff), _convert(num, r_coeff)
 
 end
 
@@ -191,7 +193,7 @@ function  gcd_noda_sasaki(p::P, q::Q;
     a, b = coeffs(p), coeffs(q)
     as =  _gcd_noda_sasaki(a,b, atol=atol,  rtol=rtol)
 
-    ⟒(P)(as, p.var)
+    _convert(p, as)
 end
 
 function _gcd_noda_sasaki(a::Vector{T}, b::Vector{S};
@@ -304,7 +306,7 @@ end
     compensated_horner(p::P, x)
     compensated_horner(ps, x)
 
-Evaluate `p(x)` using a compensation scheme of S. Graillat, Ph. Langlois, N. Louve [Compensated Horner Scheme](https://cadxfem.org/cao/Compensation-horner.pdf). Either a `Polynomial` `p` or it coefficients may be passed in.
+Evaluate `p(x)` using a compensation scheme of S. Graillat, Ph. Langlois, N. Louve [Compensated Horner Scheme](https://cadxfem.org/cao/Compensation-horner.pdf). Either a `Polynomial` `p` or its coefficients may be passed in.
 
 The Horner scheme has relative error given by
 
@@ -379,7 +381,7 @@ end
 # rule of thumb: p̂ a compute value
 # |p(x) - p̃(x)|/|p(x)| ≤ α(n)⋅u ⋅ cond(p,x), where u = finite precision of compuation (2^-p)
 function LinearAlgebra.cond(p::P, x) where {P <: Polynomials.StandardBasisPolynomial}
-    p̃ = Polynomials.:⟒(P)(abs.(coeffs(p)), p.var)
+    p̃ = map(abs, p)
     p̃(abs(x))/ abs(p(x))
 end
 

test/StandardBasis.jl

@@ -155,6 +155,11 @@ end
         @test p+q == P([1+im, 2, 3])
         @test p*q == P(im*[1,2,3])
     end
+
+    # LaurentPolynomial has an inverse for monomials
+    x = variable(LaurentPolynomial)
+    @test Polynomials.isconstant(x * inv(x))
+    @test_throws ArgumentError inv(x + x^2)
 end
 
 @testset "Divrem" begin