Support univariate Laurent divisibility over non-domains

src/sage/rings/polynomial/laurent_polynomial.pyx

@@ -2222,11 +2222,9 @@ cdef class LaurentPolynomial_univariate(LaurentPolynomial):
         r"""
         Return ``True`` if ``self`` divides ``other``.
 
-        .. NOTE::
-
-            This method is only implemented for Laurent polynomials over
-            integral domains. For rings with zero divisors, a
-            :exc:`NotImplementedError` is raised.
+        Over integral domains, this delegates to polynomial divisibility.
+        Otherwise, when supported by the polynomial backend, divisibility is
+        tested as ideal membership in the polynomial cover ring.
 
         EXAMPLES::
 
@@ -2256,16 +2254,67 @@ cdef class LaurentPolynomial_univariate(LaurentPolynomial):
         Check that :issue:`40372` is fixed::
 
             sage: R.<y> = LaurentPolynomialRing(Zmod(4))
-            sage: a = 2+y
+            sage: a = 2 + y
             sage: a.divides(a)
+            True
+            sage: a.divides(2*a)
+            True
+            sage: (2 + y^-1).divides(1)
+            True
+
+        Negative powers and nilpotent coefficients are handled without a
+        degree bound::
+
+            sage: R.<y> = LaurentPolynomialRing(Zmod(8))
+            sage: (4*y^2 + y + 4).divides(2)
+            True
+            sage: R.<y> = LaurentPolynomialRing(Zmod(2^100))
+            sage: (y - 2).divides(y^100)
+            True
+            sage: R.<y> = LaurentPolynomialRing(Zmod(6))
+            sage: (2*y).divides(3*y)
+            False
+
+        Zero and monomial-unit cases do not require a Gröbner basis::
+
+            sage: z = R.zero()
+            sage: z.divides(z), z.divides(y), y.divides(z)
+            (True, False, True)
+            sage: (y^-100000).divides(1)
+            True
+            sage: a = y^-100000 * (2 + y)
+            sage: a.divides(a * (1 + y))
+            True
+            sage: (2 + y).divides((2 + y) * y^-100000)
+            True
+
+        Unsupported polynomial backends raise instead of silently returning
+        an incorrect result::
+
+            sage: A.<e> = QQ[]
+            sage: B.<ebar> = A.quotient(e^2)
+            sage: L.<x> = LaurentPolynomialRing(B)
+            sage: f = 1 + ebar*x
+            sage: f.divides(f * (x^-2 + ebar))
             Traceback (most recent call last):
             ...
-            NotImplementedError: divisibility test not implemented for Laurent polynomials over non-integral domains
+            NotImplementedError: ideal membership in the polynomial cover is not implemented over this base ring
         """
         if self.base_ring().is_integral_domain() is True:
             p = self.polynomial_construction()[0]
             q = other.polynomial_construction()[0]
             return p.divides(q)
-        else:
-            raise NotImplementedError("divisibility test not implemented for Laurent"
-                                      " polynomials over non-integral domains")
+        if not other:
+            return True
+        if not self:
+            return False
+        if self is other or self.is_unit():
+            return True
+        p, n = self.polynomial_construction()
+        q, m = other.polynomial_construction()
+        if n == m and p == q:
+            return True
+        P = self.parent()
+        I = P._ideal_in_extended_ring((self,))
+        target = P._extended_ring(q)
+        return I.reduce(target).is_zero()

src/sage/rings/polynomial/laurent_polynomial_ideal.py

@@ -94,6 +94,10 @@ def __init__(self, ring, gens, coerce=True, hint=None) -> None:
         Ideal_generic.__init__(self, ring, gens, coerce=coerce)
         self._poly_ring = ring.polynomial_ring()
         self._poly_ideal = None  # Create only as needed
+        self._extended_ideal = None  # Create only as needed
+        self._use_polynomial_cover = (
+            isinstance(ring, LaurentPolynomialRing_univariate)
+            and ring.base_ring().is_integral_domain() is not True)
         self._saturated = False
         if hint is None:
             self._hint = self._poly_ring.zero_ideal()
@@ -200,15 +204,76 @@ def __contains__(self, f) -> bool:
             sage: I = P.ideal([x^2 + 3*x])
             sage: 1 + 3*x^-1 in I
             True
