drop OffsetArrays dependency; clean up scalar operations for Immutabl…

…ePolynomial

Project.toml

@@ -7,20 +7,20 @@ version = "1.2.0"
 [deps]
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 
+
 [compat]
 Intervals = "0.5, 1.0, 1.3"
-OffsetArrays = "1"
 RecipesBase = "0.7, 0.8, 1"
 julia = "1"
 
 [extras]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["LinearAlgebra", "SparseArrays", "SpecialFunctions", "Test"]
+test = ["LinearAlgebra", "SparseArrays", "OffsetArrays", "SpecialFunctions", "Test"]

src/Polynomials.jl

@@ -3,7 +3,6 @@ module Polynomials
 # using GenericLinearAlgebra ## remove for now. cf: https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/pull/71#issuecomment-743928205
 using LinearAlgebra
 using Intervals
-using OffsetArrays
 
 include("abstract.jl")
 include("show.jl")

src/polynomials/ChebyshevT.jl

@@ -30,16 +30,15 @@ ChebyshevT(1.0⋅T_0(x))
 struct ChebyshevT{T <: Number, X} <: AbstractPolynomial{T, X}
  coeffs::Vector{T}
  function ChebyshevT{T, X}(coeffs::AbstractVector{S}) where {T <: Number,X, S}
- length(coeffs) == 0 && return new{T,X}(zeros(T, 1))
+
  if Base.has_offset_axes(coeffs)
- # throw(ArgumentError("The `ChebyshevT` constructor does not accept `OffsetArrays`. Try `LaurentPolynomial`."))
-
  @warn "ignoring the axis offset of the coefficient vector"
- coeffs = OffsetArrays.no_offset_view(coeffs) 
  end
- last_nz = findlast(!iszero, coeffs)
- last = max(1, last_nz === nothing ? 0 : last_nz)
- return new{T,X}(convert(Vector{T}, coeffs[1:last]))
+
+ N = findlast(!iszero, coeffs)
+ isnothing(N) && return new{T,X}(zeros(T,1))
+ cs = T[coeffs[i] for i ∈ firstindex(coeffs):N]
+ new{T,X}(cs)
  end
 end
 

src/polynomials/ImmutablePolynomial.jl

@@ -63,20 +63,19 @@ function ImmutablePolynomial{T,X}(coeffs::AbstractVector{S}) where {T,X,S}
  R = promote_type(T,S)
 
  if Base.has_offset_axes(coeffs)
- # throw(ArgumentError("The `ImmutablePolynomial` constructor does not accept `OffsetArrays`. Try `LaurentPolynomial`."))
  @warn "ignoring the axis offset of the coefficient vector"
- coeffs = OffsetArrays.no_offset_view(coeffs) 
  end
- #c = OffsetArrays.no_offset_view(coeffs)
  N = findlast(!iszero, coeffs)
- N == nothing && return zero(ImmutablePolynomial{R,X})
- ImmutablePolynomial{T, X, N}(NTuple{N,T}(cᵢ for cᵢ ∈ coeffs))
+ isnothing(N) && return ImmutablePolynomial{R,X,0}(())
+ N′ = N + 1 - firstindex(coeffs)
+ cs = NTuple{N′,T}(coeffs[i] for i ∈ firstindex(coeffs):N)
+ ImmutablePolynomial{T, X, N′}(cs)
 end
 
 ## -- Tuple arguments
 function ImmutablePolynomial{T,X}(coeffs::Tuple) where {T,X}
  N = findlast(!iszero, coeffs)
- N == nothing && return zero(ImmutablePolynomial{T,X})
+ isnothing(N) && return zero(ImmutablePolynomial{T,X})
  ImmutablePolynomial{T,X,N}(NTuple{N,T}(coeffs[i] for i in 1:N))
 end
 
@@ -276,14 +275,15 @@ end
 
 function Base.:*(p::ImmutablePolynomial{T,X,N}, c::S) where {T, X,N, S <: Number}
  R = eltype(one(T)*one(S))
