From 1fe1303ac398761c0fa7365b15be95b61347b019 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 2 May 2024 14:59:41 +0200 Subject: [PATCH] Promote integer to rational for division --- docs/src/division.md | 28 +++++++++++++++-- src/division.jl | 74 ++++++++++++++++++++++++++++++++++++++++++-- test/division.jl | 10 ++++-- 3 files changed, 105 insertions(+), 7 deletions(-) diff --git a/docs/src/division.md b/docs/src/division.md index 19fd76ce..b2c6282f 100644 --- a/docs/src/division.md +++ b/docs/src/division.md @@ -1,7 +1,5 @@ # 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)``. @@ -9,10 +7,34 @@ 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 +``` diff --git a/src/division.jl b/src/division.jl index feedeb47..4cc5ab15 100644 --- a/src/division.jl +++ b/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 diff --git a/test/division.jl b/test/division.jl index 36b995ab..cce74b1e 100644 --- a/test/division.jl +++ b/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