Add generic is_unit and is_nilpotent (#1933)

src/MPoly.jl

@@ -425,16 +425,22 @@ end
 
 iszero(x::MPolyRingElem{T}) where T <: RingElement = length(x) == 0
 
-function is_unit(a::MPolyRingElem{T}) where T <: RingElement
-   if is_constant(a)
-      return is_unit(leading_coefficient(a))
-   elseif is_domain_type(elem_type(coefficient_ring(a)))
-      return false
-   elseif length(a) == 1
-      return false
-   else
-      throw(NotImplementedError(:is_unit, a))
-   end
+function is_unit(f::T) where {T <: MPolyRingElem}
+  # for constant polynomials we delegate to the coefficient ring:
+  is_constant(f) && return is_unit(constant_coefficient(f))
+  # Here deg(f) > 0; over an integral domain, non-constant polynomials are never units:
+  is_domain_type(T) && return false
+  # A polynomial over a commutative ring is a unit iff its
+  # constant term is a unit and all other coefficients are nilpotent:
+  # see e.g. <https://kconrad.math.uconn.edu/blurbs/ringtheory/polynomial-properties.pdf> for a proof.
+  for (c, expv) in zip(coefficients(f), exponent_vectors(f))
+    if is_zero(expv)
+      is_unit(c) || return false
+    else
+      is_nilpotent(c) || return false
+    end
+  end
+  return true
 end
 
 function content(a::MPolyRingElem{T}) where T <: RingElement
@@ -445,8 +451,16 @@ function content(a::MPolyRingElem{T}) where T <: RingElement
    return z
 end
 
+
+function is_nilpotent(f::T) where {T <: MPolyRingElem}
+  is_domain_type(T) && return is_zero(f)
+  return all(is_nilpotent, coefficients(f))
+end
+
+
 function is_zero_divisor(x::MPolyRingElem{T}) where T <: RingElement
-   return is_zero_divisor(content(x))
+  is_domain_type(T) && return is_zero(x)
+  return is_zero_divisor(content(x))
 end
 
 function is_zero_divisor_with_annihilator(a::MPolyRingElem{T}) where T <: RingElement

src/NCPoly.jl

@@ -783,3 +783,34 @@ polynomial_ring_only(R::T, s::Symbol; cached::Bool=true) where T<:NCRing =
 # Simplified constructor
 
 PolyRing(R::NCRing) = polynomial_ring_only(R, :x; cached=false)
+
+
+
+###############################################################################
+#
+#   is_unit & is_nilpotent
+#
+###############################################################################
+
+# ASSUMES structural interface is analogous to that for univariate polynomials
+
+# This function handles both PolyRingElem & NCPolyRingElem
+function is_unit(f::T) where {T <: PolynomialElem}
+  # constant coeff must itself be a unit
+  is_unit(constant_coefficient(f)) || return false
+  is_constant(f) && return true
+  # Here deg(f) > 0; over an integral domain, non-constant polynomials are never units:
+  is_domain_type(T) && return false
+  for i in 1:degree(f) # we have already checked coeff(f,0)
+    if !is_nilpotent(coeff(f, i))
+      return false
+    end
+  end
+  return true
+end
+
+# This function handles both PolyRingElem & NCPolyRingElem
+function is_nilpotent(f::T) where {T <: PolynomialElem}
+  is_domain_type(T) && return is_zero(f)
+  return all(is_nilpotent, coefficients(f))
+end

src/NCRings.jl

@@ -137,7 +137,7 @@ Base.literal_pow(::typeof(^), x::NCRingElem, ::Val{p}) where {p} = x^p
 ###############################################################################
 
 @doc raw"""
-    is_unit(a::T) where {T <: NCRingElem}
+    is_unit(a::T) where {T <: NCRingElement}
 
 Return true if $a$ is invertible, else return false.
 
@@ -155,6 +155,27 @@ julia> is_unit(ZZ(-1)), is_unit(ZZ(4))
 """
 function is_unit end
 
+@doc raw"""
+    is_nilpotent(a::T) where {T <: NCRingElement}
+
+Return true iff $a$ is nilpotent, i.e. a^k == 0 for some k.
+
+# Examples
+```jldoctest
+julia> R, _ = residue_ring(ZZ,720);
+
+julia> S, x = polynomial_ring(R, :x);
+
+julia> is_nilpotent(30*x), is_nilpotent(30+90*x), is_nilpotent(S(15))
+(true, true, false)
+```
+"""
+function is_nilpotent(a::T) where {T <: NCRingElement}
+  is_domain_type(T) && return is_zero(a)
+  throw(NotImplementedError(:is_nilpotent, a))
+end
+
+
 ###############################################################################
 #
 #   Characteristic

src/Poly.jl

@@ -222,17 +222,9 @@ function is_monic(a::PolynomialElem)
     return isone(leading_coefficient(a))
 end
 
-function is_unit(a::PolynomialElem)
-   if length(a) <= 1
-      return is_unit(coeff(a, 0))
-   elseif is_domain_type(elem_type(coefficient_ring(a)))
-      return false
-   elseif !is_unit(coeff(a, 0)) || is_unit(coeff(a, length(a) - 1))
-      return false
-   else
-      throw(NotImplementedError(:is_unit, a))
-   end
-end
+# function is_unit(...)      see NCPoly.jl
+# function is_nilpotent(...) see NCPoly.jl
+
 
 is_zero_divisor(a::PolynomialElem) = is_zero_divisor(content(a))
 

src/Residue.jl

@@ -85,6 +85,19 @@ function is_unit(a::ResElem)
    return isone(g)
 end
 
+function is_nilpotent(res::ResElem)
+  m = modulus(res)
+  r = data(res)
+  while true
+    g = gcd(r, m)
+    (g == m) && return true
+    is_one(g) && return false
+    m = divexact(m, g)
+    g = mod(g, m);   r = g^2  # if computation domain is limited precision integer then g = mod(g,m) guarantees that g^2 will not overflow!
+  end
+end
+
+
 # currently residue rings are only allowed over domains
 # otherwise this function would be more complicated
 is_zero_divisor(a::ResElem) = !is_unit(a)

src/algorithms/LaurentPoly.jl

@@ -138,12 +138,28 @@ function is_monomial_recursive(p::LaurentPolyRingElem)
    is_monomial_recursive(coeff(p, dr[]))
 end
 
-function is_unit(p::LaurentPolyRingElem)
-   dr = degrees_range(p)
-   length(dr) == 1 || return false
-   is_unit(coeff(p, dr[]))
+function is_unit(f::T) where {T <: LaurentPolyRingElem}
+  # **NOTE**  f.poly is not normalized so that the degree 0 coeff is non-zero
+  is_trivial(parent(f)) && return true  # coeffs in zero ring
+  unit_seen = false
+  for i in 0:degree(f.poly)
+    if is_nilpotent(coeff(f.poly, i))
+      continue
+    end
+    if unit_seen || !is_unit(coeff(f.poly, i))
+      return false
+    end
+    unit_seen = true
+  end
+  return unit_seen
 end
 
+
+function is_nilpotent(f::T) where {T <: LaurentPolyRingElem}
+  return is_nilpotent(f.poly);
+end
+
+
 ###############################################################################
 #
 #   Comparisons

src/generic/LaurentMPoly.jl

@@ -109,15 +109,28 @@ function Base.inv(a::LaurentMPolyWrap)
     return LaurentMPolyWrap(parent(a), inv(ap), neg!(ad, ad))
 end
 
-function is_unit(a::LaurentMPolyWrap)
-    (ap, ad) = _normalize(a)
-    if is_domain_type(elem_type(coefficient_ring(a))) || length(ap) <= 1
-        return is_unit(ap)
-    else
-        throw(NotImplementedError(:is_unit, a))
+function is_unit(f::T) where {T <: LaurentMPolyRingElem}
+  # **NOTE** f.mpoly is not normalized in any way
+  is_trivial(parent(f)) && return true  # coeffs in zero ring
+  unit_seen = false
+  for i in 1:length(f.mpoly)
+    if is_nilpotent(coeff(f.mpoly, i))
+      continue
     end
+    if unit_seen || !is_unit(coeff(f.mpoly, i))
+      return false
+    end
+    unit_seen = true
+  end
+  return unit_seen
 end
 
+
+function is_nilpotent(f::T) where {T <: LaurentMPolyRingElem}
+  return is_nilpotent(f.mpoly);
+end
+
+
 is_zero_divisor(p::LaurentMPolyWrap) = is_zero_divisor(p.mpoly)
 
 function is_zero_divisor_with_annihilator(p::LaurentMPolyWrap)

src/generic/LaurentPoly.jl

@@ -201,15 +201,6 @@ function Base.inv(p::LaurentPolyWrap)
    return LaurentPolyWrap(parent(p), inv(g), -p.mindeg-v)
 end
 
-function is_unit(p::LaurentPolyWrap)
-   v, g = _remove_gen(p)
-   if is_domain_type(elem_type(coefficient_ring(p))) || length(g) <= 1
-      return is_unit(g)
-   else
-      throw(NotImplementedError(:is_unit, p))
-   end
-end
-
 is_zero_divisor(p::LaurentPolyWrap) = is_zero_divisor(p.poly)
 
 function is_zero_divisor_with_annihilator(p::LaurentPolyWrap)

src/generic/imports.jl

@@ -149,6 +149,7 @@ import ..AbstractAlgebra: is_exact_type
 import ..AbstractAlgebra: is_finite
 import ..AbstractAlgebra: is_gen
 import ..AbstractAlgebra: is_monomial
+import ..AbstractAlgebra: is_nilpotent
 import ..AbstractAlgebra: is_power
 import ..AbstractAlgebra: is_square
 import ..AbstractAlgebra: is_square_with_sqrt

src/julia/Float.jl

@@ -38,7 +38,7 @@ zero(::Floats{T}) where T <: AbstractFloat = T(0)
 
 one(::Floats{T}) where T <: AbstractFloat = T(1)
 
-is_unit(a::AbstractFloat) = a != 0
+is_unit(a::AbstractFloat) = !is_zero(a)
 
 canonical_unit(a::AbstractFloat) = a
 

src/julia/Rational.jl

@@ -38,7 +38,7 @@ zero(::Rationals{T}) where T <: Integer = Rational{T}(0)
 
 one(::Rationals{T}) where T <: Integer = Rational{T}(1)
 
-is_unit(a::Rational) = a != 0
+is_unit(a::Rational) = !is_zero(a)
 
 is_zero_divisor(a::Rational) = is_zero(a)