Skip to content

Instantly share code, notes, and snippets.

@mathpup
Last active October 1, 2021 12:21
Show Gist options
  • Save mathpup/8514578 to your computer and use it in GitHub Desktop.
Save mathpup/8514578 to your computer and use it in GitHub Desktop.
Extension of Polynomial.jl to support polynomial division. Also adds some handy conversions and promotion rules.
import Base.convert, Base.promote_rule
using Polynomial
convert{T<:FloatingPoint}(::Type{Poly{T}}, p::Poly) = return Poly(convert(Vector{T}, p.a))
convert{T<:Rational}(::Type{Poly{T}}, p::Poly) = Poly(convert(Vector{T}, p.a))
convert{T<:Integer}(::Type{Poly{T}}, p::Poly) = Poly(convert(Vector{T}, p.a))
promote_rule{T<:FloatingPoint, S<:FloatingPoint}(::Type{Poly{T}}, ::Type{Poly{S}}) =
Poly{promote_type(T, S)}
promote_rule{T<:Rational, S<:Rational}(::Type{Poly{T}}, ::Type{Poly{S}}) =
Poly{promote_type(T, S)}
# degree of the polynomial
deg(p::Poly) = length(p) - 1
# Leading coefficient
lc{T}(p::Poly{T}) = p[1]
# Returns the polynomial c*x^n
function cxn{T}(c::T, n::Int)
v = zeros(T, n + 1)
v[1] = c
return Poly(v)
end
# Generic Euclidean division that should work for any two polynomials
# over field T. Since polynomial operators will be called repeatedly,
# type promotion is done once at the beginning. Returns (quot, rem).
function divrem{T, S}(a::Poly{T}, b::Poly{S}, divop = /)
(A, B) = promote(a, b)
Q = Poly([zero(T)])
R = A
while R != 0 && (delta = deg(R) - deg(B)) >= 0
quot = divop(lc(R), lc(B))
T = cxn(quot, delta)
Q = Q + T
R = R - B * T
end
return (Q, R)
end
# Test code
p0 = Poly([0])
p1 = Poly([4, 2, -3, 6, 5])
p2 = Poly([6, 2, -3, 7])
p3 = p1 * p2
println("p1 * p2 = $p3")
println()
# Division with integer coefficients gives floating-point result by default
result1 = divrem(p3, p2)
println("divrem(p3, p2) = $result1")
println()
# The optional parameter // gives rational coefficients
result2 = divrem(p3, p2, //)
println("divrem(p3, p2, //) = $result2")
@jwmerrill
Copy link

Couldn't you condense the promote and convert rules to a single pair?

convert{T}(::Type{Poly{T}}, p::Poly) = Poly(convert(Vector{T}, p.a))

promote_rule{T, S}(::Type{Poly{T}}, ::Type{Poly{S}}) =
Poly{promote_type(T, S)}

Also, I think it might be more idiomatic to pull the promotion rule into a separate definition of divrem, so

divrem{T, S}(a::Poly{T}, b::Poly{S}, divop = /) = divrem(promote(a, b)..., divop)

# Note, using S here and not T because you use T to mean something else
# later in the function body
function divrem{S}(A::Poly{S}, B::Poly{S}, divop = /)
    Q = Poly([zero(S)])
    R = A

    while R != 0 && (delta = deg(R) - deg(B)) >= 0
    quot = divop(lc(R), lc(B))
    T = cxn(quot, delta)
    Q = Q + T
    R = R - B * T
    end
    return (Q, R)
end

@jwmerrill
Copy link

It might be better to just define new / and // methods rather than defining a divrem method that takes an operator.

@vtjnash
Copy link

vtjnash commented Jan 21, 2014

@jwmerrill @mathpup Do you agree to allow me to incorporate this into Polynomial.jl (MIT licensed)?

@jwmerrill
Copy link

Fine by me.

@mathpup
Copy link
Author

mathpup commented Jan 27, 2014

I am conflicted about /, //, and polydiv, and I actually went back and forth in writing the code. To make Poly consistent with subtypes of Number, polydiv could be renamed divrem, and the / and % operators could return only the quotient and remainder, respectively. My concern with this approach is one almost always wants both quotient and remainder when dividing polynomials, and it is possible that users might naively call / and % separately.

Condensing the promotion and conversion rules seems like a good idea.

@vtjnash You have had concerns about floating point precision. Something interesting I noticed in NumPy is that in their older, no longer recommended, poly1d type, polydiv does have some tolerance checking, but their newer, recommended Polynomial type does not.

@mathpup
Copy link
Author

mathpup commented Jan 27, 2014

@vtjnash You have permission to incorporate the code PolyExt.jl into Polynomial.jl under the MIT license.

If you need a more formal declaration of that, let me know what you need.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment