17402: merge 17659; add SymbolicSeries.power_series

src/doc/en/reference/calculus/index.rst

@@ -13,6 +13,7 @@ Symbolic Calculus
    sage/symbolic/function
    sage/symbolic/function_factory
    sage/calculus/functional
+   sage/symbolic/series
    sage/symbolic/integration/integral
    sage/calculus/test_sympy
    sage/calculus/tests

src/sage/symbolic/expression.pyx

@@ -142,6 +142,7 @@ import sage.rings.integer
 import sage.rings.rational
 from sage.structure.element cimport ModuleElement, RingElement, Element
 from sage.symbolic.getitem cimport OperandsWrapper
+from sage.symbolic.series cimport SymbolicSeries
 from sage.symbolic.complexity_measures import string_length
 from sage.symbolic.function import get_sfunction_from_serial, SymbolicFunction
 from sage.rings.rational import Rational  # Used for sqrt.
@@ -1772,12 +1773,7 @@ cdef class Expression(CommutativeRingElement):
 
     def is_series(self):
         """
-        Return True if ``self`` is a series.
-
-        Series are special kinds of symbolic expressions that are
-        constructed via the :meth:`series` method. They usually have
-        an ``Order()`` term unless the series representation is exact,
-        see :meth:`is_terminating_series`.
+        Return True if ``self`` is a :class:`~sage.symbolic.series.SymbolicSeries`.
 
         OUTPUT:
 
@@ -1817,7 +1813,7 @@ cdef class Expression(CommutativeRingElement):
             sage: sum_expr.is_series()
             False
         """
-        return is_a_series(self._gobj)
+        return False
 
     def is_terminating_series(self):
         """
@@ -1849,7 +1845,7 @@ cdef class Expression(CommutativeRingElement):
             sage: exp(x).series(x,10).is_terminating_series()
             False
         """
-        return g_is_a_terminating_series(self._gobj)
+        return False
 
     cpdef bint is_polynomial(self, var):
         """
@@ -3595,6 +3591,7 @@ cdef class Expression(CommutativeRingElement):
         """
         cdef Expression symbol0 = self.coerce_in(symbol)
         cdef GEx x
+        cdef SymbolicSeries nex
         cdef int prec
         if order is None:
             from sage.misc.defaults import series_precision
@@ -3604,9 +3601,12 @@ cdef class Expression(CommutativeRingElement):
         sig_on()
         try:
             x = self._gobj.series(symbol0._gobj, prec, 0)
+            nex = <SymbolicSeries>PY_NEW(SymbolicSeries)
+            nex._parent = self._parent
+            GEx_construct_ex(&nex._gobj, x)
         finally:
             sig_off()
-        return new_Expression_from_GEx(self._parent, x)
+        return nex
 
     def residue(self, symbol):
         """
@@ -3766,9 +3766,7 @@ cdef class Expression(CommutativeRingElement):
             sage: f.series(x==1,3).truncate().expand()
             -2*x^2*cos(1) + 5/2*x^2*sin(1) + 5*x*cos(1) - 7*x*sin(1) - 3*cos(1) + 11/2*sin(1)
         """
-        if not is_a_series(self._gobj):
-            return self
-        return new_Expression_from_GEx(self._parent, series_to_poly(self._gobj))
+        return self
 
     def expand(Expression self, side=None):
         """
@@ -5239,34 +5237,17 @@ cdef class Expression(CommutativeRingElement):
             doctest:...: DeprecationWarning: coeffs is deprecated. Please use coefficients instead.
             See http://trac.sagemath.org/17438 for details.
             [[1, 1]]
-
-        Series coefficients are now handled correctly (:trac:`17399`)::
-
-            sage: s=(1/(1-x)).series(x,6); s
-            1 + 1*x + 1*x^2 + 1*x^3 + 1*x^4 + 1*x^5 + Order(x^6)
-            sage: s.coefficients()
-            [[1, 0], [1, 1], [1, 2], [1, 3], [1, 4], [1, 5]]
-            sage: s.coefficients(x, sparse=False)
-            [1, 1, 1, 1, 1, 1]
-            sage: x,y = var("x,y")
-            sage: s=(1/(1-y*x-x)).series(x,3); s
-            1 + (y + 1)*x + ((y + 1)^2)*x^2 + Order(x^3)
-            sage: s.coefficients(x, sparse=False)
-            [1, y + 1, (y + 1)^2]
         """
