gh-39889: Implement algorithm parameter for .series()
    
See #6119

Instead of deprecating it, this just add algorithm=maxima. Which is less
intrusive.

### 📝 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: #39889
Reported by: user202729
Reviewer(s): Travis Scrimshaw, user202729

Author: Release Manager

src/sage/symbolic/expression.pyx

@@ -4757,7 +4757,7 @@ cdef class Expression(Expression_abc):
         return matrix([[g.derivative(x) for x in self.arguments()]
                        for g in self.gradient()])
 
-    def series(self, symbol, order=None):
+    def series(self, symbol, order=None, algorithm='ginac'):
         r"""
         Return the power series expansion of ``self`` in terms of the
         given variable to the given order.
@@ -4771,6 +4771,10 @@ cdef class Expression(Expression_abc):
         - ``order`` -- integer; if nothing given, it is set
           to the global default (``20``), which can be changed
           using :func:`set_series_precision`
+        - ``algorithm`` -- string (default: ``'ginac'``); one of the following:
+
+          * ``'ginac'``
+          * ``'maxima'``
 
         OUTPUT: a power series
 
@@ -4869,7 +4873,20 @@ cdef class Expression(Expression_abc):
 
             sage: ((1 - x)^-x).series(x, 8)
             1 + 1*x^2 + 1/2*x^3 + 5/6*x^4 + 3/4*x^5 + 33/40*x^6 + 5/6*x^7 + Order(x^8)
+
+        Try different algorithms::
+
+            sage: ((1 - x)^-x).series(x, 8, algorithm="maxima")
+            1 + 1*x^2 + 1/2*x^3 + 5/6*x^4 + 3/4*x^5 + 33/40*x^6 + 5/6*x^7 + Order(x^8)
+            sage: ((1 - x)^-x).series(x, 8, algorithm="ginac")
+            1 + 1*x^2 + 1/2*x^3 + 5/6*x^4 + 3/4*x^5 + 33/40*x^6 + 5/6*x^7 + Order(x^8)
         """
+        if algorithm == "maxima":
+            # call series() again to convert the result (a rational function in the symbol)
+            # to a SymbolicSeries with the correct order
+            return self.taylor(symbol, 0, order-1).series(symbol, order, algorithm="ginac")
+        if algorithm != "ginac":
+            raise ValueError("invalid algorithm")
         cdef Expression symbol0 = self.coerce_in(symbol)
         cdef GEx x
         cdef SymbolicSeries nex
@@ -4981,6 +4998,10 @@ cdef class Expression(Expression_abc):
 
            - ``(x, a)``, ``(y, b)``, ``n`` -- variables with points, degree of polynomial
 
+        .. SEEALSO::
+
+            :meth:`series`
+
         EXAMPLES::
 
             sage: var('a, x, z')