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")
@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