Merge pull request #72 from JuliaAlgebra/fix-0.7

Towards v0.7 compatibility

Author: rdeits

REQUIRE

@@ -1 +1,2 @@
 julia 0.6
+Compat 0.49

src/MultivariatePolynomials.jl

@@ -6,7 +6,20 @@ module MultivariatePolynomials
 
 import Base: *, +, -, /, ^, ==,
     promote_rule, convert, show, isless, size, getindex,
-    one, zero, transpose, isapprox, @pure, dot, copy
+    one, zero, isapprox, @pure, copy
+using Compat
+import Compat.LinearAlgebra: dot, norm
+
+# The ' operator lowers to `transpose()` in v0.6 and to
+# `adjoint()` in v0.7+. 
+# TOOD: remove this switch when dropping v0.6 support
+@static if VERSION <= v"0.7.0-DEV.3351"
+    import Base: transpose
+    const adjoint_operator = transpose
+else
+    import Compat.LinearAlgebra: adjoint
+    const adjoint_operator = adjoint
+end
 
 export AbstractPolynomialLike, AbstractTermLike, AbstractMonomialLike
 """

src/comparison.jl

@@ -7,12 +7,12 @@ Base.iszero(t::AbstractPolynomial) = iszero(nterms(t))
 
 # See https://github.com/blegat/MultivariatePolynomials.jl/issues/22
 # avoids the call to be transfered to eqconstant
-(==)(α::Void, x::APL) = false
-(==)(x::APL, α::Void) = false
+(==)(α::Nothing, x::APL) = false
+(==)(x::APL, α::Nothing) = false
 (==)(α::Dict, x::APL) = false
 (==)(x::APL, α::Dict) = false
-(==)(α::Void, x::RationalPoly) = false
-(==)(x::RationalPoly, α::Void) = false
+(==)(α::Nothing, x::RationalPoly) = false
+(==)(x::RationalPoly, α::Nothing) = false
 (==)(α::Dict, x::RationalPoly) = false
 (==)(x::RationalPoly, α::Dict) = false
 
@@ -100,7 +100,7 @@ end
 
 # α could be a JuMP affine expression
 isapproxzero(α; ztol::Real=0.) = false
-function isapproxzero(α::Number; ztol::Real=Base.rtoldefault(α, α))
+function isapproxzero(α::Number; ztol::Real=Base.rtoldefault(α, α, 0))
     abs(α) <= ztol
 end
 

src/conversion.jl

@@ -1,8 +1,8 @@
 export variable
 
-function convertconstant end
+convertconstant(::Type{P}, α) where P<:APL = P(α)
 convert(::Type{P}, α) where P<:APL = convertconstant(P, α)
-convert(::Type{P}, p::APL) where P<:AbstractPolynomial = convert(P, polynomial(p))
+convert(::Type{P}, p::APL) where P<:AbstractPolynomial = P(polynomial(p))
 
 Base.convert(::Type{Any}, p::APL) = p
 # Conversion polynomial -> scalar
@@ -11,12 +11,16 @@ function Base.convert(::Type{S}, p::APL) where {S}
     for t in terms(p)
         if !isconstant(t)
             # The polynomial is not constant
-            throw(InexactError())
+            throw(InexactError(:convert, S, p))
         end
         s += S(coefficient(t))
     end
     s
 end
+
+# Fix ambiguity caused by above conversions
+Base.convert(::Type{P}, p::APL) where P<:APL = P(p)
+
 Base.convert(::Type{PT}, p::PT) where {PT<:APL} = p
 function Base.convert(::Type{MT}, t::AbstractTerm) where {MT<:AbstractMonomial}
     if isone(coefficient(t))

src/differentiation.jl

@@ -2,15 +2,20 @@
 export differentiate
 
 """
-    differentiate(p::AbstractPolynomialLike, v::AbstractVariable, deg::Int=1)
+    differentiate(p::AbstractPolynomialLike, v::AbstractVariable, deg::Union{Int, Val}=1)
 
 Differentiate `deg` times the polynomial `p` by the variable `v`.
 
-    differentiate(p::AbstractPolynomialLike, vs, deg::Int=1)
 
