Issue 566a (#570)

* laurent and rational type issue

* use Laurent for promotion

Author: jverzani

src/polynomial-container-types/mutable-dense-laurent-polynomial.jl

@@ -204,6 +204,10 @@ function offset_vector_combine(op, p::MutableDenseLaurentPolynomial{B,T,X}, q::M
 
 end
 
+# multiply by xⁿ
+_shift(p::P, n::Int) where {P <: MutableDenseLaurentPolynomial} =
+    P(Val(false), p.coeffs, firstindex(p) + n)
+
 function Base.numerator(q::MutableDenseLaurentPolynomial{B,T,X}) where {B,T,X}
     p = chop(q)
     o = firstindex(p)

src/polynomials/standard-basis/laurent-polynomial.jl

@@ -99,6 +99,10 @@ function evalpoly(c, p::LaurentPolynomial{T,X}) where {T,X}
     EvalPoly.evalpoly(c, p.coeffs) * c^p.order[]
 end
 
+# ---
+
+
+
 # scalar add
 function scalar_add(c::S, p::LaurentPolynomial{T,X}) where {S, T, X}
     R = promote_type(T,S)

src/promotions.jl

@@ -15,11 +15,11 @@ Base.promote_rule(::Type{P}, ::Type{Q}) where {B,T,S,X,
                                                Q<:AbstractLaurentUnivariatePolynomial{B,S,X}} = MutableDenseLaurentPolynomial{B,promote_type(T,S),X}
 
 
-
 # Methods to ensure that matrices of polynomials behave as desired
 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}
+                                                   LaurentPolynomial{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}} =

src/rational-functions/common.jl

