Promote integer to rational for division

docs/src/division.md

@@ -1,18 +1,40 @@
 # Division
 
-The `gcd` and `lcm` functions of `Base` have been implemented for monomials, you have for example `gcd(x^2*y^7*z^3, x^4*y^5*z^2)` returning `x^2*y^5*z^2` and `lcm(x^2*y^7*z^3, x^4*y^5*z^2)` returning `x^4*y^7*z^3`.
-
 Given two polynomials, ``p`` and ``d``, there are unique ``r`` and ``q`` such that ``p = q d + r`` and the leading term of ``d`` does not divide the leading term of ``r``.
 You can obtain ``q`` using the `div` function and ``r`` using the `rem` function.
 The `divrem` function returns ``(q, r)``.
 
 Given a polynomial ``p`` and divisors ``d_1, \ldots, d_n``, one can find ``r`` and ``q_1, \ldots, q_n`` such that ``p = q_1 d_1 + \cdots + q_n d_n + r`` and none of the leading terms of ``q_1, \ldots, q_n`` divide the leading term of ``r``.
 You can obtain the vector ``[q_1, \ldots, q_n]`` using `div(p, d)` where ``d = [d_1, \ldots, d_n]`` and ``r`` using the `rem` function with the same arguments.
 The `divrem` function returns ``(q, r)``.
-
 ```@docs
+divrem
+div
+rem
 divides
 div_multiple
+```
+
+Note that the coefficients of the polynomials need to be a field for `div`,
+`rem` and `divrem` to work.
+If the coefficient type is not a field, it is promoted to a field using [`promote_to_field`](@ref).
+```@docs
+promote_to_field
+```
+Alternatively, [`pseudo_rem`](@ref) or [`pseudo_divrem`](@ref) can be used
+instead as they do not require the coefficient type to be a field.
+```@docs
 pseudo_rem
+pseudo_divrem
 rem_or_pseudo_rem
 ```
+
+## Greatest Common Divisor (GCD)
+
+The Greatest Common Divisor (GCD) and Least Common Multiple (LCM) can be
+obtained for integers respectively with the `gcd` and `lcm` functions.
+The same functions can be used with monomials and polynomials:
+```@docs
+gcd
+lcm
+```

src/division.jl

@@ -26,9 +26,36 @@ function divides(t1::AbstractTermLike, t2::AbstractTermLike)
 end
 divides(t1::AbstractVariable, t2::AbstractVariable) = t1 == t2
 
+"""
+    gcd(m1::AbstractMonomialLike, m2::AbstractMonomialLike)
+
+Return the largest monomial `m` such that both `divides(m, m1)`
+and `divides(m, m2)` are `true`.
+
+```@example
+julia> @polyvar x y z;
+
+julia> gcd(x^2*y^7*z^3, x^4*y^5*z^2)
+x²y⁵z²
+```
+"""
 function Base.gcd(m1::AbstractMonomialLike, m2::AbstractMonomialLike)
     return map_exponents(min, m1, m2)
 end
+
+"""
+    lcm(m1::AbstractMonomialLike, m2::AbstractMonomialLike)
+
+Return the smallest monomial `m` such that both `divides(m1, m)`
+and `divides(m2, m)` are `true`.
+
+```@example
+julia> @polyvar x y z;
+
+julia> lcm(x^2*y^7*z^3, x^4*y^5*z^2)
+x^4*y^7*z^3
+```
+"""
 function Base.lcm(m1::AbstractMonomialLike, m2::AbstractMonomialLike)
     return map_exponents(max, m1, m2)
 end
@@ -37,8 +64,25 @@ struct Field end
 struct UniqueFactorizationDomain end
 const UFD = UniqueFactorizationDomain
 
