Rational function (#335)

* add rational-function type; fix bug with ngcd; plot recipe

* fix bugs with ngcd; add tests; rename to lowest_terms

* use VERSION check for rational functions

Project.toml

@@ -2,7 +2,7 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "2.0.8"
+version = "2.0.9"
 
 [deps]
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"

src/Polynomials.jl

@@ -29,6 +29,16 @@ include("polynomials/multroot.jl")
 
 include("polynomials/ChebyshevT.jl")
 
+# Rational functions
+if VERSION >= v"1.2.0"
+    include("rational-functions/common.jl")
+    include("rational-functions/rational-function.jl")
+    include("rational-functions/fit.jl")
+    #include("rational-transfer-function.jl")
+    include("rational-functions/plot-recipes.jl")
+end
+
+
 # compat; opt-in with `using Polynomials.PolyCompat`
 include("polynomials/Poly.jl")
 

src/common.jl

@@ -1019,6 +1019,7 @@ function Base.isapprox(p1::AbstractPolynomial{T,X},
     rtol::Real = (Base.rtoldefault(T, S, 0)),
                        atol::Real = 0,) where {T,X,S,Y}
     assert_same_variable(p1, p2)
+    (hasnan(p1) || hasnan(p2)) && return false  # NaN poisons comparisons
     # copy over from abstractarray.jl
     Δ  = norm(p1-p2)
     if isfinite(Δ)

src/polynomials/multroot.jl

@@ -79,6 +79,9 @@ is misidentified.
 function multroot(p::Polynomials.StandardBasisPolynomial{T}; verbose=false,
                   kwargs...) where {T}
 
+    # degenerate case, constant
+    degree(p) == 0 && return (values=T[], multiplicities=Int[], κ=NaN, ϵ=NaN)
+
     # degenerate case, all zeros
     if (nz = findfirst(!iszero, coeffs(p))) == length(coeffs(p))
         return (values=zeros(T,1), multiplicities=[nz-1], κ=NaN, ϵ=NaN)
@@ -202,7 +205,7 @@ function pejorative_root(p, zs::Vector{S}, ls::Vector{Int};
 The multiplicity count may be in error: the initial guess for the roots failed
 to converge to a pejorative root.
 """)
-        return(zₘs)
+        return(zₖs)
     end
 
 end

src/polynomials/ngcd.jl

@@ -19,6 +19,7 @@ function ngcd(p::P, q::Q,
         a,b = divrem(p,q)
         return ngcd(q,b, args...; λ=100, kwargs...)
     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)
@@ -28,15 +29,30 @@ function ngcd(p::P, q::Q,
     Polynomials.assert_same_variable(p,q)
 
     R = promote_type(float(T), float(S))
+    𝑷 = Polynomials.constructorof(promote_type(P,Q)){R,X}
+
     ps = R[pᵢ for pᵢ ∈ coeffs(p)]
     qs = R[qᵢ for qᵢ ∈ coeffs(q)]
-    p′ = PnPolynomial(ps)
-    q′ = PnPolynomial(qs)
 
+    # cancel zeros
+    nz = min(findfirst(!iszero, ps), findfirst(!iszero, qs))
+    if nz == length(qs)
+        x = variable(p)
+        u = x^(nz-1)
+        v,w = 𝑷(ps[nz:end]), 𝑷(qs[nz:end])
+        return (u=u, v=v, w=w, Θ=NaN, κ=NaN)
+    end
+
+    ## call ngcd
+    p′ = PnPolynomial{R,X}(ps[nz:end])
+    q′ = PnPolynomial{R,X}(qs[nz:end])
     out = NGCD.ngcd(p′, q′, args...; kwargs...)
 
-    𝑷 = Polynomials.constructorof(promote_type(P,Q)){R,X} 
+    𝑷 = Polynomials.constructorof(promote_type(P,Q)){R,X}
     u,v,w = convert.(𝑷, (out.u,out.v,out.w))
+    if nz > 1
+        u *= variable(u)^(nz-1)
+    end
     (u=u,v=v,w=w, Θ=out.Θ, κ = out.κ)
 
 end