gh-36030: Implement algebra_containment from Singular (issue #34502)
    
<!-- ^^^^^
Please provide a concise, informative and self-explanatory title.
Don't put issue numbers in there, do this in the PR body below.
For example, instead of "Fixes #1234" use "Introduce new method to
calculate 1+1"
-->

This implements subalgebra containment for polynomials: whether a
polynomial is contained in the subalgebra determined by given
generators.

<!-- Describe your changes here in detail -->

<!-- Why is this change required? What problem does it solve? -->
<!-- If this PR resolves an open issue, please link to it here. For
example "Fixes #12345". -->
<!-- If your change requires a documentation PR, please link it
appropriately. -->

This fixes #34502. It provides several algorithms for the calculution,
using Singular's algebra_containment, Singular's inSubring, or Sage's
Groebner basis tools.


### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->
<!-- If your change requires a documentation PR, please link it
appropriately -->
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
<!-- Feel free to remove irrelevant items. -->

- [X] The title is concise, informative, and self-explanatory.
- [X] The description explains in detail what this PR is about.
- [X] I have linked a relevant issue or discussion.
- [X] I have created tests covering the changes.
- [X] I have updated the documentation accordingly.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on
- #12345: short description why this is a dependency
- #34567: ...
-->

<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
    
URL: #36030
Reported by: John H. Palmieri
Reviewer(s): John H. Palmieri, Kwankyu Lee

Author: Release Manager

src/sage/rings/polynomial/multi_polynomial_libsingular.pyx

@@ -5757,6 +5757,83 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
         else:
             return self * self.denominator()
 
+    def in_subalgebra(self, J, algorithm="algebra_containment"):
+        """
+        Return whether this polynomial is contained in the subalgebra
+        generated by ``J``
+
+        INPUT:
+
+        - ``J`` -- list of elements of the parent polynomial ring
+
+        - ``algorithm`` -- can be ``"algebra_containment"`` (the default),
+          ``"inSubring"``, or ``"groebner"``.
+
+          - ``"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`.
+
+        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()
+        algorithm = algorithm.lower()
+        from sage.libs.singular.function import singular_function, lib as singular_lib
+        singular_lib('algebra.lib')
+        if algorithm == "algebra_containment":
+            execute = singular_function('execute')
+            try:
+                get_printlevel = singular_function('get_printlevel')
+            except NameError:
+                execute('proc get_printlevel {return (printlevel);}')
+                get_printlevel = singular_function('get_printlevel')
+            # It's fairly verbose unless printlevel is -1.
+            saved_printlevel = get_printlevel()
+            execute('printlevel=-1')
+            contains = singular_function('algebra_containment')(self, R.ideal(J)) == 1
+            execute('printlevel={}'.format(saved_printlevel))
+            return contains
+        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):
     """