Extension of Polynomial.jl to support polynomial division. Also adds some handy conversions and promotion rules.

PolyExt.jl

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

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

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

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

Fine by me.

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.

@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.