Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
""" Power Series
Sage provides an implementation of dense and sparse power series over any Sage base ring. This is the base class of the implementations of univariate and multivariate power series ring elements in Sage (see also :doc:`power_series_poly`, :doc:`multi_power_series_ring_element`).
AUTHORS:
- William Stein - David Harvey (2006-09-11): added solve_linear_de() method - Robert Bradshaw (2007-04): sqrt, rmul, lmul, shifting - Robert Bradshaw (2007-04): Cython version - Simon King (2012-08): use category and coercion framework, :trac:`13412`
EXAMPLES::
sage: R.<x> = PowerSeriesRing(ZZ) sage: TestSuite(R).run() sage: R([1,2,3]) 1 + 2*x + 3*x^2 sage: R([1,2,3], 10) 1 + 2*x + 3*x^2 + O(x^10) sage: f = 1 + 2*x - 3*x^3 + O(x^4); f 1 + 2*x - 3*x^3 + O(x^4) sage: f^10 1 + 20*x + 180*x^2 + 930*x^3 + O(x^4) sage: g = 1/f; g 1 - 2*x + 4*x^2 - 5*x^3 + O(x^4) sage: g * f 1 + O(x^4)
In Python (as opposed to Sage) create the power series ring and its generator as follows::
sage: R = PowerSeriesRing(ZZ, 'x') sage: x = R.gen() sage: parent(x) Power Series Ring in x over Integer Ring
EXAMPLES:
This example illustrates that coercion for power series rings is consistent with coercion for polynomial rings.
::
sage: poly_ring1.<gen1> = PolynomialRing(QQ) sage: poly_ring2.<gen2> = PolynomialRing(QQ) sage: huge_ring.<x> = PolynomialRing(poly_ring1)
The generator of the first ring gets coerced in as itself, since it is the base ring.
::
sage: huge_ring(gen1) gen1
The generator of the second ring gets mapped via the natural map sending one generator to the other.
::
sage: huge_ring(gen2) x
With power series the behavior is the same.
::
sage: power_ring1.<gen1> = PowerSeriesRing(QQ) sage: power_ring2.<gen2> = PowerSeriesRing(QQ) sage: huge_power_ring.<x> = PowerSeriesRing(power_ring1) sage: huge_power_ring(gen1) gen1 sage: huge_power_ring(gen2) x """
#***************************************************************************** # Copyright (C) 2005 William Stein <wstein@gmail.com> # 2017 Vincent Delecroix <20100.delecroix@gmail.com> # # Distributed under the terms of the GNU General Public License (GPL) # # This code is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # The full text of the GPL is available at: # # http://www.gnu.org/licenses/ #***************************************************************************** from __future__ import absolute_import
import operator
from .infinity import infinity, is_Infinite from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing import sage.rings.polynomial.polynomial_element import sage.misc.misc import sage.arith.all as arith import sage.misc.latex from . import rational_field from . import integer_ring from .integer import Integer from sage.rings.finite_rings.integer_mod_ring import IntegerModRing from sage.misc.superseded import deprecated_function_alias, deprecation from warnings import warn
from sage.categories.fields import Fields _Fields = Fields()
from sage.misc.derivative import multi_derivative
Polynomial = sage.rings.polynomial.polynomial_element.Polynomial_generic_dense
from sage.structure.element cimport AlgebraElement, RingElement, ModuleElement, Element from sage.structure.richcmp cimport richcmp
def is_PowerSeries(x): """ Return True if ``x`` is an instance of a univariate or multivariate power series.
EXAMPLES::
sage: R.<x> = PowerSeriesRing(ZZ) sage: from sage.rings.power_series_ring_element import is_PowerSeries sage: is_PowerSeries(1+x^2) True sage: is_PowerSeries(x-x) True sage: is_PowerSeries(0) False sage: var('x') x sage: is_PowerSeries(1+x^2) False """
cdef class PowerSeries(AlgebraElement): """ A power series. Base class of univariate and multivariate power series. The following methods are available with both types of objects. """ def __init__(self, parent, prec, is_gen=False): """ Initialize a power series. Not for public use. It gets called by the ``PowerSeries_poly`` and ``MPowerSeries`` constructors.
EXAMPLES::
sage: PowerSeriesRing(CC, 'q') Power Series Ring in q over Complex Field with 53 bits of precision sage: T = PowerSeriesRing(GF(3),5,'t'); T Multivariate Power Series Ring in t0, t1, t2, t3, t4 over Finite Field of size 3 """
def __hash__(self): """ Compute a hash of self.
EXAMPLES::
sage: R.<x> = PowerSeriesRing(ZZ) sage: (1+x^2).__hash__() # random 15360174650385709 """
def __reduce__(self): """ EXAMPLES::
sage: K.<t> = PowerSeriesRing(QQ, 5) sage: f = 1 + t - t^3 + O(t^12) sage: loads(dumps(f)) == f True """ return make_element_from_parent_v0, (self._parent, self.polynomial(), self.prec())
def is_sparse(self): """ EXAMPLES::
sage: R.<t> = PowerSeriesRing(ZZ) sage: t.is_sparse() False sage: R.<t> = PowerSeriesRing(ZZ, sparse=True) sage: t.is_sparse() True """
def is_dense(self): """ EXAMPLES::
sage: R.<t> = PowerSeriesRing(ZZ) sage: t.is_dense() True sage: R.<t> = PowerSeriesRing(ZZ, sparse=True) sage: t.is_dense() False """
def is_gen(self): """ Return True if this is the generator (the variable) of the power series ring.
EXAMPLES::
sage: R.<t> = QQ[[]] sage: t.is_gen() True sage: (1 + 2*t).is_gen() False
Note that this only returns True on the actual generator, not on something that happens to be equal to it.
::
sage: (1*t).is_gen() False sage: 1*t == t True """
def _im_gens_(self, codomain, im_gens): """ Return the image of this series under the map that sends the generators to ``im_gens``. This is used internally for computing homomorphisms.
EXAMPLES::
sage: R.<t> = QQ[[]] sage: f = 1 + t + t^2 sage: f._im_gens_(ZZ, [3]) 13 """
cpdef base_extend(self, R): """ Return a copy of this power series but with coefficients in R.
The following coercion uses ``base_extend`` implicitly::
sage: R.<t> = ZZ[['t']] sage: (t - t^2) * Mod(1, 3) t + 2*t^2 """
def change_ring(self, R): """ Change if possible the coefficients of self to lie in R.
EXAMPLES::
sage: R.<T> = QQ[[]]; R Power Series Ring in T over Rational Field sage: f = 1 - 1/2*T + 1/3*T^2 + O(T^3) sage: f.base_extend(GF(5)) Traceback (most recent call last): ... TypeError: no base extension defined sage: f.change_ring(GF(5)) 1 + 2*T + 2*T^2 + O(T^3) sage: f.change_ring(GF(3)) Traceback (most recent call last): ... ZeroDivisionError: Inverse does not exist.
We can only change the ring if there is a ``__call__`` coercion defined. The following succeeds because ``ZZ(K(4))`` is defined.
::
sage: K.<a> = NumberField(cyclotomic_polynomial(3), 'a') sage: R.<t> = K[['t']] sage: (4*t).change_ring(ZZ) 4*t
This does not succeed because ``ZZ(K(a+1))`` is not defined.
::
sage: K.<a> = NumberField(cyclotomic_polynomial(3), 'a') sage: R.<t> = K[['t']] sage: ((a+1)*t).change_ring(ZZ) Traceback (most recent call last): ... TypeError: Unable to coerce a + 1 to an integer """
cpdef _richcmp_(self, right, int op): r""" Comparison of self and ``right``.
We say two approximate power series are equal if they agree for all coefficients up to the *minimum* of the precisions of each. Thus, e.g., `f = 1 + q + O(q^2)` is equal to `g = 1 + O(q)`.
This is how PARI defines equality of power series, but not how Magma defines equality. (Magma would declare `f` and `g` unequal.) The PARI/Sage convention is consistent with the idea that comparison should take place after coercing both elements into a common parent. Hence, in the above example `f` is truncated to `f + O(q)`, which is considered to be equal to `g`, even though the coefficients of `q` are unknown for both series in that comparison.
Comparison is done in dictionary order from lowest degree to highest degree coefficients. This is different than polynomial comparison.
EXAMPLES::
sage: R.<q> = ZZ[[ ]]; R Power Series Ring in q over Integer Ring sage: f=1+q+O(q^2); g = 1+O(q) sage: f == g True sage: 1 - 2*q + q^2 +O(q^3) == 1 - 2*q^2 + q^2 + O(q^4) False
::
sage: R.<t> = ZZ[[]] sage: 1 + t^2 < 2 - t True sage: f = 1 + t + t^7 - 5*t^10 sage: g = 1 + t + t^7 - 5*t^10 + O(t^15) sage: f == f True sage: f < g False sage: f == g True
TESTS:
:trac:`9457` is fixed::
sage: A.<t> = PowerSeriesRing(ZZ) sage: g = t + t^3 + t^5 + O(t^6); g t + t^3 + t^5 + O(t^6) sage: [g == g.add_bigoh(i) for i in range(7)] [True, True, True, True, True, True, True] sage: A(g.polynomial()) == g True
sage: f = t + t^2 + O(t^10) sage: f == f.truncate() True """
def __call__(self, x): """ Implementations *MUST* override this in the derived class.
EXAMPLES::
sage: R.<x> = PowerSeriesRing(ZZ) sage: PowerSeries.__call__(1+x^2,x) Traceback (most recent call last): ... NotImplementedError """
def coefficients(self): """ Return the nonzero coefficients of self.
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ) sage: f = t + t^2 - 10/3*t^3 sage: f.coefficients() [1, 1, -10/3] """
def exponents(self): """ Return the exponents appearing in self with nonzero coefficients.
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ) sage: f = t + t^2 - 10/3*t^3 sage: f.exponents() [1, 2, 3] """
def list(self): """ See this method in derived classes:
- :meth:`sage.rings.power_series_poly.PowerSeries_poly.list`,
- :meth:`sage.rings.multi_power_series_ring_element.MPowerSeries.list`
Implementations *MUST* override this in the derived class.
EXAMPLES::
sage: R.<x> = PowerSeriesRing(ZZ) sage: PowerSeries.list(1+x^2) Traceback (most recent call last): ... NotImplementedError """
def polynomial(self): """ See this method in derived classes:
- :meth:`sage.rings.power_series_poly.PowerSeries_poly.polynomial`,
- :meth:`sage.rings.multi_power_series_ring_element.MPowerSeries.polynomial`
Implementations *MUST* override this in the derived class.
EXAMPLES::
sage: R.<x> = PowerSeriesRing(ZZ) sage: PowerSeries.polynomial(1+x^2) Traceback (most recent call last): ... NotImplementedError """
def __copy__(self): """ Return this power series. Power series are immutable so copy can safely just return the same polynomial.
EXAMPLES::
sage: R.<q> = ZZ[[ ]]; R Power Series Ring in q over Integer Ring sage: f = 1 + 3*q + O(q^10) sage: copy(f) is f # !!! ok since power series are immutable. True """
def base_ring(self): """ Return the base ring that this power series is defined over.
EXAMPLES::
sage: R.<t> = GF(49,'alpha')[[]] sage: (t^2 + O(t^3)).base_ring() Finite Field in alpha of size 7^2 """
def padded_list(self, n=None): """ Return a list of coefficients of self up to (but not including) `q^n`.
Includes 0's in the list on the right so that the list has length `n`.
INPUT:
- ``n`` - (optional) an integer that is at least 0. If ``n`` is not given, it will be taken to be the precision of self, unless this is ``+Infinity``, in which case we just return ``self.list()``.
EXAMPLES::
sage: R.<q> = PowerSeriesRing(QQ) sage: f = 1 - 17*q + 13*q^2 + 10*q^4 + O(q^7) sage: f.list() [1, -17, 13, 0, 10] sage: f.padded_list(7) [1, -17, 13, 0, 10, 0, 0] sage: f.padded_list(10) [1, -17, 13, 0, 10, 0, 0, 0, 0, 0] sage: f.padded_list(3) [1, -17, 13] sage: f.padded_list() [1, -17, 13, 0, 10, 0, 0] sage: g = 1 - 17*q + 13*q^2 + 10*q^4 sage: g.list() [1, -17, 13, 0, 10] sage: g.padded_list() [1, -17, 13, 0, 10] sage: g.padded_list(10) [1, -17, 13, 0, 10, 0, 0, 0, 0, 0] """ else: raise ValueError("n must be at least 0") else:
def prec(self): """ The precision of `...+O(x^r)` is by definition `r`.
EXAMPLES::
sage: R.<t> = ZZ[[]] sage: (t^2 + O(t^3)).prec() 3 sage: (1 - t^2 + O(t^100)).prec() 100 """
def precision_absolute(self): """ Return the absolute precision of this series.
By definition, the absolute precision of `...+O(x^r)` is `r`.
EXAMPLES::
sage: R.<t> = ZZ[[]] sage: (t^2 + O(t^3)).precision_absolute() 3 sage: (1 - t^2 + O(t^100)).precision_absolute() 100 """
def precision_relative(self): """ Return the relative precision of this series, that is the difference between its absolute precision and its valuation.
By convention, the relative precision of `0` (or `O(x^r)` for any `r`) is `0`.
EXAMPLES::
sage: R.<t> = ZZ[[]] sage: (t^2 + O(t^3)).precision_relative() 1 sage: (1 - t^2 + O(t^100)).precision_relative() 100 sage: O(t^4).precision_relative() 0 """ else:
def _repr_(self): """ Return the string representation of this power series.
EXAMPLES::
sage: R.<t> = ZZ[[]] sage: (t^2 + O(t^3))._repr_() 't^2 + O(t^3)'
::
sage: R.<t> = QQ[[]] sage: 1 / (1+2*t +O(t^5)) 1 - 2*t + 4*t^2 - 8*t^3 + 16*t^4 + O(t^5)
::
sage: R.<t> = PowerSeriesRing(QQ, sparse=True) sage: 1 / (1+2*t +O(t^5)) 1 - 2*t + 4*t^2 - 8*t^3 + 16*t^4 + O(t^5) sage: -13/2 * t^3 + 5*t^5 + O(t^10) -13/2*t^3 + 5*t^5 + O(t^10) """ else:
x = "(%s)"%x else: else: else: # end
bigoh = "O(1)" else: return bigoh
def _latex_(self): r""" Return the latex representation of this power series.
EXAMPLES::
sage: R.<t> = QQ[[]] sage: f = -1/2 * t + 2/3*t^2 - 9/7 * t^15 + O(t^20); f -1/2*t + 2/3*t^2 - 9/7*t^15 + O(t^20) sage: latex(f) -\frac{1}{2}t + \frac{2}{3}t^{2} - \frac{9}{7}t^{15} + O(t^{20}) """ if self.prec() is infinity: return "0" else: return "0 + \\cdots" x = "\\left(%s\\right)"%x else: var = "" else: s += repr(x)
bigoh = "O(1)" bigoh = "O(%s)"%(X,) else: return bigoh
def truncate(self, prec=infinity): """ The polynomial obtained from power series by truncation.
EXAMPLES::
sage: R.<I> = GF(2)[[]] sage: f = 1/(1+I+O(I^8)); f 1 + I + I^2 + I^3 + I^4 + I^5 + I^6 + I^7 + O(I^8) sage: f.truncate(5) I^4 + I^3 + I^2 + I + 1 """
cdef _inplace_truncate(self, long prec): return self.truncate(prec)
def add_bigoh(self, prec): r""" Return the power series of precision at most ``prec`` got by adding `O(q^\text{prec})` to `f`, where `q` is the variable.
EXAMPLES::
sage: R.<A> = RDF[[]] sage: f = (1+A+O(A^5))^5; f 1.0 + 5.0*A + 10.0*A^2 + 10.0*A^3 + 5.0*A^4 + O(A^5) sage: f.add_bigoh(3) 1.0 + 5.0*A + 10.0*A^2 + O(A^3) sage: f.add_bigoh(5) 1.0 + 5.0*A + 10.0*A^2 + 10.0*A^3 + 5.0*A^4 + O(A^5) """
def __getitem__(self,n): r""" Return the coefficient of `t^n` in this power series, where `t` is the indeterminate of the power series ring.
If `n` is negative return 0. If `n` is beyond the precision, raise an IndexError.
EXAMPLES::
sage: R.<m> = CDF[[]] sage: f = CDF(pi)^2 + m^3 + CDF(e)*m^4 + O(m^10); f # abs tol 5e-16 9.869604401089358 + 0.0*m + 0.0*m^2 + 1.0*m^3 + 2.718281828459045*m^4 + O(m^10) sage: f[-5] 0.0 sage: f[0] 9.869604401089358 sage: f[4] # abs tol 5e-16 2.718281828459045 sage: f[9] 0.0 sage: f[10] Traceback (most recent call last): ... IndexError: coefficient not known sage: f[1000] Traceback (most recent call last): ... IndexError: coefficient not known """ if n<0: return self.base_ring()(0) c = self.list() if n >= len(c): if self._prec > n: return self.base_ring()(0) else: raise IndexError("coefficient not known") return c[n]
def common_prec(self, f): r""" Return minimum precision of `f` and ``self``.
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ)
::
sage: f = t + t^2 + O(t^3) sage: g = t + t^3 + t^4 + O(t^4) sage: f.common_prec(g) 3 sage: g.common_prec(f) 3
::
sage: f = t + t^2 + O(t^3) sage: g = t^2 sage: f.common_prec(g) 3 sage: g.common_prec(f) 3
::
sage: f = t + t^2 sage: f = t^2 sage: f.common_prec(g) +Infinity """
cdef common_prec_c(self, PowerSeries f): else:
def _mul_prec(self, RingElement right_r): else: else: # sp != infinity else: # endif
def __nonzero__(self): """ Return True if this power series is not equal to 0.
EXAMPLES::
sage: R.<q> = ZZ[[ ]]; R Power Series Ring in q over Integer Ring sage: f = 1 + 3*q + O(q^10) sage: f.is_zero() False sage: (0 + O(q^2)).is_zero() True sage: R(0).is_zero() True sage: (0 + O(q^1000)).is_zero() True """
def is_unit(self): """ Return True if this power series is invertible.
A power series is invertible precisely when the constant term is invertible.
EXAMPLES::
sage: R.<t> = PowerSeriesRing(ZZ) sage: (-1 + t - t^5).is_unit() True sage: (3 + t - t^5).is_unit() False
AUTHORS:
- David Harvey (2006-09-03) """
def inverse(self): """ Return the inverse of self, i.e., self^(-1).
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ, sparse=True) sage: t.inverse() t^-1 sage: type(_) <type 'sage.rings.laurent_series_ring_element.LaurentSeries'> sage: (1-t).inverse() 1 + t + t^2 + t^3 + t^4 + t^5 + t^6 + t^7 + t^8 + ... """
def valuation_zero_part(self): r""" Factor self as `q^n \cdot (a_0 + a_1 q + \cdots)` with `a_0` nonzero. Then this function returns `a_0 + a_1 q + \cdots` .
.. NOTE::
This valuation zero part need not be a unit if, e.g., `a_0` is not invertible in the base ring.
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ) sage: ((1/3)*t^5*(17-2/3*t^3)).valuation_zero_part() 17/3 - 2/9*t^3
In this example the valuation 0 part is not a unit::
sage: R.<t> = PowerSeriesRing(ZZ, sparse=True) sage: u = (-2*t^5*(17-t^3)).valuation_zero_part(); u -34 + 2*t^3 sage: u.is_unit() False sage: u.valuation() 0 """ raise ValueError("power series has no valuation 0 part") else:
cpdef _div_(self, denom_r): """ EXAMPLES::
sage: k.<t> = QQ[[]] sage: t/t 1 sage: (t/(t^3 + 1)) * (t^3 + 1) t + O(t^21) sage: (t^5/(t^2 - 2)) * (t^2 -2 ) t^5 + O(t^25) """ raise ZeroDivisionError("Can't divide by something indistinguishable from 0")
# Algorithm: Cancel common factors of q from top and bottom, # then invert the denominator. We do the cancellation first # because we can only invert a unit (and remain in the ring # of power series).
else:
def __mod__(self, other): """ EXAMPLES::
sage: R.<T> = Qp(7)[[]] sage: f = (48*67 + 46*67^2)*T + (1 + 42*67 + 5*67^3)*T^2 + O(T^3) sage: f % 67 T^2 + O(T^3) """ raise NotImplementedError("Mod on power series ring elements not defined except modulo an integer.")
def shift(self, n): r""" Return this power series multiplied by the power `t^n`. If `n` is negative, terms below `t^n` will be discarded. Does not change this power series.
.. NOTE::
Despite the fact that higher order terms are printed to the right in a power series, right shifting decreases the powers of `t`, while left shifting increases them. This is to be consistent with polynomials, integers, etc.
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ['y'], 't', 5) sage: f = ~(1+t); f 1 - t + t^2 - t^3 + t^4 + O(t^5) sage: f.shift(3) t^3 - t^4 + t^5 - t^6 + t^7 + O(t^8) sage: f >> 2 1 - t + t^2 + O(t^3) sage: f << 10 t^10 - t^11 + t^12 - t^13 + t^14 + O(t^15) sage: t << 29 t^30
AUTHORS:
- Robert Bradshaw (2007-04-18) """
def __lshift__(self, n):
def __rshift__(self, n):
def is_monomial(self): """ Return True if this element is a monomial. That is, if self is `x^n` for some non-negative integer `n`.
EXAMPLES::
sage: k.<z> = PowerSeriesRing(QQ, 'z') sage: z.is_monomial() True sage: k(1).is_monomial() True sage: (z+1).is_monomial() False sage: (z^2909).is_monomial() True sage: (3*z^2909).is_monomial() False """
def is_square(self): """ Return True if this function has a square root in this ring, e.g., there is an element `y` in ``self.parent()`` such that `y^2` equals ``self``.
ALGORITHM: If the base ring is a field, this is true whenever the power series has even valuation and the leading coefficient is a perfect square.
For an integral domain, it attempts the square root in the fraction field and tests whether or not the result lies in the original ring.
EXAMPLES::
sage: K.<t> = PowerSeriesRing(QQ, 't', 5) sage: (1+t).is_square() True sage: (2+t).is_square() False sage: (2+t.change_ring(RR)).is_square() True sage: t.is_square() False sage: K.<t> = PowerSeriesRing(ZZ, 't', 5) sage: (1+t).is_square() False sage: f = (1+t)^100 sage: f.is_square() True """ else:
def sqrt(self, prec=None, extend=False, all=False, name=None): r""" Return a square root of self.
INPUT:
- ``prec`` - integer (default: None): if not None and the series has infinite precision, truncates series at precision prec.
- ``extend`` - bool (default: False); if True, return a square root in an extension ring, if necessary. Otherwise, raise a ValueError if the square root is not in the base power series ring. For example, if ``extend`` is True the square root of a power series with odd degree leading coefficient is defined as an element of a formal extension ring.
- ``name`` - string; if ``extend`` is True, you must also specify the print name of the formal square root.
- ``all`` - bool (default: False); if True, return all square roots of self, instead of just one.
ALGORITHM: Newton's method
.. MATH::
x_{i+1} = \frac{1}{2}( x_i + \mathrm{self}/x_i )
EXAMPLES::
sage: K.<t> = PowerSeriesRing(QQ, 't', 5) sage: sqrt(t^2) t sage: sqrt(1+t) 1 + 1/2*t - 1/8*t^2 + 1/16*t^3 - 5/128*t^4 + O(t^5) sage: sqrt(4+t) 2 + 1/4*t - 1/64*t^2 + 1/512*t^3 - 5/16384*t^4 + O(t^5) sage: u = sqrt(2+t, prec=2, extend=True, name = 'alpha'); u alpha sage: u^2 2 + t sage: u.parent() Univariate Quotient Polynomial Ring in alpha over Power Series Ring in t over Rational Field with modulus x^2 - 2 - t sage: K.<t> = PowerSeriesRing(QQ, 't', 50) sage: sqrt(1+2*t+t^2) 1 + t sage: sqrt(t^2 +2*t^4 + t^6) t + t^3 sage: sqrt(1 + t + t^2 + 7*t^3)^2 1 + t + t^2 + 7*t^3 + O(t^50) sage: sqrt(K(0)) 0 sage: sqrt(t^2) t
::
sage: K.<t> = PowerSeriesRing(CDF, 5) sage: v = sqrt(-1 + t + t^3, all=True); v [1.0*I - 0.5*I*t - 0.125*I*t^2 - 0.5625*I*t^3 - 0.2890625*I*t^4 + O(t^5), -1.0*I + 0.5*I*t + 0.125*I*t^2 + 0.5625*I*t^3 + 0.2890625*I*t^4 + O(t^5)] sage: [a^2 for a in v] [-1.0 + 1.0*t + 0.0*t^2 + 1.0*t^3 + O(t^5), -1.0 + 1.0*t + 0.0*t^2 + 1.0*t^3 + O(t^5)]
A formal square root::
sage: K.<t> = PowerSeriesRing(QQ, 5) sage: f = 2*t + t^3 + O(t^4) sage: s = f.sqrt(extend=True, name='sqrtf'); s sqrtf sage: s^2 2*t + t^3 + O(t^4) sage: parent(s) Univariate Quotient Polynomial Ring in sqrtf over Power Series Ring in t over Rational Field with modulus x^2 - 2*t - t^3 + O(t^4)
TESTS::
sage: R.<x> = QQ[[]] sage: (x^10/2).sqrt() Traceback (most recent call last): ... ValueError: unable to take the square root of 1/2
AUTHORS:
- Robert Bradshaw
- William Stein """ return [ans] else:
raise NotImplementedError('all roots not implemented over a non-integral domain')
# TODO, fix underlying element sqrt() s = u[0].sqrt() if not formal_sqrt: a = self.parent()([s], self.prec()) if all: return [a, -a] else: return a
raise ValueError("the square root generates an extension, so you must specify the name of the square root") if not self.base_ring().is_integral_domain(): raise NotImplementedError('all roots not implemented over a non-integral domain') return [a, -a] else: else: raise ValueError("power series does not have a square root since it has odd valuation.")
else: pr = prec else:
a = a.change_ring(R)
else:
def square_root(self): """ Return the square root of self in this ring. If this cannot be done then an error will be raised.
This function succeeds if and only if ``self``. :meth:`.is_square`
EXAMPLES::
sage: K.<t> = PowerSeriesRing(QQ, 't', 5) sage: (1+t).square_root() 1 + 1/2*t - 1/8*t^2 + 1/16*t^3 - 5/128*t^4 + O(t^5) sage: (2+t).square_root() Traceback (most recent call last): ... ValueError: Square root does not live in this ring. sage: (2+t.change_ring(RR)).square_root() 1.41421356237309 + 0.353553390593274*t - 0.0441941738241592*t^2 + 0.0110485434560398*t^3 - 0.00345266983001244*t^4 + O(t^5) sage: t.square_root() Traceback (most recent call last): ... ValueError: Square root not defined for power series of odd valuation. sage: K.<t> = PowerSeriesRing(ZZ, 't', 5) sage: f = (1+t)^20 sage: f.square_root() 1 + 10*t + 45*t^2 + 120*t^3 + 210*t^4 + O(t^5) sage: f = 1+t sage: f.square_root() Traceback (most recent call last): ... ValueError: Square root does not live in this ring.
AUTHORS:
- Robert Bradshaw """ else:
def nth_root(self, n, prec=None): r""" Return the ``n``-th root of this power series.
INPUT:
- ``n`` -- integer
- ``prec`` -- integer (optional) - precision of the result. Though, if this series has finite precision, then the result can not have larger precision.
EXAMPLES::
sage: R.<x> = QQ[[]] sage: (1+x).nth_root(5) 1 + 1/5*x - 2/25*x^2 + ... + 12039376311816/2384185791015625*x^19 + O(x^20)
sage: (1 + x + O(x^5)).nth_root(5) 1 + 1/5*x - 2/25*x^2 + 6/125*x^3 - 21/625*x^4 + O(x^5)
Check that the results are consistent with taking log and exponential::
sage: R.<x> = PowerSeriesRing(QQ, default_prec=100) sage: p = (1 + 2*x - x^4)**200 sage: p1 = p.nth_root(1000, prec=100) sage: p2 = (p.log()/1000).exp() sage: p1.prec() == p2.prec() == 100 True sage: p1.polynomial() == p2.polynomial() True
Positive characteristic::
sage: R.<u> = GF(3)[[]] sage: p = 1 + 2 * u^2 sage: p.nth_root(4) 1 + 2*u^2 + u^6 + 2*u^8 + u^12 + 2*u^14 + O(u^20) sage: p.nth_root(4)**4 1 + 2*u^2 + O(u^20) """ else:
def cos(self, prec=infinity): r""" Apply cos to the formal power series.
INPUT:
- ``prec`` -- Integer or ``infinity``. The degree to truncate the result to.
OUTPUT:
A new power series.
EXAMPLES:
For one variable::
sage: t = PowerSeriesRing(QQ, 't').gen() sage: f = (t + t**2).O(4) sage: cos(f) 1 - 1/2*t^2 - t^3 + O(t^4)
For several variables::
sage: T.<a,b> = PowerSeriesRing(ZZ,2) sage: f = a + b + a*b + T.O(3) sage: cos(f) 1 - 1/2*a^2 - a*b - 1/2*b^2 + O(a, b)^3 sage: f.cos() 1 - 1/2*a^2 - a*b - 1/2*b^2 + O(a, b)^3 sage: f.cos(prec=2) 1 + O(a, b)^2
If the power series has a non-zero constant coefficient `c`, one raises an error::
sage: g = 2+f sage: cos(g) Traceback (most recent call last): ... ValueError: Can only apply cos to formal power series with zero constant term.
If no precision is specified, the default precision is used::
sage: T.default_prec() 12 sage: cos(a) 1 - 1/2*a^2 + 1/24*a^4 - 1/720*a^6 + 1/40320*a^8 - 1/3628800*a^10 + O(a, b)^12 sage: a.cos(prec=5) 1 - 1/2*a^2 + 1/24*a^4 + O(a, b)^5 sage: cos(a + T.O(5)) 1 - 1/2*a^2 + 1/24*a^4 + O(a, b)^5
TESTS::
sage: cos(a^2 + T.O(5)) 1 - 1/2*a^4 + O(a, b)^5 """
'series with zero constant term.')
def sin(self, prec=infinity): r""" Apply sin to the formal power series.
INPUT:
- ``prec`` -- Integer or ``infinity``. The degree to truncate the result to.
OUTPUT:
A new power series.
EXAMPLES:
For one variable::
sage: t = PowerSeriesRing(QQ, 't').gen() sage: f = (t + t**2).O(4) sage: sin(f) t + t^2 - 1/6*t^3 + O(t^4)
For several variables::
sage: T.<a,b> = PowerSeriesRing(ZZ,2) sage: f = a + b + a*b + T.O(3) sage: sin(f) a + b + a*b + O(a, b)^3 sage: f.sin() a + b + a*b + O(a, b)^3 sage: f.sin(prec=2) a + b + O(a, b)^2
If the power series has a non-zero constant coefficient `c`, one raises an error::
sage: g = 2+f sage: sin(g) Traceback (most recent call last): ... ValueError: Can only apply sin to formal power series with zero constant term.
If no precision is specified, the default precision is used::
sage: T.default_prec() 12 sage: sin(a) a - 1/6*a^3 + 1/120*a^5 - 1/5040*a^7 + 1/362880*a^9 - 1/39916800*a^11 + O(a, b)^12 sage: a.sin(prec=5) a - 1/6*a^3 + O(a, b)^5 sage: sin(a + T.O(5)) a - 1/6*a^3 + O(a, b)^5
TESTS::
sage: sin(a^2 + T.O(5)) a^2 + O(a, b)^5 """
'series with zero constant term.')
def O(self, prec): r""" Return this series plus `O(x^\text{prec})`. Does not change self.
EXAMPLES::
sage: R.<x> = PowerSeriesRing(ZZ) sage: p = 1 + x^2 + x^10; p 1 + x^2 + x^10 sage: p.O(15) 1 + x^2 + x^10 + O(x^15) sage: p.O(5) 1 + x^2 + O(x^5) sage: p.O(-5) Traceback (most recent call last): ... ValueError: prec (= -5) must be non-negative """
def solve_linear_de(self, prec = infinity, b = None, f0 = None): r""" Obtain a power series solution to an inhomogeneous linear differential equation of the form:
.. MATH::
f'(t) = a(t) f(t) + b(t).
INPUT:
- ``self`` - the power series `a(t)`
- ``b`` - the power series `b(t)` (default is zero)
- ``f0`` - the constant term of `f` ("initial condition") (default is 1)
- ``prec`` - desired precision of result (this will be reduced if either a or b have less precision available)
OUTPUT: the power series `f`, to indicated precision
ALGORITHM: A divide-and-conquer strategy; see the source code. Running time is approximately `M(n) \log n`, where `M(n)` is the time required for a polynomial multiplication of length `n` over the coefficient ring. (If you're working over something like `\QQ`, running time analysis can be a little complicated because the coefficients tend to explode.)
.. NOTE::
- If the coefficient ring is a field of characteristic zero, then the solution will exist and is unique.
- For other coefficient rings, things are more complicated. A solution may not exist, and if it does it may not be unique. Generally, by the time the nth term has been computed, the algorithm will have attempted divisions by `n!` in the coefficient ring. So if your coefficient ring has enough 'precision', and if your coefficient ring can perform divisions even when the answer is not unique, and if you know in advance that a solution exists, then this function will find a solution (otherwise it will probably crash).
AUTHORS:
- David Harvey (2006-09-11): factored functionality out from exp() function, cleaned up precision tests a bit
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ, default_prec=10)
::
sage: a = 2 - 3*t + 4*t^2 + O(t^10) sage: b = 3 - 4*t^2 + O(t^7) sage: f = a.solve_linear_de(prec=5, b=b, f0=3/5) sage: f 3/5 + 21/5*t + 33/10*t^2 - 38/15*t^3 + 11/24*t^4 + O(t^5) sage: f.derivative() - a*f - b O(t^4)
::
sage: a = 2 - 3*t + 4*t^2 sage: b = b = 3 - 4*t^2 sage: f = a.solve_linear_de(b=b, f0=3/5) Traceback (most recent call last): ... ValueError: cannot solve differential equation to infinite precision
::
sage: a.solve_linear_de(prec=5, b=b, f0=3/5) 3/5 + 21/5*t + 33/10*t^2 - 38/15*t^3 + 11/24*t^4 + O(t^5) """
# reduce precision to whatever is available from a and b raise ValueError("prec (=%s) must be a non-negative integer" % prec)
def exp(self, prec=None): r""" Return exp of this power series to the indicated precision.
INPUT:
- ``prec`` - integer; default is ``self.parent().default_prec``
ALGORITHM: See :meth:`.solve_linear_de`.
.. NOTE::
- Screwy things can happen if the coefficient ring is not a field of characteristic zero. See :meth:`.solve_linear_de`.
AUTHORS:
- David Harvey (2006-09-08): rewrote to use simplest possible "lazy" algorithm.
- David Harvey (2006-09-10): rewrote to use divide-and-conquer strategy.
- David Harvey (2006-09-11): factored functionality out to solve_linear_de().
- Sourav Sen Gupta, David Harvey (2008-11): handle constant term
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ, default_prec=10)
Check that `\exp(t)` is, well, `\exp(t)`::
sage: (t + O(t^10)).exp() 1 + t + 1/2*t^2 + 1/6*t^3 + 1/24*t^4 + 1/120*t^5 + 1/720*t^6 + 1/5040*t^7 + 1/40320*t^8 + 1/362880*t^9 + O(t^10)
Check that `\exp(\log(1+t))` is `1+t`::
sage: (sum([-(-t)^n/n for n in range(1, 10)]) + O(t^10)).exp() 1 + t + O(t^10)
Check that `\exp(2t + t^2 - t^5)` is whatever it is::
sage: (2*t + t^2 - t^5 + O(t^10)).exp() 1 + 2*t + 3*t^2 + 10/3*t^3 + 19/6*t^4 + 8/5*t^5 - 7/90*t^6 - 538/315*t^7 - 425/168*t^8 - 30629/11340*t^9 + O(t^10)
Check requesting lower precision::
sage: (t + t^2 - t^5 + O(t^10)).exp(5) 1 + t + 3/2*t^2 + 7/6*t^3 + 25/24*t^4 + O(t^5)
Can't get more precision than the input::
sage: (t + t^2 + O(t^3)).exp(10) 1 + t + 3/2*t^2 + O(t^3)
Check some boundary cases::
sage: (t + O(t^2)).exp(1) 1 + O(t) sage: (t + O(t^2)).exp(0) O(t^0)
Handle nonzero constant term (fixes :trac:`4477`)::
sage: R.<x> = PowerSeriesRing(RR) sage: (1 + x + x^2 + O(x^3)).exp() 2.71828182845905 + 2.71828182845905*x + 4.07742274268857*x^2 + O(x^3)
::
sage: R.<x> = PowerSeriesRing(ZZ) sage: (1 + x + O(x^2)).exp() Traceback (most recent call last): ... ArithmeticError: exponential of constant term does not belong to coefficient ring (consider working in a larger ring)
::
sage: R.<x> = PowerSeriesRing(GF(5)) sage: (1 + x + O(x^2)).exp() Traceback (most recent call last): ... ArithmeticError: constant term of power series does not support exponentiation """
def log(self, prec=None): r""" Return log of this power series to the indicated precision.
This works only if the constant term of the power series is 1 or the base ring can take the logarithm of the constant coefficient.
INPUT:
- ``prec`` -- integer; default is ``self.parent().default_prec()``
ALGORITHM: See :meth:`solve_linear_de()`.
.. WARNING::
Screwy things can happen if the coefficient ring is not a field of characteristic zero. See :meth:`solve_linear_de()`.
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ, default_prec=10) sage: (1 + t + O(t^10)).log() t - 1/2*t^2 + 1/3*t^3 - 1/4*t^4 + 1/5*t^5 - 1/6*t^6 + 1/7*t^7 - 1/8*t^8 + 1/9*t^9 + O(t^10)
sage: t.exp().log() t + O(t^10)
sage: (1+t).log().exp() 1 + t + O(t^10)
sage: (-1 + t + O(t^10)).log() Traceback (most recent call last): ... ArithmeticError: constant term of power series is not 1
sage: R.<t> = PowerSeriesRing(RR) sage: (2+t).log().exp() 2.00000000000000 + 1.00000000000000*t + O(t^20) """
else:
def V(self, n): r""" If `f = \sum a_m x^m`, then this function returns `\sum a_m x^{nm}`.
EXAMPLES::
sage: R.<x> = PowerSeriesRing(ZZ) sage: p = 1 + x^2 + x^10; p 1 + x^2 + x^10 sage: p.V(3) 1 + x^6 + x^30 sage: (p+O(x^20)).V(3) 1 + x^6 + x^30 + O(x^60) """ else:
def valuation(self): """ Return the valuation of this power series.
This is equal to the valuation of the underlying polynomial.
EXAMPLES:
Sparse examples::
sage: R.<t> = PowerSeriesRing(QQ, sparse=True) sage: f = t^100000 + O(t^10000000) sage: f.valuation() 100000 sage: R(0).valuation() +Infinity
Dense examples::
sage: R.<t> = PowerSeriesRing(ZZ) sage: f = 17*t^100 +O(t^110) sage: f.valuation() 100 sage: t.valuation() 1 """ return self.polynomial().valuation()
def variable(self): """ Return a string with the name of the variable of this power series.
EXAMPLES::
sage: R.<x> = PowerSeriesRing(Rationals()) sage: f = x^2 + 3*x^4 + O(x^7) sage: f.variable() 'x'
AUTHORS:
- David Harvey (2006-08-08) """
def degree(self): """ Return the degree of this power series, which is by definition the degree of the underlying polynomial.
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ, sparse=True) sage: f = t^100000 + O(t^10000000) sage: f.degree() 100000 """
def derivative(self, *args): r""" The formal derivative of this power series, with respect to variables supplied in args.
Multiple variables and iteration counts may be supplied; see documentation for the global derivative() function for more details.
.. SEEALSO::
:meth:`_derivative`
EXAMPLES::
sage: R.<x> = PowerSeriesRing(QQ) sage: g = -x + x^2/2 - x^4 + O(x^6) sage: g.derivative() -1 + x - 4*x^3 + O(x^5) sage: g.derivative(x) -1 + x - 4*x^3 + O(x^5) sage: g.derivative(x, x) 1 - 12*x^2 + O(x^4) sage: g.derivative(x, 2) 1 - 12*x^2 + O(x^4) """
def __setitem__(self, n, value): """ Called when an attempt is made to change a power series.
EXAMPLES::
sage: R.<t> = ZZ[[]] sage: f = 3 - t^3 + O(t^5) sage: f[1] = 5 Traceback (most recent call last): ... IndexError: power series are immutable """
def laurent_series(self): """ Return the Laurent series associated to this power series, i.e., this series considered as a Laurent series.
EXAMPLES::
sage: k.<w> = QQ[[]] sage: f = 1+17*w+15*w^3+O(w^5) sage: parent(f) Power Series Ring in w over Rational Field sage: g = f.laurent_series(); g 1 + 17*w + 15*w^3 + O(w^5) """
def egf_to_ogf(self): r""" Returns the ordinary generating function power series, assuming self is an exponential generating function power series.
This function is known as ``serlaplace`` in PARI/GP.
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ) sage: f = t + t^2/factorial(2) + 2*t^3/factorial(3) sage: f.egf_to_ogf() t + t^2 + 2*t^3 """
def ogf_to_egf(self): r""" Returns the exponential generating function power series, assuming self is an ordinary generating function power series.
This can also be computed as ``serconvol(f,exp(t))`` in PARI/GP.
EXAMPLES::
sage: R.<t> = PowerSeriesRing(QQ) sage: f = t + t^2 + 2*t^3 sage: f.ogf_to_egf() t + 1/2*t^2 + 1/3*t^3 """
ogf = deprecated_function_alias(15705, egf_to_ogf) egf = deprecated_function_alias(15705, ogf_to_egf)
def __pari__(self): """ Return a PARI representation of this series.
There are currently limits to the possible base rings over which this function works. See the documentation for ``sage.rings.polynomial.polynomial_element.Polynomial.__pari__``
EXAMPLES::
sage: k.<w> = QQ[[]] sage: f = 1+17*w+15*w^3+O(w^5) sage: pari(f) # indirect doctest 1 + 17*w + 15*w^3 + O(w^5) sage: pari(1 - 19*w + w^5) # indirect doctest w^5 - 19*w + 1 sage: R.<x> = Zmod(6)[[]] sage: pari(1 + x + 8*x^3 + O(x^8)) # indirect doctest Mod(1, 6) + Mod(1, 6)*x + Mod(2, 6)*x^3 + O(x^8)
TESTS::
sage: pari(1 + O(x^1)) Mod(1, 6) + O(x) sage: pari(O(x^1)) O(x) sage: pari(O(x^0)) O(x^0) """
def _solve_linear_de(R, N, L, a, b, f0): r""" Internal function used by PowerSeries.solve_linear_de().
INPUT:
- ``R`` - a PolynomialRing
- ``N`` - integer = 0
- ``L`` - integer = 1
- ``a`` - list of coefficients of `a`, any length, all coefficients should belong to base ring of R.
- ``b`` - list of coefficients of `b`, length at least `L` (only first `L` coefficients are used), all coefficients should belong to base ring of R.
- ``f0`` - constant term of `f` (only used if `N == 0`), should belong to base ring of R.
OUTPUT: List of coefficients of `f` (length exactly `L`), where `f` is a solution to the linear inhomogeneous differential equation:
.. MATH::
(t^N f)' = a t^N f + t^{N-1} b + O(t^{N+L-1}).
The constant term of `f` is taken to be f0 if `N == 0`, and otherwise is determined by `N f_0 = b_0`.
ALGORITHM: The algorithm implemented here is inspired by the "fast lazy multiplication algorithm" described in the paper "Lazy multiplication of formal power series" by J van der Hoeven (1997 ISSAC conference).
Our situation is a bit simpler than the one described there, because only one of the series being multiplied needs the lazy treatment; the other one is known fully in advance. I have reformulated the algorithm as an explicit divide-and-conquer recursion. Possibly it is slower than the one described by van der Hoeven, by a constant factor, but it seems easier to implement. Also, because we repeatedly split in half starting at the top level, instead of repeatedly doubling in size from the bottom level, it's easier to avoid problems with "overshooting" in the last iteration.
The idea is to split the problem into two instances with `L` about half the size. Take `L'` to be the ceiling of `(L/2)`. First recursively find `g` modulo `t^{L'}` such that
.. MATH::
(t^N g)' = a t^N g + t^{N-1} b + O(t^{N+L'-1}).
Next we want to find `h` modulo `t^{L-L'}` such that `f = g + t^{L'} h` is a solution of the original equation. We can find a suitable `h` by recursively solving another differential equation of the same form, namely
.. MATH::
(t^{N+L'} h)' = a t^{N+L'} h + c t^{N+L'-1} + O(t^{N+L'-1}),
where `c` is determined by
.. MATH::
(t^N g)' - a t^N g - t^{N-1} b = -c t^{N+L'-1} + O(t^{N+L-1}).
Once `g` is known, `c` can be recovered from this relation by computing the coefficients of `t^j` of the product `a g` for `j` in the range `L'-1 \leq j < L-1`.
At the bottom of the recursion we have `L = 1`, in which case the solution is simply given by `f_0 = b_0/N` (or by the supplied initial condition `f_0` when `N == 0`).
The algorithm has to do one multiplication of length `N`, two of length `N/2`, four of length `N/4`, etc, hence giving runtime `O(M(N) \log N)`.
AUTHORS:
- David Harvey (2006-09-11): factored out of PowerSeries.exp(). """ # base case else:
# todo: perhaps next loop could be made more efficient
def make_powerseries_poly_v0(parent, f, prec, is_gen): # This is only used to unpickle old pickles. The new pickling # works differently! from . import power_series_poly return power_series_poly.PowerSeries_poly(parent, f, prec, 0, is_gen)
def make_element_from_parent_v0(parent, *args): # This is only used to unpickle old pickles. The new pickling # works differently! return parent(*args) |