- iszero(c) && return zero(ImmutablePolynomial{R,X})
+ iszero(p[end]*c) && return ImmutablePolynomial{R,X}(p.coeffs .* c)
  ImmutablePolynomial{R,X,N}(p.coeffs .* c)
 end
 
 function Base.:/(p::ImmutablePolynomial{T,X,N}, c::S) where {T,X,N,S <: Number}
  R = eltype(one(T)/one(S))
- isinf(c) && return zero(ImmutablePolynomial{R,X})
- ImmutablePolynomial{R,X,N}(p.coeffs ./ c)
+ cs = p.coeffs ./ c
+ iszero(cs[end]) && return ImmutablePolynomial{R,X}(cs)
+ ImmutablePolynomial{R,X,N}(cs)
 end
 
 Base.:-(p::ImmutablePolynomial{T,X,N}) where {T,X,N} = ImmutablePolynomial{T,X,N}(.-p.coeffs)

src/polynomials/Polynomial.jl

@@ -34,14 +34,16 @@ struct Polynomial{T <: Number, X} <: StandardBasisPolynomial{T, X}
  coeffs::Vector{T}
  function Polynomial{T, X}(coeffs::AbstractVector{S}) where {T <: Number, X, S}
  if Base.has_offset_axes(coeffs)
- #throw(ArgumentError("The `Polynomial` constructor does not accept `OffsetArrays`. Try `LaurentPolynomial`."))
  @warn "ignoring the axis offset of the coefficient vector"
- coeffs = OffsetArrays.no_offset_view(coeffs) 
  end
- length(coeffs) == 0 && return new{T,X}(zeros(T, 1))
- last_nz = findlast(!iszero, coeffs)
- last = max(1, last_nz === nothing ? 0 : last_nz)
- return new{T, X}(convert(Vector{T}, coeffs[1:last]))
+ N = findlast(!iszero, coeffs)
+ isnothing(N) && return new{T,X}(zeros(T,1))
+ cs = T[coeffs[i] for i ∈ firstindex(coeffs):N]
+ new{T,X}(cs)
+ end
+ # non-copying alternative assuming !iszero(coeffs[end])
+ function Polynomial{T, X}(checked::Val{false}, coeffs::AbstractVector{T}) where {T <: Number, X}
+ new{T, X}(coeffs)
  end
 end
 
@@ -73,68 +75,64 @@ julia> p.(0:3)
 
 # scalar _,* faster than standard-basis/common versions
 function Base.:+(p::P, c::S) where {T, X, P <: Polynomial{T, X}, S<:Number}
- R = promote_type(T, S)
- as = convert(Vector{R}, copy(coeffs(p)))
- as[1] += c
- return Polynomial{R, X}(as)
+ R = promote_type(T, S)
+ Q = Polynomial{R,X}
+ as = convert(Vector{R}, copy(coeffs(p)))
+ as[1] += c
+ iszero(as[end]) ? Q(as) : Q(Val(false), as)
 end
 
 function Base.:*(p::P, c::S) where {T, X, P <: Polynomial{T,X} , S <: Number}
- as = [aᵢ * c for aᵢ ∈ coeffs(p)]
- Polynomial{promote_type(T,S),X}(as)
+ R = promote_type(T,S)
+ Q = Polynomial{R, X}
+ as = R[aᵢ * c for aᵢ ∈ coeffs(p)]
+ iszero(as[end]) ? Q(as) : Q(Val(false), as)
 end
 
-
-
 function Base.:+(p1::Polynomial{T}, p2::Polynomial{S}) where {T, S}
+ isconstant(p1) && return p2 + p1[0]
+ isconstant(p2) && return p1 + p2[0]
+ X, Y = indeterminate(p1), indeterminate(p2)
+ X != Y && error("Polynomials must have same variable")
  n1, n2 = length(p1), length(p2)
- if n1 > 1 && n2 > 1
- indeterminate(p1) != indeterminate(p2) && error("Polynomials must have same variable")
- end
  R = promote_type(T,S)
+
  c = zeros(R, max(n1, n2))
