gh-39555: Added option to return explicit polynomial in .in_subalgebra
    
<!-- ^ Please provide a concise and informative title. -->
<!-- ^ Don't put issue numbers in the title, do this in the PR
description below. -->
<!-- ^ For example, instead of "Fixes #12345" use "Introduce new method
to calculate 1 + 2". -->
<!-- v Describe your changes below in detail. -->
<!-- v Why is this change required? What problem does it solve? -->
<!-- v If this PR resolves an open issue, please link to it here. For
example, "Fixes #12345". -->

Fixes #39225.

Previously the function in_subalgebra in
src/sage/rins/polynomial/multi_polynomial_libsingular.pyx returned True
or False, depending on whether the polynomial was present in the
provided subalgebra. When the parameter 'algorithm' was set to
'groebner' it would do this by generating a polynomial, but this
polynomial would not be made available to the user.

This function was enhanced by adding a key-word only argument
'certificate' to the function which defaults to False. When the value
provided for this argument evaluates to False the function behaves as
before. When the value provided evaluates to True and the Groebner
algorithm is used the function returns the polynomial generated by the
algorithm when it would have otherwise returned True, and returns None
when it would have otherwise returned False.

Additionally, if the value provided for the 'certificate' argument is a
string which evaluates to True then that string will be used to name the
variables in the generated polynomial.

The existing example for this function has been modified to test this
new functionality.

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

- [X] The title is concise and informative.
- [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 and checked the documentation
preview.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on. For example,
-->
<!-- - #12345: short description why this is a dependency -->
<!-- - #34567: ... -->
    
URL: #39555
Reported by: Caleb Van't Land
Reviewer(s): Caleb Van't Land, Travis Scrimshaw

Author: Release Manager

src/sage/rings/polynomial/multi_polynomial_libsingular.pyx

@@ -6073,7 +6073,7 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
         else:
             return self * self.denominator()
 
-    def in_subalgebra(self, J, algorithm='algebra_containment'):
+    def in_subalgebra(self, J, algorithm='algebra_containment', *, certificate=False):
         """
         Return whether this polynomial is contained in the subalgebra
         generated by ``J``
@@ -6104,6 +6104,12 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
             and only if the remainder involves only the new variables
             `z_i`.
 
+        - ``certificate`` --  boolean (default: ``False``) or string; if``True``
+          and ``algorithm='groebner'``, return the polynomial generated by
+          the algorithm if such a polynomial exists, otherwise return ``None``;
+          if a nonempty string, then this is used as the name of the variables
+          in the returned polynomial, otherwise the name will be ``'newgens'``
+
         EXAMPLES::
 
             sage: P.<x,y,z> = QQ[]
@@ -6115,6 +6121,8 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
             True
             sage: (a^2).in_subalgebra(J, algorithm='groebner')
             True
+            sage: (a^2).in_subalgebra(J, algorithm='groebner', certificate='x')
+            x0^2*x1^2
             sage: (a + a^2).in_subalgebra(J)
             True
         """
@@ -6135,14 +6143,22 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
         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))]
+            var_name = "newgens"
+            if (certificate and isinstance(certificate, str)):
+                var_name = certificate
+            new_gens = [f"{var_name}{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)
+            if certificate:
+                if set(z.variables()).issubset(new_gens):
+                    T = PolynomialRing(R.base_ring(), tuple(new_gens))
+                    return T(z)
+                return None
             return set(z.variables()).issubset(new_gens)
         else:
             raise ValueError("unknown algorithm")