-Differentiate `deg` times the polynomial `p` by the variables of the vector or tuple of variable `vs` and return an array of dimension `deg`.
+    differentiate(p::AbstractPolynomialLike, vs, deg::Union{Int, Val}=1)
 
-    differentiate(p::AbstractArray{<:AbstractPolynomialLike, N}, vs, deg::Int=1) where N
+Differentiate `deg` times the polynomial `p` by the variables of the vector or
+tuple of variable `vs` and return an array of dimension `deg`. It is recommended
+to pass `deg` as a `Val` instance when the degree is known at compile time, e.g.
+`differentiate(p, v, Val{2}())` instead of `differentiate(p, x, 2)`, as this
+will help the compiler infer the return type.
+
+    differentiate(p::AbstractArray{<:AbstractPolynomialLike, N}, vs, deg::Union{Int, Val}=1) where N
 
 Differentiate the polynomials in `p` by the variables of the vector or tuple of variable `vs` and return an array of dimension `N+deg`.
 
@@ -19,63 +24,84 @@ Differentiate the polynomials in `p` by the variables of the vector or tuple of
 ```julia
 p = 3x^2*y + x + 2y + 1
 differentiate(p, x) # should return 6xy + 1
+differentiate(p, x, Val{1}()) # equivalent to the above
 differentiate(p, (x, y)) # should return [6xy+1, 3x^2+1]
 ```
 """
 function differentiate end
 
 # Fallback for everything else
-_diff_promote_op(::Type{T}, ::Type{<:AbstractVariable}) where T = T
 differentiate(α::T, v::AbstractVariable) where T = zero(T)
-
-_diff_promote_op(::Type{<:AbstractVariable}, ::Type{<:AbstractVariable}) = Int
 differentiate(v1::AbstractVariable, v2::AbstractVariable) = v1 == v2 ? 1 : 0
-
-_diff_promote_op(::Type{TT}, ::Type{<:AbstractVariable}) where {T, TT<:AbstractTermLike{T}} = changecoefficienttype(TT, Base.promote_op(*, T, Int))
 differentiate(t::AbstractTermLike, v::AbstractVariable) = coefficient(t) * differentiate(monomial(t), v)
-
-_diff_promote_op(::Type{PT}, ::Type{<:AbstractVariable}) where {T, PT<:APL{T}} = polynomialtype(PT, Base.promote_op(*, T, Int))
 # The polynomial function will take care of removing the zeros
 differentiate(p::APL, v::AbstractVariable) = polynomial(differentiate.(terms(p), v), SortedState())
-
 differentiate(p::RationalPoly, v::AbstractVariable) = (differentiate(p.num, v) * p.den - p.num * differentiate(p.den, v)) / p.den^2
 
 const ARPL = Union{APL, RationalPoly}
 
-_vec_diff_promote_op(::Type{PT}, ::AbstractVector{VT}) where {PT, VT} = _diff_promote_op(PT, VT)
-_vec_diff_promote_op(::Type{PT}, ::NTuple{N, VT}) where {PT, N, VT}   = _diff_promote_op(PT, VT)
-_vec_diff_promote_op(::Type{PT}, ::VT, xs...) where {PT, VT}          = _diff_promote_op(PT, VT)
-_vec_diff_promote_op(::Type{PT}, xs::Tuple) where PT = _vec_diff_promote_op(PT, xs...)
-
-# even if I annotate with ::Array{_diff_promote_op(T, PolyVar{C}), N+1}, it cannot detect the type since it seems to be unable to determine the dimension N+1 :(
-function differentiate(ps::AbstractArray{PT, N}, xs::Union{AbstractArray, Tuple}) where {N, PT<:ARPL}
-    qs = Array{_vec_diff_promote_op(PT, xs), N+1}(length(xs), size(ps)...)
-    for (i, x) in enumerate(xs)
-        for j in linearindices(ps)
-            J = ind2sub(ps, j)
-            qs[i, J...] = differentiate(ps[J...], x)
-        end
-    end
-    qs
+function differentiate(ps::AbstractArray{PT}, xs::AbstractArray) where {PT <: ARPL}
+    differentiate.(reshape(ps, (1, size(ps)...)), reshape(xs, :))
+end
+
+function differentiate(ps::AbstractArray{PT}, xs::Tuple) where {PT <: ARPL}
+    differentiate.(reshape(ps, (1, size(ps)...)), xs)
 end
 
+
+# TODO: this signature is probably too wide and creates the potential
+# for stack overflows
 differentiate(p::ARPL, xs) = [differentiate(p, x) for x in xs]
 
 # differentiate(p, [x, y]) with TypedPolynomials promote x to a Monomial
 differentiate(p::ARPL, m::AbstractMonomial) = differentiate(p, variable(m))
 
-# In Julia v0.5, Base.promote_op returns Any for PolyVar, Monomial and MatPolynomial
-# Even on Julia v0.6 and Polynomial, Base.promote_op returns Any...
-_diff_promote_op(::Type{PT}, ::Type{VT}) where {PT, VT} = Base.promote_op(differentiate, PT, VT)
-_diff_promote_op(::Type{MT}, ::Type{<:AbstractVariable}) where {MT<:AbstractMonomialLike} = termtype(MT, Int)
-
-function differentiate(p, x, deg::Int)
+# The `R` argument indicates a desired result type. We use this in order
+# to attempt to preserve type-stability even though the value of `deg` cannot
+# be known at compile time. For scalar `p` and `x`, we set R to be the type
+# of differentiate(p, x) to give a stable result type regardless of `deg`. For
+# vectors p and/or x this is impossible (since differentiate may return an array),
+# so we just set `R` to `Any`
+function (_differentiate_recursive(p, x, deg::Int, ::Type{R})::R) where {R}
     if deg < 0
         throw(DomainError())
     elseif deg == 0
-        # Need the conversion with promote_op to be type stable for PolyVar, Monomial and MatPolynomial
-        return convert(_diff_promote_op(typeof(p), typeof(x)), p)
+        return p
     else
         return differentiate(differentiate(p, x), x, deg-1)
     end
 end
+
+differentiate(p, x, deg::Int) = _differentiate_recursive(p, x, deg, Base.promote_op(differentiate, typeof(p), typeof(x)))
+differentiate(p::AbstractArray, x,                              deg::Int) = _differentiate_recursive(p, x, deg, Any)
+differentiate(p,                x::Union{AbstractArray, Tuple}, deg::Int) = _differentiate_recursive(p, x, deg, Any)
+differentiate(p::AbstractArray, x::Union{AbstractArray, Tuple}, deg::Int) = _differentiate_recursive(p, x, deg, Any)
+
+
+# This is alternative, Val-based interface for nested differentiation.
+# It has the advantage of not requiring an conversion or calls to
+# Base.promote_op, while maintaining type stability for any argument
+# type.
+differentiate(p, x, ::Val{0}) = p
+differentiate(p, x, ::Val{1}) = differentiate(p, x)
+
+@static if VERSION < v"v0.7.0-"
+    # Marking this @pure helps julia v0.6 figure this out
+    Base.@pure _reduce_degree(::Val{N}) where {N} = Val{N - 1}()
+    function differentiate(p, x, deg::Val{N}) where N
+        if N < 0
+            throw(DomainError(deg))
+        else
+            differentiate(differentiate(p, x), x, _reduce_degree(deg))
+        end
+    end
+else
+    # In Julia v0.7 and above, we can remove the _reduce_degree trick
+    function differentiate(p, x, deg::Val{N}) where N
+        if N < 0
+            throw(DomainError(deg))
+        else
+            differentiate(differentiate(p, x), x, Val{N - 1}())
+        end
+    end
+end

src/division.jl

@@ -67,7 +67,7 @@ function Base.divrem(f::APL{T}, g::AbstractVector{<:APL{S}}; kwargs...) where {T
     lt = leadingterm.(g)
     rg = removeleadingterm.(g)
     lm = monomial.(lt)
-    useful = IntSet(eachindex(g))
+    useful = BitSet(eachindex(g))
     while !iszero(rf)
         ltf = leadingterm(rf)
         if isapproxzero(ltf; kwargs...)

src/monomial.jl

@@ -71,9 +71,10 @@ degree(v::AbstractVariable, var::AbstractVariable) = (v == var ? 1 : 0)
 #_deg(v::AbstractVariable) = 0
 #_deg(v::AbstractVariable, power, powers...) = v == power[1] ? power[2] : _deg(v, powers...)
 #degree(m::AbstractMonomial, v::AbstractVariable) = _deg(v, powers(t)...)
+
 function degree(m::AbstractMonomial, v::AbstractVariable)
-    i = findfirst(variables(m), v)
-    if i == nothing || iszero(i)
+    i = findfirst(equalto(v), variables(m))
+    if i === nothing || iszero(i)
         0
     else
         exponents(m)[i]

src/monovec.jl

@@ -20,7 +20,7 @@ Calling `monovec` on ``[xy, x, xy, x^2y, x]`` should return ``[x^2y, xy, x]``.
 """
 function monovec(X::AbstractVector{MT}) where {MT<:AbstractMonomial}
     Y = sort(X, rev=true)
