From 87f02dbebefc7dd71fe5fc132405252fc3efd447 Mon Sep 17 00:00:00 2001 From: john verzani Date: Fri, 10 May 2024 10:59:19 -0400 Subject: [PATCH] Issue 566a (#570) * laurent and rational type issue * use Laurent for promotion --- .../mutable-dense-laurent-polynomial.jl | 4 ++ .../standard-basis/laurent-polynomial.jl | 4 ++ src/promotions.jl | 4 +- src/rational-functions/common.jl | 21 +------ src/rational-functions/rational-function.jl | 63 +++++++++++++++---- test/rational-functions.jl | 21 ++++++- 6 files changed, 80 insertions(+), 37 deletions(-) diff --git a/src/polynomial-container-types/mutable-dense-laurent-polynomial.jl b/src/polynomial-container-types/mutable-dense-laurent-polynomial.jl index a8513ae7..04dfc948 100644 --- a/src/polynomial-container-types/mutable-dense-laurent-polynomial.jl +++ b/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) diff --git a/src/polynomials/standard-basis/laurent-polynomial.jl b/src/polynomials/standard-basis/laurent-polynomial.jl index 55934913..504ad0fc 100644 --- a/src/polynomials/standard-basis/laurent-polynomial.jl +++ b/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) diff --git a/src/promotions.jl b/src/promotions.jl index 50881568..d188dcad 100644 --- a/src/promotions.jl +++ b/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}} = diff --git a/src/rational-functions/common.jl b/src/rational-functions/common.jl index 472de718..dae8bb5a 100644 --- a/src/rational-functions/common.jl +++ b/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 diff --git a/src/rational-functions/rational-function.jl b/src/rational-functions/rational-function.jl index 918c7195..6c039f92 100644 --- a/src/rational-functions/rational-function.jl +++ b/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 diff --git a/test/rational-functions.jl b/test/rational-functions.jl index 391af0e9..31078cdd 100644 --- a/test/rational-functions.jl +++ b/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]