+"""
+    promote_to_field(::Type{T})
+
+Promote the type `T` to a field. For instance, `promote_to_field(T)` returns
+`Rational{T}` if `T` is an integer and `promote_to_field(T)` returns `RationalPoly{T}`
+if `T` is a polynomial.
+"""
+function promote_to_field end
+
+function promote_to_field(::Type{T}) where {T<:Integer}
+    return Rational{T}
+end
+function promote_to_field(::Type{T}) where {T<:_APL}
+    return RationalPoly{T,T}
+end
+promote_to_field(::Type{T}) where {T} = T
+
 algebraic_structure(::Type{<:Integer}) = UFD()
-algebraic_structure(::Type{<:AbstractPolynomialLike}) = UFD()
+algebraic_structure(::Type{<:_APL}) = UFD()
 # `Rational`, `AbstractFloat`, JuMP expressions, etc... are fields
 algebraic_structure(::Type) = Field()
 _field_absorb(::UFD, ::UFD) = UFD()
@@ -152,6 +196,26 @@ function Base.rem(f::_APL, g::Union{_APL,AbstractVector{<:_APL}}; kwargs...)
     return divrem(f, g; kwargs...)[2]
 end
 
+"""
+    pseudo_divrem(f::_APL{S}, g::_APL{T}, algo) where {S,T}
+
+Return the pseudo divisor and remainder of `f` modulo `g` as defined in [Knu14, Algorithm R, p. 425].
+
+When the coefficient type is not a field, it is not always possible to carry a
+division. For instance, the division of `f = 3x + 1` by `g = 2x + 1` cannot be done over
+integers. On the other hand, one can write `2f = 3g - 1`.
+In general, the *pseudo* division of `f` by `g` is:
+```math
+l f(x) = q(x) g(x) + r(x)
+```
+where `l` is a power of the leading coefficient of `g` some constant.
+
+See also [`pseudo_rem`](@ref).
+
+[Knu14] Knuth, D.E., 2014.
+*Art of computer programming, volume 2: Seminumerical algorithms.*
+Addison-Wesley Professional. Third edition.
+"""
 function pseudo_divrem(f::_APL{S}, g::_APL{T}, algo) where {S,T}
     return _pseudo_divrem(
         algebraic_structure(MA.promote_operation(-, S, T)),
@@ -189,6 +253,8 @@ end
 
 Return the pseudo remainder of `f` modulo `g` as defined in [Knu14, Algorithm R, p. 425].
 
+See [`pseudo_divrem`](@ref) for more details.
+
 [Knu14] Knuth, D.E., 2014.
 *Art of computer programming, volume 2: Seminumerical algorithms.*
 Addison-Wesley Professional. Third edition.
@@ -381,7 +447,11 @@ function MA.promote_operation(
     ::Type{P},
     ::Type{Q},
 ) where {T,S,P<:_APL{T},Q<:_APL{S}}
-    U = MA.promote_operation(/, T, S)
+    U = MA.promote_operation(
+        /,
+        promote_to_field(T),
+        promote_to_field(S),
+    )
     # `promote_type(P, Q)` is needed for TypedPolynomials in case they use different variables
     return polynomial_type(promote_type(P, Q), MA.promote_operation(-, U, U))
 end

test/division.jl

@@ -54,8 +54,14 @@ end
 
 function divrem_test()
     Mod.@polyvar x y
-    @test (@inferred div(x * y^2 + 1, x * y + 1)) == y
-    @test (@inferred rem(x * y^2 + 1, x * y + 1)) == -y + 1
+    p = x * y^2 + 1
+    q = x * y + 1
+    @test typeof(div(p, q)) == MA.promote_operation(div, typeof(p), typeof(q))
+    @test typeof(rem(p, q)) == MA.promote_operation(rem, typeof(p), typeof(q))
+    @test coefficient_type(div(p, q)) == Rational{Int}
+    @test coefficient_type(rem(p, q)) == Rational{Int}
+    @test (@inferred div(p, q)) == y
+    @test (@inferred rem(p, q)) == -y + 1
     @test (@inferred div(x * y^2 + x, y)) == x * y
     @test (@inferred rem(x * y^2 + x, y)) == x
     @test (@inferred rem(x^4 + x^3 + (1 + 1e-10) * x^2 + 1, x^2 + x + 1)) == 1