+
+        For univariate Laurent polynomials over non-integral base rings,
+        membership is computed in the polynomial cover ring when supported::
+
+            sage: P.<x> = LaurentPolynomialRing(Zmod(8))
+            sage: 2 in P.ideal(4*x^2 + x + 4)
+            True
+            sage: 1 in P.ideal(x^-1)
+            True
+            sage: 2 + x in P.ideal(2*x^-1 + 1)
+            True
+            sage: 2 + x in P.ideal((2 + x)*x^-100000)
+            True
+            sage: 2*x^-3 in P.ideal(2*x)
+            True
+            sage: 3*x in P.ideal(2*x)
+            False
+            sage: 1 in P.ideal(x^-100000)
+            True
+
+        The zero ideal and zero ring are handled without a Gröbner basis::
+
+            sage: I = P.ideal([])
+            sage: 0 in I, 1 in I
+            (True, False)
+            sage: Z.<z> = LaurentPolynomialRing(Zmod(1))
+            sage: Z.zero() in Z.ideal([])
+            True
+
+        Unsupported polynomial backends raise instead of silently returning
+        an incorrect result::
+
+            sage: A.<e> = QQ[]
+            sage: B.<ebar> = A.quotient(e^2)
+            sage: L.<x> = LaurentPolynomialRing(B)
+            sage: f = 1 + ebar*x
+            sage: f * (x^-2 + ebar) in L.ideal(f)
+            Traceback (most recent call last):
+            ...
+            NotImplementedError: ideal membership in the polynomial cover is not implemented over this base ring
         """
-        if not f or f in self.gens():
+        R = self.ring()
+        if not self._use_polynomial_cover:
+            if not f or f in self.gens():
+                return True
+            f = R(f)
+            if isinstance(R, LaurentPolynomialRing_univariate):
+                g = f.__reduce__()[1][1]
+            else:
+                g = f.__reduce__()[1][0]
+            return g in self.polynomial_ideal()
+
+        f = R(f)
+        generators = self.gens()
+        if not f or any(f is g for g in generators):
             return True
-        f = self.ring()(f)
-        if isinstance(self.ring(), LaurentPolynomialRing_univariate):
-            g = f.__reduce__()[1][1]
-        else:
-            g = f.__reduce__()[1][0]
-        return (g in self.polynomial_ideal())
+        f_polynomial, f_shift = f.polynomial_construction()
+        for g in generators:
+            g_polynomial, g_shift = g.polynomial_construction()
+            if f_shift == g_shift and f_polynomial == g_polynomial:
+                return True
+        if not any(generators):
+            return False
+        if self._extended_ideal is None:
+            if any(g.is_unit() for g in generators):
+                self._extended_ideal = R._extended_ring.unit_ideal()
+            else:
+                self._extended_ideal = R._ideal_in_extended_ring(generators)
+        target = R._extended_ring(f_polynomial)
+        return self._extended_ideal.reduce(target).is_zero()
 
     def gens_reduced(self) -> tuple:
         """

src/sage/rings/polynomial/laurent_polynomial_ring_base.py

@@ -71,6 +71,27 @@ def __init__(self, R) -> None:
         self._extended_ring_ideal = ER.ideal([ER.gen(2*i) * ER.gen(2*i+1) - 1
                                               for i in range(self._n)])
 
+    def _ideal_in_extended_ring(self, generators):
+        """
+        Return the ideal generated by univariate ``generators`` in the
+        polynomial cover.
+
+        The inversion relations defining this Laurent polynomial ring are
+        included among the generators.  Laurent monomial factors, which are
+        units, are removed to keep exponents in the cover as small as possible.
+        """
+        from sage.rings.polynomial.multi_polynomial_libsingular import (
+            MPolynomialRing_libsingular,
+        )
+
+        if not isinstance(self._extended_ring, MPolynomialRing_libsingular):
+            raise NotImplementedError("ideal membership in the polynomial cover "
+                                      "is not implemented over this base ring")
+        generators = tuple(self._extended_ring(g.polynomial_construction()[0])
+                           for g in generators)
+        return self._extended_ring.ideal(
+            generators + tuple(self._extended_ring_ideal.gens()))
+
     def ngens(self):
         """
         Return the number of generators of ``self``.