+        f = self._maxima_()
+        maxima = f.parent()
+        maxima._eval_line('load(coeflist)')
         if x is None:
             x = self.default_variable()
-        if is_a_series(self._gobj):
-            l = [[self.coefficient(x, d), d] for d in xrange(self.degree(x))]
-        else:
-            f = self._maxima_()
-            maxima = f.parent()
-            maxima._eval_line('load(coeflist)')
-            G = f.coeffs(x)
-            from sage.calculus.calculus import symbolic_expression_from_maxima_string
-            S = symbolic_expression_from_maxima_string(repr(G))
-            l = S[1:]
-
+        x = self.parent().var(repr(x))
+        G = f.coeffs(x)
+        from sage.calculus.calculus import symbolic_expression_from_maxima_string
+        S = symbolic_expression_from_maxima_string(repr(G))
+        l = S[1:]
         if sparse is True:
             return l
         else:
@@ -5856,11 +5837,11 @@ cdef class Expression(CommutativeRingElement):
             sage: g.parent()
             Power Series Ring in theta over Rational Field
         """
+        from sage.rings.all import PowerSeriesRing
         v = self.variables()
         if len(v) != 1:
             raise ValueError("self must be a polynomial in one variable but it is in the variables %s" % tuple([v]))
         f = self.polynomial(base_ring)
-        from sage.rings.all import PowerSeriesRing
         R = PowerSeriesRing(base_ring, names=f.parent().variable_names())
         return R(f, f.degree()+1)
 
@@ -10968,3 +10949,4 @@ cdef operators compatible_relation(operators lop, operators rop) except <operato
        return greater
     else:
         raise TypeError("incompatible relations")
+

src/sage/symbolic/ginac.pxd

@@ -134,6 +134,7 @@ cdef extern from "ginac_wrap.h":
     bint is_a_series "is_a<pseries>" (GEx e)
     # you must ensure that is_a_series(e) is true before calling this:
     bint g_is_a_terminating_series(GEx e) except +
+    GEx g_series_var(GEx e) except +
 
     # Relations
     ctypedef enum operators "relational::operators":

src/sage/symbolic/ginac_wrap.h

@@ -76,6 +76,13 @@ bool g_is_a_terminating_series(const ex& e) {
     return false;
 }
 
+ex g_series_var(const ex& e) {
+    if (is_a<pseries>(e)) {
+        return (ex_to<pseries>(e)).get_var();
+    }
+    return NULL;
+}
+
 relational::operators relational_operator(const ex& e) {
     // unsafe cast -- be damn sure the input is a relational.
     return (ex_to<relational>(e)).the_operator();

src/sage/symbolic/series.pxd

@@ -0,0 +1,4 @@
+from sage.symbolic.expression cimport Expression
+
+cdef class SymbolicSeries(Expression):
+    pass

src/sage/symbolic/series.pyx

@@ -0,0 +1,229 @@
+"""
+Symbolic Series
+
+Symbolic series are special kinds of symbolic expressions that are
+constructed via the
+:meth:`Expression.series <sage.symbolic.expression.Expression.series>`
+method.
+They usually have an ``Order()`` term unless the series representation
+is exact, see
+:meth:`~sage.symbolic.series.SymbolicSeries.is_terminating_series`.
+
+For series over general rings see
+:class:`power series <sage.rings.power_series_poly.PowerSeries_poly>`
+and
+:class:`Laurent series<sage.rings.laurent_series_ring_element.LaurentSeries>`.
+
+EXAMPLES:
+
+We expand a polynomial in `x` about 0, about `1`, and also truncate
+it back to a polynomial::
+
+    sage: var('x,y')
+    (x, y)
+    sage: f = (x^3 - sin(y)*x^2 - 5*x + 3); f
+    x^3 - x^2*sin(y) - 5*x + 3
+    sage: g = f.series(x, 4); g
+    3 + (-5)*x + (-sin(y))*x^2 + 1*x^3
+    sage: g.truncate()
+    x^3 - x^2*sin(y) - 5*x + 3
+    sage: g = f.series(x==1, 4); g
+    (-sin(y) - 1) + (-2*sin(y) - 2)*(x - 1) + (-sin(y) + 3)*(x - 1)^2 + 1*(x - 1)^3
+    sage: h = g.truncate(); h
+    (x - 1)^3 - (x - 1)^2*(sin(y) - 3) - 2*(x - 1)*(sin(y) + 1) - sin(y) - 1
+    sage: h.expand()
+    x^3 - x^2*sin(y) - 5*x + 3
+
+We compute another series expansion of an analytic function::
+
+    sage: f = sin(x)/x^2
+    sage: f.series(x,7)
+    1*x^(-1) + (-1/6)*x + 1/120*x^3 + (-1/5040)*x^5 + Order(x^7)
+    sage: f.series(x==1,3)
+    (sin(1)) + (cos(1) - 2*sin(1))*(x - 1) + (-2*cos(1) + 5/2*sin(1))*(x - 1)^2 + Order((x - 1)^3)
+    sage: f.series(x==1,3).truncate().expand()
+    -2*x^2*cos(1) + 5/2*x^2*sin(1) + 5*x*cos(1) - 7*x*sin(1) - 3*cos(1) + 11/2*sin(1)
+
+Following the GiNaC tutorial, we use John Machin's amazing
+formula `\pi = 16 \mathrm{tan}^{-1}(1/5) - 4 \mathrm{tan}^{-1}(1/239)`
+to compute digits of `\pi`. We expand the arc tangent around 0 and insert
+the fractions 1/5 and 1/239.
+
+::
+
+    sage: x = var('x')
+    sage: f = atan(x).series(x, 10); f
+    1*x + (-1/3)*x^3 + 1/5*x^5 + (-1/7)*x^7 + 1/9*x^9 + Order(x^10)
+    sage: float(16*f.subs(x==1/5) - 4*f.subs(x==1/239))
+    3.1415926824043994
+"""
+########################################################################
+#       Copyright (C) 2015 Ralf Stephan <ralf@ark.in-berlin.de>
+#
+#  Distributed under the terms of the GNU General Public License (GPL)
+#   as published by the Free Software Foundation; either version 2 of
+#   the License, or (at your option) any later version.
+#                   http://www.gnu.org/licenses/
+########################################################################
+
+include "sage/ext/interrupt.pxi"
+include "sage/ext/stdsage.pxi"
+include "sage/ext/cdefs.pxi"
+include "sage/ext/python.pxi"
+
+from ginac cimport *
+from sage.symbolic.expression cimport Expression, new_Expression_from_GEx
+
+cdef class SymbolicSeries(Expression):
+    def __init__(self, SR):
+        Expression.__init__(self, SR, 0)
+        self._parent = SR
+
+    def is_series(self):
+        """
+        Return True.
+
+        EXAMPLES::
+
+            sage: exp(x).series(x,10).is_series()
+            True
+        """
+        return True
+
+    def is_terminating_series(self):
+        """
+        Return True if ``self`` is without order term.
+
+        A series is terminating if it can be represented exactly,
+        without requiring an order term.
+
+        OUTPUT:
+
+        Boolean. Whether ``self`` has no order term.
+
+        EXAMPLES::
+
+            sage: (x^5+x^2+1).series(x,10)
+            1 + 1*x^2 + 1*x^5
+            sage: (x^5+x^2+1).series(x,10).is_series()
+            True
+            sage: (x^5+x^2+1).series(x,10).is_terminating_series()
+            True
+            sage: SR(5).is_terminating_series()
+            False
+            sage: exp(x).series(x,10).is_terminating_series()
+            False
+        """
+        return g_is_a_terminating_series((<Expression>self)._gobj)
+
+    def truncate(self):
+        """
+        Given a power series or expression, return the corresponding
+        expression without the big oh.
+
+        INPUT:
+
+        - ``self`` -- a series as output by the :meth:`series` command.
+
+        OUTPUT:
+
+        A symbolic expression.
+
+        EXAMPLES::
+
+            sage: f = sin(x)/x^2
+            sage: f.truncate()
+            sin(x)/x^2
+            sage: f.series(x,7)
+            1*x^(-1) + (-1/6)*x + 1/120*x^3 + (-1/5040)*x^5 + Order(x^7)
+            sage: f.series(x,7).truncate()
+            -1/5040*x^5 + 1/120*x^3 - 1/6*x + 1/x
+            sage: f.series(x==1,3).truncate().expand()
+            -2*x^2*cos(1) + 5/2*x^2*sin(1) + 5*x*cos(1) - 7*x*sin(1) - 3*cos(1) + 11/2*sin(1)
+        """
+        return new_Expression_from_GEx(self._parent, series_to_poly(self._gobj))
+
+    def default_variable(self):
+        """
+        Return the expansion variable of this symbolic series.
+
+        EXAMPLES::
+
+            sage: s=(1/(1-x)).series(x,3); s
+            1 + 1*x + 1*x^2 + Order(x^3)
+            sage: s.default_variable()
+            x
+        """
+        cdef GEx x = g_series_var(self._gobj)
+        cdef Expression ex = new_Expression_from_GEx(self._parent, x)
+        return ex
+
+    def coefficients(self, x=None, sparse=True):
+        r"""
+        Return the coefficients of this symbolic series as a polynomial in x.
+
+        INPUT:
+
+        -  ``x`` -- optional variable.
+
+        OUTPUT:
+
+        Depending on the value of ``sparse``,
+
+        - A list of pairs ``(expr, n)``, where ``expr`` is a symbolic
+          expression and ``n`` is a power (``sparse=True``, default)
+
+        - A list of expressions where the ``n``-th element is the coefficient of
+          ``x^n`` when self is seen as polynomial in ``x`` (``sparse=False``).
+
+        EXAMPLES::
+
+            sage: s=(1/(1-x)).series(x,6); s
+            1 + 1*x + 1*x^2 + 1*x^3 + 1*x^4 + 1*x^5 + Order(x^6)
+            sage: s.coefficients()
+            [[1, 0], [1, 1], [1, 2], [1, 3], [1, 4], [1, 5]]
+            sage: s.coefficients(x, sparse=False)
+            [1, 1, 1, 1, 1, 1]
+            sage: x,y = var("x,y")
+            sage: s=(1/(1-y*x-x)).series(x,3); s
+            1 + (y + 1)*x + ((y + 1)^2)*x^2 + Order(x^3)
+            sage: s.coefficients(x, sparse=False)
+            [1, y + 1, (y + 1)^2]
+
+        """
+        if x is None:
+            x = self.default_variable()
+        l = [[self.coefficient(x, d), d] for d in xrange(self.degree(x))]
+        if sparse is True:
+            return l
+        else:
+            from sage.rings.integer_ring import ZZ
+            if any(not c[1] in ZZ for c in l):
+                raise ValueError("Cannot return dense coefficient list with noninteger exponents.")
+            val = l[0][1]
+            if val < 0:
+                raise ValueError("Cannot return dense coefficient list with negative valuation.")
+            deg = l[-1][1]
+            ret = [ZZ(0)] * int(deg+1)
+            for c in l:
+                ret[c[1]] = c[0]
+            return ret
+
+    def power_series(self, base_ring):
+        """
+        Return algebraic power series associated to this symbolic
+        series. The coefficients must be coercible to the base ring.
+
+        EXAMPLES::
+
+            sage: ex=(gamma(1-x)).series(x,3); ex
+            1 + (euler_gamma)*x + (1/2*euler_gamma^2 + 1/12*pi^2)*x^2 + Order(x^3)
+            sage: g=ex.power_series(SR); g
+            1 + euler_gamma*x + (1/2*euler_gamma^2 + 1/12*pi^2)*x^2 + O(x^3)
+            sage: g.parent()
+            Power Series Ring in x over Symbolic Ring
+        """
+        from sage.rings.all import PowerSeriesRing
+        R = PowerSeriesRing(base_ring, names=str(self.default_variable()))
+        return R(self.list(), self.degree(self.default_variable()))
+