From d800b2831f6158697258e9508d0fb647e51d834a Mon Sep 17 00:00:00 2001 From: jverzani Date: Sun, 16 Mar 2025 10:15:45 -0400 Subject: [PATCH 1/3] one part of #602 --- src/rational-functions/common.jl | 1 + test/rational-functions.jl | 12 ++++++++++++ 2 files changed, 13 insertions(+) diff --git a/src/rational-functions/common.jl b/src/rational-functions/common.jl index f82ccb48..d254d9aa 100644 --- a/src/rational-functions/common.jl +++ b/src/rational-functions/common.jl @@ -459,6 +459,7 @@ function roots(pq::AbstractRationalFunction{T}; kwargs...) where {T} pq′ = lowest_terms(pq; method=method, kwargs...) den = numerator(pq′) + (hasnan(den) || any(isinf, den)) && throw(ArgumentError("Reduced expression has numeric issues")) mmethod = something(multroot_method, default_multroot_method(T)) mr = Multroot.multroot(den; method=mmethod) (zs=mr.values, multiplicities = mr.multiplicities) diff --git a/test/rational-functions.jl b/test/rational-functions.jl index 7af60ba6..4de70cd6 100644 --- a/test/rational-functions.jl +++ b/test/rational-functions.jl @@ -124,6 +124,18 @@ end end @test norm(numerator(lowest_terms(d - pq)), Inf) <= sqrt(eps()) + + ## issue #602 + s = Polynomial([0,1],:s) + r = (15s^14 + 1e-16s^15)//(s) + @test_throws ArgumentError roots(r) + + r = (15s^14 + 1e-16s^15)//(1s + 14s^14) + num = numerator(Polynomials.lowest_terms(r)) + @test length(roots(num)) == 14 + @test_throws DimensionMismatch Polynomials.Multroot.multroot(num) # where to fix + + end @testset "As matrix elements" begin From 6619d47835bfa2facbcd93a365dcbdd5f3902724 Mon Sep 17 00:00:00 2001 From: jverzani Date: Sun, 16 Mar 2025 10:55:16 -0400 Subject: [PATCH 2/3] add deflation to ngcd --- src/polynomials/ngcd.jl | 7 +++++++ test/rational-functions.jl | 4 +++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/polynomials/ngcd.jl b/src/polynomials/ngcd.jl index bdb3e494..e2a93d84 100644 --- a/src/polynomials/ngcd.jl +++ b/src/polynomials/ngcd.jl @@ -13,6 +13,13 @@ function ngcd(p::P, q::Q, kwargs...) where {T,X,P<:StandardBasisPolynomial{T,X}, S,Y,Q<:StandardBasisPolynomial{S,Y}} + # deflate leading zeroes when present in top and bottom + i,j = findfirst.(!iszero, (p,q)) + if !isnothing(i) && !isnothing(j) + k = min(i,j) + p,q = P(coeffs(p)[k+1:end]), Q(coeffs(q)[k+1:end]) + end + # easy cases degree(p) < 0 && return (u=q, v=p, w=one(q), θ=NaN, κ=NaN) degree(p) == 0 && return (u=one(q), v=p, w=q, θ=NaN, κ=NaN) diff --git a/test/rational-functions.jl b/test/rational-functions.jl index 4de70cd6..d38da278 100644 --- a/test/rational-functions.jl +++ b/test/rational-functions.jl @@ -128,7 +128,9 @@ end ## issue #602 s = Polynomial([0,1],:s) r = (15s^14 + 1e-16s^15)//(s) - @test_throws ArgumentError roots(r) + num = numerator(Polynomials.lowest_terms(r)) + @test length(roots(num)) == 14 + @test_throws DimensionMismatch Polynomials.Multroot.multroot(num) # where to fix r = (15s^14 + 1e-16s^15)//(1s + 14s^14) num = numerator(Polynomials.lowest_terms(r)) From 88bb3b6674b22104ab709793931e32a330784016 Mon Sep 17 00:00:00 2001 From: jverzani Date: Sun, 16 Mar 2025 11:36:23 -0400 Subject: [PATCH 3/3] move where deflation happens --- src/polynomials/ngcd.jl | 39 ++++++++++++++++++++------------------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/src/polynomials/ngcd.jl b/src/polynomials/ngcd.jl index e2a93d84..32ed322b 100644 --- a/src/polynomials/ngcd.jl +++ b/src/polynomials/ngcd.jl @@ -13,12 +13,6 @@ function ngcd(p::P, q::Q, kwargs...) where {T,X,P<:StandardBasisPolynomial{T,X}, S,Y,Q<:StandardBasisPolynomial{S,Y}} - # deflate leading zeroes when present in top and bottom - i,j = findfirst.(!iszero, (p,q)) - if !isnothing(i) && !isnothing(j) - k = min(i,j) - p,q = P(coeffs(p)[k+1:end]), Q(coeffs(q)[k+1:end]) - end # easy cases degree(p) < 0 && return (u=q, v=p, w=one(q), θ=NaN, κ=NaN) @@ -26,20 +20,11 @@ function ngcd(p::P, q::Q, degree(q) < 0 && return (u=one(q), v=p, w=zero(q), θ=NaN, κ=NaN) degree(q) == 0 && return (u=one(p), v=p, w=q, Θ=NaN, κ=NaN) + if (degree(q) > degree(p)) u,w,v,Θ,κ = ngcd(q,p,args...; kwargs...) return (u=u,v=v,w=w, Θ=Θ, κ=κ) end - if degree(p) > 5*(1+degree(q)) - a,b = divrem(p,q) - return ngcd(q, b, args...; λ=100, kwargs...) - end - - # other easy cases - p ≈ q && return (u=p,v=one(p), w=one(p), θ=NaN, κ=NaN) - Polynomials.assert_same_variable(p,q) - - R = promote_type(float(T)) 𝑷 = Polynomials.constructorof(P){R,X} @@ -52,14 +37,30 @@ function ngcd(p::P, q::Q, if nz == length(qs) x = variable(p) u = x^(nz-1) - v,w = 𝑷(ps[nz:end]), 𝑷(qs[nz:end]) + ps,qs = ps[nz:end], qs[nz:end] + v,w = 𝑷(ps), 𝑷(qs) return (u=u, v=v, w=w, Θ=NaN, κ=NaN) + elseif 1 <= nz < length(qs) + ps,qs = ps[nz:end], qs[nz:end] + p,q = P(ps), Q(qs) end + if degree(p) > 5*(1+degree(q)) + a,b = divrem(p,q) + return ngcd(q, b, args...; λ=100, kwargs...) + end + + # other easy cases + p ≈ q && return (u=p,v=one(p), w=one(p), θ=NaN, κ=NaN) + Polynomials.assert_same_variable(p,q) + + @show ps, qs + ## call ngcd P′ = PnPolynomial - p′ = P′{R,X}(ps[nz:end]) - q′ = P′{R,X}(qs[nz:end]) + p′ = P′{R,X}(ps) + q′ = P′{R,X}(qs) + out = NGCD.ngcd(p′, q′, args...; kwargs...) ## convert to original polynomial type 𝑷 = Polynomials.constructorof(P){R,X}