Skip to content
This repository has been archived by the owner on Jan 30, 2023. It is now read-only.

Commit

Permalink
16203: conversion from SR.series to PowerSeries
Browse files Browse the repository at this point in the history
  • Loading branch information
rwst committed Jan 23, 2015
1 parent 9e13c1b commit 3ece975
Show file tree
Hide file tree
Showing 8 changed files with 264 additions and 14 deletions.
1 change: 1 addition & 0 deletions src/doc/en/reference/calculus/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
20 changes: 20 additions & 0 deletions src/sage/rings/power_series_ring.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.<x> = PowerSeriesRing(QQ)
sage: R(s)
1 + x + x^2 + O(x^3)
sage: ex=(gamma(1-y)).series(y,3)
sage: R.<y> = 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`)::
Expand All @@ -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
Expand All @@ -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):
Expand Down
31 changes: 18 additions & 13 deletions src/sage/symbolic/expression.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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 = <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):
"""
Calculate the residue of ``self`` with respect to ``symbol``.
Expand Down Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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::
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -10888,3 +10892,4 @@ cdef operators compatible_relation(operators lop, operators rop) except <operato
return greater
else:
raise TypeError("incompatible relations")

1 change: 1 addition & 0 deletions src/sage/symbolic/ginac.pxd
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,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":
Expand Down
7 changes: 7 additions & 0 deletions src/sage/symbolic/ginac_wrap.h
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down
4 changes: 3 additions & 1 deletion src/sage/symbolic/ring.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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 = {}

Expand Down
4 changes: 4 additions & 0 deletions src/sage/symbolic/series.pxd
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
from sage.symbolic.expression cimport Expression

cdef class SymbolicSeries(Expression):
pass
210 changes: 210 additions & 0 deletions src/sage/symbolic/series.pyx
Original file line number Diff line number Diff line change
@@ -0,0 +1,210 @@
"""
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)

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 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

0 comments on commit 3ece975

Please sign in to comment.