- if n1 > 1 && n2 > 1
- if n1 >= n2
- c .= p1.coeffs
- for i = eachindex(p2.coeffs)
+ if n1 >= n2
+ c .= p1.coeffs
+ for i in eachindex(p2.coeffs)
  c[i] += p2.coeffs[i]
- end
- else
- c .= p2.coeffs
- for i = eachindex(p1.coeffs)
- c[i] += p1.coeffs[i]
- end
  end
- return Polynomial{R, indeterminate(p1)}(c)
- elseif n1 <= 1
- c .= p2.coeffs
- c[1] += p1[0]
- return Polynomial{R, indeterminate(p2)}(c)
- else 
- c .= p1.coeffs
- c[1] += p2[0]
- return Polynomial{R, indeterminate(p1)}(c)
+ else
+ c .= p2.coeffs
+ for i in eachindex(p1.coeffs)
+ c[i] += p1.coeffs[i]
+ end
  end
+
+ Q = Polynomial{R,X}
+ return iszero(c[end]) ? Q(c) : Q(Val(false), c)
+
 end
 
 
 function Base.:*(p1::Polynomial{T}, p2::Polynomial{S}) where {T,S}
 
  n, m = length(p1)-1, length(p2)-1 # not degree, so pNULL works
+ X, Y = indeterminate(p1), indeterminate(p2)
+ R = promote_type(T, S)
  if n > 0 && m > 0
- indeterminate(p1) != indeterminate(p2) && error("Polynomials must have same variable")
- R = promote_type(T, S)
+ X != Y && error("Polynomials must have same variable")
  c = zeros(R, m + n + 1)
  for i in 0:n, j in 0:m
  @inbounds c[i + j + 1] += p1[i] * p2[j]
  end
- return Polynomial{R, indeterminate(p1)}(c)
+ Q = Polynomial{R,X}
+ return iszero(c[end]) ? Q(c) : Q(Val(false), c)
  elseif n <= 0
- cs = p2.coeffs * p1[0]
- return Polynomial{eltype(cs), indeterminate(p2)}(cs)
+ return Polynomial{R, Y}(p2.coeffs * p1[0])
  else
- cs = p1.coeffs * p2[0]
- return Polynomial{eltype(cs), indeterminate(p1)}(cs)
+ return Polynomial{R, X}(p1.coeffs * p2[0])
  end
 
 end

src/polynomials/SparsePolynomial.jl

@@ -48,6 +48,9 @@ struct SparsePolynomial{T <: Number, X} <: StandardBasisPolynomial{T, X}
  end
  new{T, X}(c)
  end
+ function SparsePolynomial{T,X}(checked::Val{false}, coeffs::AbstractDict{Int, T}) where {T <: Number, X}
+ new{T,X}(convert(Dict{Int,S}, coeffs))
+ end
 end
 
 @register SparsePolynomial
@@ -61,17 +64,13 @@ function SparsePolynomial(coeffs::AbstractDict{Int, T}, var::SymbolLike=:x) wher
 end
 
 function SparsePolynomial{T,X}(coeffs::AbstractVector{S}) where {T <: Number, X, S}
- firstindex(coeffs) >= 0 || throw(ArgumentError("Use the `LaurentPolynomial` type for arrays with a negative first index"))
 
  if Base.has_offset_axes(coeffs)
- # throw(ArgumentError("The `SparsePolynomial` constructor does not accept `OffsetArrays`. Try `LaurentPolynomial`."))
-
  @warn "ignoring the axis offset of the coefficient vector"
- coeffs = OffsetArrays.no_offset_view(coeffs) 
  end
 
- c = OffsetArrays.no_offset_view(coeffs) # ensure 1-based indexing
- p = Dict{Int,T}(i - 1 => v for (i,v) in pairs(c))
+ offset = firstindex(coeffs)
+ p = Dict{Int,T}(k - offset => v for (k,v) ∈ pairs(coeffs))
  return SparsePolynomial{T,X}(p)
 end
 

test/Project.toml

@@ -0,0 +1,7 @@
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"