Issue 566a

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 @@
 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 @@
 function Base.:https://(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]