@@ -49,21 +49,6 @@ function Base.convert(::Type{PQ}, p::Number) where {PQ <: AbstractRationalFuncti
    rational_function(PQ, p * one(P), one(P))
 end
 
-function Base.convert(::Type{PQ}, q::LaurentPolynomial) where {PQ <: AbstractRationalFunction}
-    m = firstindex(q)
-    if m >= 0
-        p = convert(Polynomial, q)
-        return convert(PQ, p)
-    else
-        z = variable(q)
-        zᵐ = z^(-m)
-        p = convert(Polynomial, zᵐ * q)
-        return rational_function(PQ, p, convert(Polynomial, zᵐ))
-    end
-
-end
-
-
 function Base.convert(::Type{PQ}, p::P) where {PQ <: AbstractRationalFunction, P<:AbstractPolynomial}
     T′ = _eltype(_eltype((PQ)))
     T =  isnothing(T′) ? eltype(p) : T′
@@ -74,9 +59,6 @@ function Base.convert(::Type{PQ}, p::P) where {PQ <: AbstractRationalFunction, P
     rational_function(PQ, 𝐩, one(𝐩))
 end
 
-
-
-
 function Base.convert(::Type{P}, pq::PQ) where {P<:AbstractPolynomial, PQ<:AbstractRationalFunction}
     p,q = pqs(pq)
     isconstant(q) || throw(ArgumentError("Can't convert rational function with non-constant denominator to a polynomial."))
@@ -97,7 +79,8 @@ function Base.promote_rule(::Type{PQ}, ::Type{PQ′}) where {T,X,P,PQ <: Abstrac
     assert_same_variable(X,X′)
     PQ_, PQ′_ = constructorof(PQ), constructorof(PQ′)
     𝑷𝑸 = PQ_ == PQ′ ? PQ_ : RationalFunction
-    𝑷 = constructorof(typeof(variable(P)+variable(P′)));  𝑷 = Polynomial
+    𝑷 = constructorof(typeof(variable(P)+variable(P′)));
+    #𝑷 = Polynomial
     𝑻 = promote_type(T,T′)
     𝑷𝑸{𝑻,X,𝑷{𝑻,X}}
 end

src/rational-functions/rational-function.jl

@@ -47,25 +47,51 @@ julia> derivative(pq)
 struct RationalFunction{T, X, P<:AbstractPolynomial{T,X}} <: AbstractRationalFunction{T,X,P}
     num::P
     den::P
-    function RationalFunction(p::P, q::P) where {T,X, P<:AbstractPolynomial{T,X}}
+    function RationalFunction{T,X,P}(p, q) where {T,X, P<:AbstractPolynomial{T,X}}
         new{T,X,P}(p, q)
     end
-    function RationalFunction(p::P, q::T) where {T,X, P<:AbstractPolynomial{T,X}}
-        new{T,X,P}(p, q*one(P))
-    end
-    function RationalFunction(p::T, q::Q) where {T,X, Q<:AbstractPolynomial{T,X}}
-        new{T,X,Q}(p*one(Q), q)
+
+end
+
+function RationalFunction(p::P, q::P) where {T,X, P<:AbstractPolynomial{T,X}}
+    RationalFunction{T,X,P}(p, q)
+end
+
+RationalFunction(p::AbstractPolynomial{T,X}, q::AbstractPolynomial{S,X}) where {T,S,X} =
+    RationalFunction(promote(p,q)...)
+
+function RationalFunction(p::P, q::T) where {T,X, P<:AbstractPolynomial{T,X}}
+    RationalFunction(p, (q * one(p)))
+end
+function RationalFunction(p::T, q::Q) where {T,X, Q<:AbstractPolynomial{T,X}}
+    RationalFunction(p * one(q),  q)
+end
+
+function RationalFunction(p::P,q::P) where {T, X, P <: LaurentPolynomial{T,X}}
+
+    m,n = firstindex(p), firstindex(q)
+    p′,q′ = _shift(p, -m), _shift(q, -n)
+    if m-n ≥ 0
+        return RationalFunction{T,X,P}(_shift(p′, m-n), q′)
+    else
+        return RationalFunction{T,X,P}(p′, _shift(q′, n-m))
     end
 end
 
-RationalFunction(p,q)  = RationalFunction(convert(LaurentPolynomial,p), convert(LaurentPolynomial,q))
-function RationalFunction(p::LaurentPolynomial,q::LaurentPolynomial)
-    𝐩 = convert(RationalFunction, p)
-    𝐪 = convert(RationalFunction, q)
-    𝐩 // 𝐪
+# RationalFunction(p,q)  = RationalFunction(convert(LaurentPolynomial,p), convert(LaurentPolynomial,q))
+
+
+# special case Laurent
+function lowest_terms(pq::PQ; method=:numerical, kwargs...) where {T,X,
+                                                                   P<:LaurentPolynomial{T,X}, #StandardBasisPolynomial{T,X},
+                                                                   PQ<:AbstractRationalFunction{T,X,P}}
+    p,q = pqs(pq)
+    p′,q′ = convert(Polynomial, p), convert(Polynomial,q)
+    u,v,w = uvw(p′,q′; method=method, kwargs...)
+    v′,w′ = convert(LaurentPolynomial, v), convert(LaurentPolynomial, w)
+    rational_function(PQ, v′/w′[end], w′/w′[end])
 end
-RationalFunction(p::LaurentPolynomial,q::Number) = convert(RationalFunction, p) // q
-RationalFunction(p::Number,q::LaurentPolynomial) = q // convert(RationalFunction, p)
+
 
 RationalFunction(p::AbstractPolynomial) = RationalFunction(p,one(p))
 
@@ -76,3 +102,14 @@ RationalFunction(p::AbstractPolynomial) = RationalFunction(p,one(p))
 function Base.://(p::AbstractPolynomial,q::AbstractPolynomial)
     RationalFunction(p,q)
 end
+
+
+# promotion
+function Base.promote(pq::RationalFunction{T,X,P},
+                      rs::RationalFunction{S,X,Q}) where {T,S,X,P<:AbstractPolynomial{T,X},
+                                                          Q<:AbstractPolynomial{S,X}}
+    𝑃 = promote_type(P, Q)
+    p,q = pq
+    r,s = rs
+    (convert(𝑃,p) // convert(𝑃,q), convert(𝑃,r) // convert(𝑃,s))
+end

test/rational-functions.jl

@@ -11,7 +11,7 @@ using LinearAlgebra
     # constructor
     @test p // q isa RationalFunction
     @test p // r isa RationalFunction
-    @test_throws ArgumentError r // s
+    @test_throws MethodError r // s
     @test RationalFunction(p) == p // one(p)
 
     pq = p // t # promotes to type of t
@@ -76,6 +76,21 @@ using LinearAlgebra
     @test numerator(p) == p * variable(p)^2
     @test denominator(p) == convert(Polynomial, variable(p)^2)
 
+    # issue 566
+    q = LaurentPolynomial([1], -1)
+    p = LaurentPolynomial([1], 1)
+    @test degree(numerator(q // p)) == 0 # q/p = 1/x^2
+    @test degree(denominator(q // p)) == 2
+
+    @test degree(numerator(p // q)) == 2 # p/q = x^2 / 1
+    @test degree(denominator(p // q)) == 0
+
+    @test degree(numerator(q // q^2)) == 1
+    @test degree(denominator(q // q^2)) == 0
+
+    @test degree(numerator(q^2 // q)) == 0
+    @test degree(denominator(q^2 // q)) == 1
+
 end
 
 @testset "zeros, poles, residues" begin
@@ -112,7 +127,7 @@ end
     p, q = Polynomial([1,2], :x), Polynomial([1,2],:y)
     pp = p // (p-1)
     PP = typeof(pp)
-
+    PP′ = RationalFunction{Int64, :x, LaurentPolynomial{Int64, :x}}
     r, s = SparsePolynomial([1,2], :x), SparsePolynomial([1,2],:y)
     rr = r // (r-1)
 
@@ -124,7 +139,7 @@ end
     # @test eltype(eltype(eltype([im, p, pp]))) == Complex{Int}
 
     ## test mixed types promote polynomial type
-    @test eltype([pp rr p r]) == PP
+    @test eltype([pp rr p r]) == PP′ # promotes to LaurentPolynomial
 
     ## test non-constant polys of different symbols throw error
     @test_throws ArgumentError [pp, q]