Implement algebra_containment from Singular (issue #34502)

src/sage/rings/polynomial/multi_polynomial_libsingular.pyx

@@ -5744,6 +5744,78 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
         else:
             return self * self.denominator()
 
+    def in_subalgebra(self, J, algorithm=None):
+        """
+        Return whether this polynomial is contained in the subalgebra
+        generated by ``J``
+
+        INPUT:
+
+        - ``J`` -- iterable of elements of the parent polynomial ring
+
+        - ``algorithm`` -- can be ``None`` (the default), "algebra_containment",
+          "inSubring", or "groebner".
+
+          - ``None`` (the default): use "algebra_containment".
+
+          - "algebra_containment": use Singular's ``algebra_containment`` function,
+            https://www.singular.uni-kl.de/Manual/4-2-1/sing_1247.htm#SEC1328. The
+            Singular documentation suggests that this is frequently faster than the
+            next option.
+
+          - "inSubring": use Singular's ``inSubring`` function,
+            https://www.singular.uni-kl.de/Manual/4-2-0/sing_1240.htm#SEC1321.
+
+          - "groebner": use the algorithm described in Singular's
+            documentation, but within Sage: if the subalgebra
+            generators are `y_1`, ..., `y_m`, then create a new
+            polynomial algebra with the old generators along with new
+            ones: `z_1`, ..., `z_m`. Create the ideal `(z_1 - y_1,
+            ..., z_m - y_m)`, and reduce the polynomial modulo this
+            ideal. The polynomial is contained in the subalgebra if
+            and only if the remainder involves only the new variables
+            `z_i`.
+
+        Note: the default algorithm may print some intermediate information.
+
+        EXAMPLES::
+
+            sage: P.<x,y,z> = QQ[]
+            sage: J = [x^2 + y^2, x^2 + z^2]
+            sage: (y^2).in_subalgebra(J)
+            ...
+            False
+            sage: a = (x^2 + y^2) * (x^2 + z^2)
+            sage: a.in_subalgebra(J, algorithm='inSubring')
+            True
+            sage: (a^2).in_subalgebra(J, algorithm='groebner')
+            True
+            sage: (a + a^2).in_subalgebra(J)
+            ...
+            True
+        """
+        R = self.parent()
+        if algorithm is not None:
+            algorithm = algorithm.lower()
+        from sage.libs.singular.function import singular_function, lib as singular_lib
+        singular_lib('algebra.lib')
+        if algorithm is None or algorithm == "algebra_containment":
+            return singular_function('algebra_containment')(self, R.ideal(J)) == 1
+        elif algorithm == "insubring":
+            return singular_function('inSubring')(self, R.ideal(J))[0] == 1
+        elif algorithm == "groebner":
+            new_gens = [f"newgens{i}" for i in range(len(J))]
+            ngens = len(new_gens)
+            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+            S = PolynomialRing(R.base_ring(), R.gens() + tuple(new_gens),
+                               order=f'degrevlex({len(R.gens())}),degrevlex({ngens})')
+            new_gens = S.gens()[-ngens:]
+            I = S.ideal([g - S(j) for (g,j) in zip(new_gens, J)])
+            z = S(self).reduce(I)
+            return set(z.variables()).issubset(new_gens)
+        else:
+            raise ValueError("unknown algorithm")
+
 
 def unpickle_MPolynomial_libsingular(MPolynomialRing_libsingular R, d):
     """