-    dups = find(i -> Y[i] == Y[i-1], 2:length(Y))
+    dups = findall(i -> Y[i] == Y[i-1], 2:length(Y))
     deleteat!(Y, dups)
     Y
 end
@@ -54,7 +54,7 @@ Calling `sortmonovec` on ``[xy, x, xy, x^2y, x]`` should return ``([4, 1, 2], [x
 """
 function sortmonovec(X::AbstractVector{MT}) where {MT<:AbstractMonomial}
     σ = sortperm(X, rev=true)
-    dups = find(i -> X[σ[i]] == X[σ[i-1]], 2:length(σ))
+    dups = findall(i -> X[σ[i]] == X[σ[i-1]], 2:length(σ))
     deleteat!(σ, dups)
     σ, X[σ]
 end

src/operators.jl

@@ -99,14 +99,14 @@ for op in [:+, :-]
     end
 end
 
-Base.transpose(v::AbstractVariable) = v
-Base.transpose(m::AbstractMonomial) = m
-Base.transpose(t::T) where {T <: AbstractTerm} = transpose(coefficient(t)) * monomial(t)
-Base.transpose(p::AbstractPolynomialLike) = polynomial(map(transpose, terms(p)))
+adjoint_operator(v::AbstractVariable) = v
+adjoint_operator(m::AbstractMonomial) = m
+adjoint_operator(t::T) where {T <: AbstractTerm} = adjoint_operator(coefficient(t)) * monomial(t)
+adjoint_operator(p::AbstractPolynomialLike) = polynomial(map(adjoint_operator, terms(p)))
 
-Base.dot(p1::AbstractPolynomialLike, p2::AbstractPolynomialLike) = p1' * p2
-Base.dot(x, p::AbstractPolynomialLike) = x' * p
-Base.dot(p::AbstractPolynomialLike, x) = p' * x
+dot(p1::AbstractPolynomialLike, p2::AbstractPolynomialLike) = p1' * p2
+dot(x, p::AbstractPolynomialLike) = x' * p
+dot(p::AbstractPolynomialLike, x) = p' * x
 
 # Amazingly, this works! Thanks, StaticArrays.jl!
 """
@@ -115,3 +115,6 @@ The element type of the vector will be Monomial{vars, length(vars)}.
 """
 Base.vec(vars::Tuple{Vararg{AbstractVariable}}) = [vars...]
 # vec(vars::Tuple{Vararg{<:TypedVariable}}) = SVector(vars)
+
+# https://github.com/JuliaLang/julia/pull/23332
+^(x::AbstractPolynomialLike, p::Integer) = Base.power_by_squaring(x, p)

src/polynomial.jl

@@ -4,7 +4,7 @@ export mindegree, maxdegree, extdegree
 export leadingterm, leadingcoefficient, leadingmonomial
 export removeleadingterm, removemonomials, monic
 
-Base.norm(p::AbstractPolynomialLike, r::Int=2) = norm(coefficients(p), r)
+norm(p::AbstractPolynomialLike, r::Int=2) = norm(coefficients(p), r)
 
 changecoefficienttype(::Type{TT}, ::Type{T}) where {TT<:AbstractTermLike, T} = termtype(TT, T)
 changecoefficienttype(::Type{PT}, ::Type{T}) where {PT<:AbstractPolynomial, T} = polynomialtype(PT, T)