From 3ece975a59933edd7e742dbdc665ab85e1f9dafb Mon Sep 17 00:00:00 2001 From: Ralf Stephan Date: Fri, 23 Jan 2015 08:47:14 +0100 Subject: [PATCH] 16203: conversion from SR.series to PowerSeries --- src/doc/en/reference/calculus/index.rst | 1 + src/sage/rings/power_series_ring.py | 20 +++ src/sage/symbolic/expression.pyx | 31 ++-- src/sage/symbolic/ginac.pxd | 1 + src/sage/symbolic/ginac_wrap.h | 7 + src/sage/symbolic/ring.pyx | 4 +- src/sage/symbolic/series.pxd | 4 + src/sage/symbolic/series.pyx | 210 ++++++++++++++++++++++++ 8 files changed, 264 insertions(+), 14 deletions(-) create mode 100644 src/sage/symbolic/series.pxd create mode 100644 src/sage/symbolic/series.pyx diff --git a/src/doc/en/reference/calculus/index.rst b/src/doc/en/reference/calculus/index.rst index ecb2c931748..5ddeef1bce2 100644 --- a/src/doc/en/reference/calculus/index.rst +++ b/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 diff --git a/src/sage/rings/power_series_ring.py b/src/sage/rings/power_series_ring.py index 9fce0a68972..5ecf248045e 100644 --- a/src/sage/rings/power_series_ring.py +++ b/src/sage/rings/power_series_ring.py @@ -694,6 +694,19 @@ def _element_constructor_(self, f, prec=infinity, check=True): 1 - x + x^2 - x^3 + x^4 + O(x^5) sage: PowerSeriesRing(PowerSeriesRing(QQ,'x'),'x')(x).coefficients() [x] + + Conversion from symbolic series:: + + sage: x,y = var('x,y') + sage: s=(1/(1-x)).series(x,3); s + 1 + 1*x + 1*x^2 + Order(x^3) + sage: R. = PowerSeriesRing(QQ) + sage: R(s) + 1 + x + x^2 + O(x^3) + sage: ex=(gamma(1-y)).series(y,3) + sage: R. = PowerSeriesRing(SR) + sage: R(ex) + 1 + euler_gamma*y + (1/2*euler_gamma^2 + 1/12*pi^2)*y^2 + O(y^3) Laurent series with non-negative valuation are accepted (see :trac:`6431`):: @@ -708,6 +721,7 @@ def _element_constructor_(self, f, prec=infinity, check=True): ArithmeticError: self is a not a power series """ + from sage.symbolic.series import SymbolicSeries if isinstance(f, power_series_ring_element.PowerSeries) and f.parent() is self: if prec >= f.prec(): return f @@ -724,6 +738,12 @@ def _element_constructor_(self, f, prec=infinity, check=True): num = self.element_class(self, f.numerator(), prec, check=check) den = self.element_class(self, f.denominator(), prec, check=check) return self.coerce(num/den) + elif isinstance(f, SymbolicSeries): + if str(f.series_variable()) is self.variable_name(): + return self.element_class(self, f.list(), + f.degree(f.series_variable()), check=check) + else: + raise TypeError("Can only convert series into ring with same variable name.") return self.element_class(self, f, prec, check=check) def construction(self): diff --git a/src/sage/symbolic/expression.pyx b/src/sage/symbolic/expression.pyx index ee7ab44d604..261c650ee07 100644 --- a/src/sage/symbolic/expression.pyx +++ b/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. @@ -1757,12 +1758,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: @@ -1802,7 +1798,7 @@ cdef class Expression(CommutativeRingElement): sage: sum_expr.is_series() False """ - return is_a_series(self._gobj) + return False def is_terminating_series(self): """ @@ -1834,7 +1830,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): """ @@ -3564,13 +3560,17 @@ cdef class Expression(CommutativeRingElement): """ cdef Expression symbol0 = self.coerce_in(symbol) cdef GEx x + cdef SymbolicSeries nex sig_on() try: x = self._gobj.series(symbol0._gobj, order, 0) + nex = 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): """ Calculate the residue of ``self`` with respect to ``symbol``. @@ -3729,9 +3729,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): """ @@ -5183,6 +5181,8 @@ cdef class Expression(CommutativeRingElement): [[2*a^2 + 1, 0], [-2*sqrt(2)*a + 1, 1], [1, 2]] sage: p.coefficients(x, sparse=False) [2*a^2 + 1, -2*sqrt(2)*a + 1, 1] + + TESTS: The behaviour is undefined with noninteger or negative exponents:: @@ -5257,6 +5257,10 @@ cdef class Expression(CommutativeRingElement): -2*sqrt(2)*a*x + 2*a^2 + x^2 + x + 1 sage: p.list(a) [x^2 + x + 1, -2*sqrt(2)*x, 2] + 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.list() + [1, 1, 1, 1, 1, 1] """ return self.coefficients(x=x, sparse=False) @@ -10888,3 +10892,4 @@ cdef operators compatible_relation(operators lop, operators rop) except " (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": diff --git a/src/sage/symbolic/ginac_wrap.h b/src/sage/symbolic/ginac_wrap.h index ee8f89ea450..520a551c628 100644 --- a/src/sage/symbolic/ginac_wrap.h +++ b/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(e)) { + return (ex_to(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(e)).the_operator(); diff --git a/src/sage/symbolic/ring.pyx b/src/sage/symbolic/ring.pyx index 13b9318afce..467546015a3 100644 --- a/src/sage/symbolic/ring.pyx +++ b/src/sage/symbolic/ring.pyx @@ -33,7 +33,9 @@ from sage.structure.parent_base import ParentWithBase from sage.rings.ring cimport CommutativeRing from sage.categories.morphism cimport Morphism -from sage.rings.all import RR, CC +from sage.rings.real_mpfr import RR +from sage.rings.complex_field import ComplexField +CC = ComplexField() pynac_symbol_registry = {} diff --git a/src/sage/symbolic/series.pxd b/src/sage/symbolic/series.pxd new file mode 100644 index 00000000000..ef15506ad16 --- /dev/null +++ b/src/sage/symbolic/series.pxd @@ -0,0 +1,4 @@ +from sage.symbolic.expression cimport Expression + +cdef class SymbolicSeries(Expression): + pass \ No newline at end of file diff --git a/src/sage/symbolic/series.pyx b/src/sage/symbolic/series.pyx new file mode 100644 index 00000000000..d3b16ef0fba --- /dev/null +++ b/src/sage/symbolic/series.pyx @@ -0,0 +1,210 @@ +""" +Symbolic Series + +Symbolic series are special kinds of symbolic expressions that are +constructed via the +:meth:`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 ` +and +:class:`Laurent series`. + +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 +# +# 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) + + 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((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 series_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.series_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 +