diff --git a/docs/src/division.md b/docs/src/division.md
index 19fd76ce..d126550b 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,41 @@ 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.
+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
+AbstractUnivariateGCDAlgorithm
+GeneralizedEuclideanAlgorithm
+SubresultantAlgorithm
+```
+Internal functions of the `gcd` algorithm:
+```@docs
+isolate_variable
+primitive_univariate_gcd!
+univariate_gcd
+content
+primitive_part
+primitive_part_content
+```
diff --git a/src/division.jl b/src/division.jl
index feedeb47..889064b9 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
@@ -152,6 +179,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 +236,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.