# source:sage/rings/power_series_ring_element.pyx@5428:e4e8a90e85d2

Revision 5428:e4e8a90e85d2, 63.3 KB checked in by Robert Bradshaw <robertwb@…>, 6 years ago (diff)

optimize sqrt some more

Line
1"""
2Power Series
3
4SAGE provides an implementation of dense and sparse power series over
5any SAGE base ring.
6
7AUTHOR:
8   -- William Stein
9   -- David Harvey (2006-09-11): added solve_linear_de() method
10   -- Robert Bradshaw (2007-04): sqrt, rmul, lmul, shifting
11   -- Robert Bradshaw (2007-04): SageX version
12
13EXAMPLE:
14    sage: R.<x> = PowerSeriesRing(ZZ)
15    sage: R([1,2,3])
16    1 + 2*x + 3*x^2
17    sage: R([1,2,3], 10)
18    1 + 2*x + 3*x^2 + O(x^10)
19    sage: f = 1 + 2*x - 3*x^3 + O(x^4); f
20    1 + 2*x - 3*x^3 + O(x^4)
21    sage: f^10
22    1 + 20*x + 180*x^2 + 930*x^3 + O(x^4)
23    sage: g = 1/f; g
24    1 - 2*x + 4*x^2 - 5*x^3 + O(x^4)
25    sage: g * f
26    1 + O(x^4)
27
28In Python (as opposted to SAGE) create the power series ring and its
29generator as follows:
30
31    sage: R, x = objgen(PowerSeriesRing(ZZ, 'x'))
32    sage: parent(x)
33    Power Series Ring in x over Integer Ring
34
35EXAMPLE: COERCION
36This example illustrates that coercion for power series rings
37is consistent with coercion for polynomial rings.
38
39    sage: poly_ring1.<gen1> = PolynomialRing(QQ)
40    sage: poly_ring2.<gen2> = PolynomialRing(QQ)
41    sage: huge_ring.<x> = PolynomialRing(poly_ring1)
42
43The generator of the first ring gets coerced in as itself,
44since it is the base ring.
45    sage: huge_ring(gen1)
46    gen1
47
48The generator of the second ring gets mapped via the
49natural map sending one generator to the other.
50    sage: huge_ring(gen2)
51    x
52
53With power series the behavior is the same.
54    sage: power_ring1.<gen1> = PowerSeriesRing(QQ)
55    sage: power_ring2.<gen2> = PowerSeriesRing(QQ)
56    sage: huge_power_ring.<x> = PowerSeriesRing(power_ring1)
57    sage: huge_power_ring(gen1)
58    gen1
59    sage: huge_power_ring(gen2)
60    x
61"""
62
63#*****************************************************************************
64#       Copyright (C) 2005 William Stein <wstein@gmail.com>
65#
67#
68#    This code is distributed in the hope that it will be useful,
69#    but WITHOUT ANY WARRANTY; without even the implied warranty of
70#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
71#    General Public License for more details.
72#
73#  The full text of the GPL is available at:
74#
76#*****************************************************************************
77
78include "../ext/stdsage.pxi"
79
80import operator
81
82from infinity import infinity, is_Infinite
83from sage.rings.polynomial.polynomial_ring import PolynomialRing
84import sage.rings.polynomial.polynomial_element
85import power_series_ring
86import sage.misc.misc
87import ring_element
88import arith
89import sage.misc.latex
90import sage.structure.coerce
91import rational_field, integer_ring
92from integer import Integer
93from integer_mod_ring import IntegerModRing
94import sage.libs.pari.all
95from sage.libs.all import PariError
96from sage.misc.functional import sqrt, log
97from sage.rings.arith import integer_ceil as ceil
98
99from sage.rings.ring import is_Field
100
101Polynomial = sage.rings.polynomial.polynomial_element.Polynomial_generic_dense
102
103from sage.structure.element cimport AlgebraElement, RingElement, ModuleElement, Element
104
105def is_PowerSeries(x):
106    return isinstance(x, PowerSeries)
107
108cdef class PowerSeries(AlgebraElement):
109    """
110    A power series.
111    """
112
113    def __init__(self, parent, prec, is_gen=False):
114        """
115        Initialize a power series.
116        """
117        AlgebraElement.__init__(self, parent)
118        self.__is_gen = is_gen
119        if not (prec is infinity):
120            prec = int(prec)
121        self._prec = prec
122
123    def __hash__(self):
124        return hash(self.polynomial())
125
126    def __reduce__(self):
127        """
128        EXAMPLES:
129            sage: K.<t> = PowerSeriesRing(QQ, 5)
130            sage: f = 1 + t - t^3 + O(t^12)
132            True
133        """
134        return make_element_from_parent_v0, (self._parent, self.polynomial(), self.prec())
135
136    def is_sparse(self):
137        """
138        EXAMPLES:
139            sage: R.<t> = PowerSeriesRing(ZZ)
140            sage: t.is_sparse()
141            False
142            sage: R.<t> = PowerSeriesRing(ZZ, sparse=True)
143            sage: t.is_sparse()
144            True
145        """
146        return self._parent.is_sparse()
147
148    def is_dense(self):
149        """
150        EXAMPLES:
151            sage: R.<t> = PowerSeriesRing(ZZ)
152            sage: t.is_dense()
153            True
154            sage: R.<t> = PowerSeriesRing(ZZ, sparse=True)
155            sage: t.is_dense()
156            False
157        """
158        return self._parent.is_dense()
159
160    def is_gen(self):
161        """
162        Returns True if this the generator (the variable) of the power series ring.
163
164        EXAMPLES:
165            sage: R.<t> = QQ[[]]
166            sage: t.is_gen()
167            True
168            sage: (1 + 2*t).is_gen()
169            False
170
171        Note that this only returns true on the actual generator, not on
172        something that happens to be equal to it.
173            sage: (1*t).is_gen()
174            False
175            sage: 1*t == t
176            True
177        """
178        return bool(self.__is_gen)
179
180    def _im_gens_(self, codomain, im_gens):
181        """
182        Returns the image of this series under the map that sends the generators
183        to im_gens.   This is used internally for computing homomorphisms.
184
185        EXAMPLES:
186            sage: R.<t> = QQ[[]]
187            sage: f = 1 + t + t^2
188            sage: f._im_gens_(ZZ, [3])
189            13
190        """
191        return codomain(self(im_gens[0]))
192
193    def base_extend(self, R):
194        """
195        Return a copy of this power series but with coefficients in R.
196
197        The following coercion uses base_extend implicitly:
198            sage: R.<t> = ZZ[['t']]
199            sage: (t - t^2) * Mod(1, 3)
200            t + 2*t^2
201        """
202        S = self._parent.base_extend(R)
203        return S(self)
204
205    def change_ring(self, R):
206        """
207        Change if possible the coefficients of self to lie in R.
208
209        EXAMPLES:
210            sage: R.<T> = QQ[[]]; R
211            Power Series Ring in T over Rational Field
212            sage: f = 1 - 1/2*T + 1/3*T^2 + O(T^3)
213            sage: f.base_extend(GF(5))
214            Traceback (most recent call last):
215            ...
216            TypeError: no base extension defined
217            sage: f.change_ring(GF(5))
218            1 + 2*T + 2*T^2 + O(T^3)
219            sage: f.change_ring(GF(3))
220            Traceback (most recent call last):
221            ...
222            ZeroDivisionError: Inverse does not exist.
223
224        We can only change irng if there is a __call__ coercion defined.
225        The following succeeds because ZZ(K(4)) is defined.
226
227            sage: K.<a> = NumberField(cyclotomic_polynomial(3), 'a')
228            sage: R.<t> = K[['t']]
229            sage: (4*t).change_ring(ZZ)
230            4*t
231
232        This does not succeed because ZZ(K(a+1)) is not defined.
233            sage: K.<a> = NumberField(cyclotomic_polynomial(3), 'a')
234            sage: R.<t> = K[['t']]
235            sage: ((a+1)*t).change_ring(ZZ)
236            Traceback (most recent call last):
237            ...
238            TypeError: Unable to coerce a + 1 to an integer
239
240        """
241        S = self._parent.change_ring(R)
242        return S(self)
243
244    def __cmp__(left, right):
245        return (<Element>left)._cmp(right)
246
247    cdef int _cmp_c_impl(self, Element right) except -2:
248        r"""
249        Comparison of self and right.
250
251        We say two approximate power series are equal, if they agree
252        for all coefficients up to the *minimum* of the precisions of
253        each.  Thus, e.g., $f=1+q+O(q^2)$ is equal to $g=1+O(q)$.
254        This is how PARI defines equality of power series, but not how
255        MAGMA defines equality.  (MAGMA would declare f and g
256        unequal.)  I side with PARI, because even if $g=1+q+O(q^2)$,
257        we don't really know whether f equals g, since we don't know
258        the coefficients of $q^2$.
259
260        Comparison is done in dictionary order from lowest degree to
261        highest degree coefficients (this is different than
262        polynomials).
263
264        EXAMPLES:
265            sage: R.<q> = ZZ[[ ]]; R
266            Power Series Ring in q over Integer Ring
267            sage: f=1+q+O(q^2); g = 1+O(q)
268            sage: f == g
269            True
270            sage: 1 - 2*q + q^2 +O(q^3) == 1 - 2*q^2 + q^2 + O(q^4)
271            False
272        """
273        # A very common case throughout code
274        if PY_TYPE_CHECK(right, int):
275            return self.is_zero()
276
277        prec = self.common_prec(right)
278        x = self.list()
279        y = right.list()
280        if not (prec is infinity):
281            x = x[:prec]
282            y = y[:prec]
283        return cmp(x,y)
284
285    def __call__(self, x):   # you *MUST* overrride this in the derived class
286        raise NotImplementedError
287
288    def list(self):          # you *MUST* overrride this in the derived class
289        raise NotImplementedError
290
291    def polynomial(self):          # you *MUST* overrride this in the derived class
292        raise NotImplementedError
293
294    def __setitem__(self, n, value):   # you *MUST* overrride this in the derived class
295        raise NotImplementedError
296
297    def __copy__(self):
298        """
299        Return this power series.  Power series are immutable so copy
300        can safely just return the same polynomial.
301
302        EXAMPLES:
303            sage: R.<q> = ZZ[[ ]]; R
304            Power Series Ring in q over Integer Ring
305            sage: f = 1 + 3*q + O(q^10)
306            sage: copy(f) is f       # !!! ok since power series are immutable.
307            True
308        """
309        return self
310
311    def base_ring(self):
312        """
313        Return the base ring that this power series is defined over.
314
315        EXAMPLES:
316            sage: R.<t> = GF(49,'alpha')[[]]
317            sage: (t^2 + O(t^3)).base_ring()
318            Finite Field in alpha of size 7^2
319        """
320        return self._parent.base_ring()
321
323        """
324        Return list of coefficients of self up to (but not include $q^n$).
325
326        Includes 0's in the list on the right so that the list has
327        length $n$.
328
329        EXAMPLES:
330            sage: R.<q> = PowerSeriesRing(QQ)
331            sage: f = 1 - 17*q + 13*q^2 + 10*q^4 + O(q^7)
332            sage: f.list()
333            [1, -17, 13, 0, 10]
335            [1, -17, 13, 0, 10, 0, 0]
337            [1, -17, 13, 0, 10, 0, 0, 0, 0, 0]
339            [1, -17, 13]
340        """
341        v = self.list()
342        if len(v) < n:
343            z = self._parent.base_ring()(0)
344            return v + [z]*(n - len(v))
345        else:
346            return v[:int(n)]
347
348    def prec(self):
349        """
350        The precision of $...+O(x^r)$ is by definition $r$.
351
352        EXAMPLES:
353            sage: R.<t> = ZZ[[]]
354            sage: (t^2 + O(t^3)).prec()
355            3
356            sage: (1 - t^2 + O(t^100)).prec()
357            100
358        """
359        return self._prec
360
361    def _repr_(self):
362        """
363        Return string represenation of this power series.
364
365        EXAMPLES:
366            sage: R.<t> = ZZ[[]]
367            sage: (t^2 + O(t^3))._repr_()
368            't^2 + O(t^3)'
369
370            sage: R.<t> = QQ[[]]
371            sage: 1 / (1+2*t +O(t^5))
372            1 - 2*t + 4*t^2 - 8*t^3 + 16*t^4 + O(t^5)
373
374            sage: R.<t> = PowerSeriesRing(QQ, sparse=True)
375            sage: 1 / (1+2*t +O(t^5))
376            1 - 2*t + 4*t^2 - 8*t^3 + 16*t^4 + O(t^5)
377            sage: -13/2 * t^3  + 5*t^5 + O(t^10)
378            -13/2*t^3 + 5*t^5 + O(t^10)
379
380        """
381        if self.is_zero():
382            if self.prec() is infinity:
383                return "0"
384            else:
385                return "O(%s^%s)"%(self._parent.variable_name(),self.prec())
386
387        atomic_repr = self._parent.base_ring().is_atomic_repr()
388        X = self._parent.variable_name()
389
390        s = " "
391        if self.is_sparse():
392            f = self.polynomial()
393            m = f.degree() + 1
394            d = f._dict_unsafe()
395            coeffs = list(d.iteritems())
396            coeffs.sort()
397            for (n, x) in coeffs:
398                x = repr(x)
399                if x != '0':
400                    if s != ' ':
401                        s += " + "
402                    if not atomic_repr and n > 0 and (x.find("+") != -1 or x.find("-") != -1):
403                        x = "(%s)"%x
404                    if n > 1:
405                        var = "*%s^%s"%(X,n)
406                    elif n==1:
407                        var = "*%s"%X
408                    else:
409                        var = ""
410                    s += "%s%s"%(x,var)
411        else:
412            v = self.list()
413            m = len(v)
414            first = True
415            for n in xrange(m):
416                x = v[n]
417                x = repr(x)
418                if x != '0':
419                    if not first:
420                        s += " + "
421                    if not atomic_repr and n > 0 and (x[1:].find("+") != -1 or x[1:].find("-") != -1):
422                        x = "(%s)"%x
423                    if n > 1:
424                        var = "*%s^%s"%(X,n)
425                    elif n==1:
426                        var = "*%s"%X
427                    else:
428                        var = ""
429                    s += "%s%s"%(x,var)
430                    first = False
431        # end
432
433        if atomic_repr:
434            s = s.replace(" + -", " - ")
435        s = s.replace(" 1*"," ")
436        s = s.replace(" -1*", " -")
437        if not (self._prec is infinity):
438            if self._prec == 0:
439                bigoh = "O(1)"
440            elif self._prec == 1:
441                bigoh = "O(%s)"%self._parent.variable_name()
442            else:
443                bigoh = "O(%s^%s)"%(self._parent.variable_name(),self._prec)
444            if s==" ":
445                return bigoh
446            s += " + %s"%bigoh
447        return s[1:]
448
449    def _latex_(self):
450        r"""
451        Return latex representation of this power series.
452
453        EXAMPLES:
454            sage: R.<t> = QQ[[]]
455            sage: f = -1/2 * t + 2/3*t^2 + -9/7 * t^15 + O(t^20); f
456            -1/2*t + 2/3*t^2 - 9/7*t^15 + O(t^20)
457            sage: latex(f)
458            -\frac{1}{2}t + \frac{2}{3}t^{2} - \frac{9}{7}t^{15} + O(\text{t}^{20})
459        """
460        if self.is_zero():
461            if self.prec() is infinity:
462                return "0"
463            else:
464                return "0 + \\cdots"
465        s = " "
466        v = self.list()
467        m = len(v)
468        X = self._parent.variable_name()
469        atomic_repr = self._parent.base_ring().is_atomic_repr()
470        first = True
471        for n in xrange(m):
472            x = v[n]
473            x = sage.misc.latex.latex(x)
474            if x != '0':
475                if not first:
476                    s += " + "
477                if not atomic_repr and n > 0 and (x[1:].find("+") != -1 or x[1:].find("-") != -1):
478                    x = "\\left(%s\\right)"%x
479                if n > 1:
480                    var = "%s^{%s}"%(X,n)
481                elif n==1:
482                    var = "%s"%X
483                else:
484                    var = ""
485                if n > 0:
486                    s += "%s|%s"%(x,var)
487                else:
488                    s += repr(x)
489                first = False
490
491        if atomic_repr:
492            s = s.replace(" + -", " - ")
493        s = s.replace(" -1|", " -")
494        s = s.replace(" 1|"," ")
495        s = s.replace("|","")
496        if not (self._prec is infinity):
497            if self._prec == 0:
498                bigoh = "O(1)"
499            elif self._prec == 1:
500                bigoh = "O(%s)"%sage.misc.latex.latex(self._parent.variable_name())
501            else:
502                bigoh = "O(%s^{%s})"%(sage.misc.latex.latex(self._parent.variable_name()),self._prec)
503            if s == " ":
504                return bigoh
505            s += " + %s"%bigoh
506        return s[1:]
507
508
509    def truncate(self, prec=infinity):
510        """
511        The polynomial obtained from power series by truncation.
512
513        EXAMPLES:
514            sage: R.<I> = GF(2)[[]]
515            sage: f = 1/(1+I+O(I^8)); f
516            1 + I + I^2 + I^3 + I^4 + I^5 + I^6 + I^7 + O(I^8)
517            sage: f.truncate(5)
518            I^4 + I^3 + I^2 + I + 1
519        """
520        if prec is infinity:
521            prec = self._prec
522        a = self.list()
523        v = [a[i] for i in range(min(prec, len(a)))]
524        return self._parent._poly_ring()(v)
525
527        r"""
528        Returns the power series of precision at most prec got by
529        adding $O(q^\text{prec})$ to f, where q is the variable.
530
531        EXAMPLES:
532            sage: R.<A> = RDF[[]]
533            sage: f = (1+A+O(A^5))^5; f
534            1.0 + 5.0*A + 10.0*A^2 + 10.0*A^3 + 5.0*A^4 + O(A^5)
536            1.0 + 5.0*A + 10.0*A^2 + O(A^3)
537        """
538        if prec is infinity or prec >= self.prec():
539            return self
540        a = self.list()
541        v = [a[i] for i in range(min(prec, len(a)))]
542        return self._parent(v, prec)
543
544    def __getitem__(self,n):
545        r"""
546        Return the coefficient of $t^n$ in this power series, where
547        $t$ is the indeterminate of the power series ring.
548
549        If n is negative returns 0.  If n is beyond the precision,
550        raises an IndexError.
551
552        EXAMPLES:
553            sage: R.<m> = CDF[[]]
554            sage: f = CDF(pi)^2 + m^3 + CDF(e)*m^4 + O(m^10); f
555            9.86960440109 + 1.0*m^3 + 2.71828182846*m^4 + O(m^10)
556            sage: f[-5]
557            0
558            sage: f[0]
559            9.86960440109
560            sage: f[4]
561            2.71828182846
562            sage: f[9]
563            0
564            sage: f[10]
565            Traceback (most recent call last):
566            ...
567            IndexError: coefficient not known
568            sage: f[1000]
569            Traceback (most recent call last):
570            ...
571            IndexError: coefficient not known
572        """
573        if n<0:
574            return self.base_ring()(0)
575        c = self.list()
576        if n >= len(c):
577            if self._prec > n:
578                return self.base_ring()(0)
579            else:
580                raise IndexError, "coefficient not known"
581        return c[n]
582
583    def common_prec(self, f):
584        r"""
585        Returns minimum precision of $f$ and self.
586
587        EXAMPLES:
588            sage: R.<t> = PowerSeriesRing(QQ)
589
590            sage: f = t + t^2 + O(t^3)
591            sage: g = t + t^3 + t^4 + O(t^4)
592            sage: f.common_prec(g)
593            3
594            sage: g.common_prec(f)
595            3
596
597            sage: f = t + t^2 + O(t^3)
598            sage: g = t^2
599            sage: f.common_prec(g)
600            3
601            sage: g.common_prec(f)
602            3
603
604            sage: f = t + t^2
605            sage: f = t^2
606            sage: f.common_prec(g)
607            +Infinity
608        """
609        if self.prec() is infinity:
610            return f.prec()
611        elif f.prec() is infinity:
612            return self.prec()
613        return min(self.prec(), f.prec())
614
615    cdef common_prec_c(self, PowerSeries f):
616        if self._prec is infinity:
617            return f._prec
618        elif f._prec is infinity:
619            return self._prec
620        elif self._prec < f._prec:
621            return self._prec
622        else:
623            return f._prec
624
625    cdef RingElement _mul_c_impl(self, RingElement right_r):
626        # TODO: doctest
627        cdef PowerSeries right = <PowerSeries>right_r
628        if self.is_zero():
629            return self
630        if right.is_zero():
631            return right
632        sp = self._prec
633        rp = right._prec
634        if sp is infinity:
635            if rp is infinity:
636                prec = infinity
637            else:
638                prec = rp + self.valuation()
639        else:  # sp != infinity
640            if rp is infinity:
641                prec = sp + right.valuation()
642            else:
643                prec = min(rp + self.valuation(), sp + right.valuation())
644        # endif
645        return self._mul_(right, prec) # ???
646
647    def is_zero(self):
648        """
649        Return True if this power series equals 0.
650
651        EXAMPLES:
652            sage: R.<q> = ZZ[[ ]]; R
653            Power Series Ring in q over Integer Ring
654            sage: f = 1 + 3*q + O(q^10)
655            sage: f.is_zero()
656            False
657            sage: (0 + O(q^2)).is_zero()
658            True
659            sage: R(0).is_zero()
660            True
661            sage: (0 + O(q^1000)).is_zero()
662            True
663        """
664        return self.polynomial().is_zero()
665
666    def is_unit(self):
667        """
668        Returns whether this power series is invertible, which is the
669        case precisely when the constant term is invertible.
670
671        EXAMPLES:
672            sage: R.<t> = PowerSeriesRing(ZZ)
673            sage: (-1 + t - t^5).is_unit()
674            True
675            sage: (3 + t - t^5).is_unit()
676            False
677
678        AUTHOR: David Harvey (2006-09-03)
679        """
680        return self[0].is_unit()
681
682    def __invert__(self):
683        """
684        Inverse of the power series (i.e. a series Y such that XY = 1).
685        The first nonzero coefficient must be a unit in the coefficient ring.
686        If the valuation of the series is positive, this function will return
687        a Laurent series.
688
689        ALGORITHM:
690            Uses Newton's method. Complexity is around $O(M(n) \log n)$,
691            where $n$ is the precision and $M(n)$ is the time required to
692            multiply polynomials of length $n$.
693
694        EXAMPLES:
695            sage: R.<q> = QQ[[]]
696            sage: 1/(1+q + O(q**2))
697            1 - q + O(q^2)
698            sage: 1/(1+q)
699            1 - q + q^2 - q^3 + q^4 - q^5 + q^6 - q^7 + q^8 - q^9 + q^10 - q^11 + q^12 - q^13 + q^14 - q^15 + q^16 - q^17 + q^18 - q^19 + O(q^20)
700            sage: prec = R.default_prec(); prec
701            20
702            sage: R.set_default_prec(5)
703            sage: 1/(1+q)
704            1 - q + q^2 - q^3 + q^4 + O(q^5)
705
706            sage: 1/(q + q^2)
707             q^-1 - 1 + q - q^2 + q^3 + O(q^4)
708            sage: g = 1/(q + q^2 + O(q^5))
709            sage: g; g.parent()
710            q^-1 - 1 + q - q^2 + O(q^3)
711            Laurent Series Ring in q over Rational Field
712
713            sage: 1/g
714            q + q^2 + O(q^5)
715            sage: (1/g).parent()
716            Laurent Series Ring in q over Rational Field
717
718            sage: 1/(2 + q)
719             1/2 - 1/4*q + 1/8*q^2 - 1/16*q^3 + 1/32*q^4 + O(q^5)
720
721            sage: R.<q> = QQ[['q']]
722            sage: R.set_default_prec(5)
723            sage: f = 1 + q + q^2 + O(q^50)
724            sage: f/10
725            1/10 + 1/10*q + 1/10*q^2 + O(q^50)
726            sage: f/(10+q)
727            1/10 + 9/100*q + 91/1000*q^2 - 91/10000*q^3 + 91/100000*q^4 + O(q^5)
728
729            sage: R.<t> = PowerSeriesRing(QQ, sparse=True)
730            sage: u = 17 + 3*t^2 + 19*t^10 + O(t^12)
731            sage: v = ~u; v
732            1/17 - 3/289*t^2 + 9/4913*t^4 - 27/83521*t^6 + 81/1419857*t^8 - 1587142/24137569*t^10 + O(t^12)
733            sage: u*v
734            1 + O(t^12)
735
736        AUTHORS:
737            -- David Harvey (2006-09-09): changed to use Newton's method
738
739        """
740        if self == 1:
741            return self
742        prec = self.prec()
743        if prec is infinity and self.degree() > 0:
744            prec = self._parent.default_prec()
745        if self.valuation() > 0:
746            u = ~self.valuation_zero_part()    # inverse of unit part
747            R = self._parent.laurent_series_ring()
748            return R(u, -self.valuation())
749
751        # and then iteratively compute $x' = 2x - Ax^2$, where $A$ is the
752        # series we're trying to invert.
753
754        try:
755            first_coeff = ~self[0]
756        except ValueError, ZeroDivisionError:
757            raise ZeroDivisionError, "leading coefficient must be a unit"
758
759        if prec is infinity:
760            return self._parent(first_coeff, prec=prec)
761
762        A = self.truncate()
763        R = A.parent()     # R is the corresponding polynomial ring
764        current = R(first_coeff)
765
766        # todo: in the case that the underlying polynomial ring is
767        # implemented via NTL, the truncate() method should use NTL's
768        # methods. Currently it is very slow because it uses generic code
769        # that has to pull all the data in and out of the polynomials.
770
771        # todo: also, NTL has built-in series inversion. We should use
772        # that when available.
773
774        for next_prec in sage.misc.misc.newton_method_sizes(prec)[1:]:
775            z = current.square() * A.truncate(next_prec)
776            current = 2*current - z.truncate(next_prec)
777
778        return self._parent(current, prec=prec)
779
780        # Here is the old code, which uses a simple recursion, and is
781        # asymptotically inferior:
782        #
783        #a0 = self[0]
784        #try:
785        #    b = [~a0]
786        #except ValueError, ZeroDivisionError:
787        #    raise ZeroDivisionError, "leading coefficient must be a unit"
788
789        ##  By multiplying through we may assume that the leading coefficient
790        ##  of f=self is 1.   If f = 1 + a_1*q + a_2*q^2 + ..., then we have
791        ##  the following recursive formula for the coefficients b_n of the
792        ##  expansion of f^(-1):
793        ##        b_n = -b_0*(b_{n-1}*a_1 + b_{n-2}*a_2 + ... + b_0 a_n).
794        #if self.degree() > 0:
795        #    a = self.list()
796        #    for n in range(1,prec):
797        #        b.append(-b[0]*sum([b[n-i]*a[i] for i in range(1,n+1) if i < len(a)]))
798        #return self.parent()(b, prec=prec)
799
800    def valuation_zero_part(self):
801        r"""
802        Factor self as as $q^n\cdot (a_0 + a_1 q + \cdots)$ with $a_0$
803        nonzero.  Then this function returns $a_0 + a_1 q + \cdots$.
804
805        NOTE: this valuation zero part need not be a unit if, e.g.,
806        $a_0$ is not invertible in the base ring.
807
808        EXAMPLES:
809            sage: R.<t> = PowerSeriesRing(QQ)
810            sage: ((1/3)*t^5*(17-2/3*t^3)).valuation_zero_part()
811            17/3 - 2/9*t^3
812
813        In this example the valuation 0 part is not a unit:
814            sage: R.<t> = PowerSeriesRing(ZZ, sparse=True)
815            sage: u = (-2*t^5*(17-t^3)).valuation_zero_part(); u
816            -34 + 2*t^3
817            sage: u.is_unit()
818            False
819            sage: u.valuation()
820            0
821        """
822        if self.is_zero():
823            raise ValueError, "power series has no valuation 0 part"
824        n = self.valuation()
825        if n == 0:
826            return self
827        elif self.is_dense():
828            v = self.list()[int(n):]
829        else:
830            n = int(n)
831            v = {}
832            for k, x in self.dict().iteritems():
833                if k >= n:
834                    v[k-n] = x
835        return self._parent(v, self.prec()-n)
836
837    cdef RingElement _div_c_impl(self, RingElement denom_r):
838        """
839        EXAMPLES:
840            sage: k.<t> = QQ[[]]
841            sage: t/t
842            1
843            sage: (t/(t^3 + 1)) * (t^3 + 1)
844            t + O(t^21)
845            sage: (t^5/(t^2 - 2)) * (t^2 -2 )
846            t^5 + O(t^25)
847        """
848        denom = <PowerSeries>denom_r
849        if denom.is_zero():
850            raise ZeroDivisionError, "Can't divide by something indistinguishable from 0"
851        u = denom.valuation_zero_part()
852        inv = ~u  # inverse
853
854        v = denom.valuation()
855        if v > self.valuation():
856            R = self._parent.laurent_series_ring()
857            return R(self)/R(denom)
858
859        # Algorithm: Cancel common factors of q from top and bottom,
860        # then invert the denominator.  We do the cancellation first
861        # because we can only invert a unit (and remain in the ring
862        # of power series).
863
864        if v > 0:
865            num = self >> v
866        else:
867            num = self
868        return num*inv
869
870    def __mod__(self, other):
871        """
872        EXAMPLES:
873            sage: R.<T> = Qp(7)[[]]
874            sage: f = (48*67 + 46*67^2)*T + (1 + 42*67 + 5*67^3)*T^2 + O(T^3)
875            sage: f % 67
876            T^2 + O(T^3)
877        """
878        if isinstance(other,(int,Integer,long)):
879            return power_series_ring.PowerSeriesRing(IntegerModRing(other), self.variable())(self)
880        raise NotImplementedError, "Mod on power series ring elements not defined except modulo an integer."
881
882    def shift(self, n):
883        r"""
884        Returns this power series multiplied by the power $t^n$. If $n$
885        is negative, terms below $t^n$ will be discarded. Does not
886        change this power series.
887
888        NOTE:
889            Despite the fact that higher order terms are printed to the
890            right in a power series, right shifting decreases the powers
891            of $t$, while left shifting increases them. This is to be
892            consistant with polynomials, integers, etc.
893
894        EXAMPLES:
895            sage: R.<t> = PowerSeriesRing(QQ['y'], 't', 5)
896            sage: f = ~(1+t); f
897            1 + -t + t^2 + -t^3 + t^4 + O(t^5)
898            sage: f.shift(3)
899            t^3 + -t^4 + t^5 + -t^6 + t^7 + O(t^8)
900            sage: f >> 2
901            1 + -t + t^2 + O(t^3)
902            sage: f << 10
903            t^10 + -t^11 + t^12 + -t^13 + t^14 + O(t^15)
904            sage: t << 29
905            t^30
906
907        AUTHOR:
909        """
910        return self._parent(self.polynomial().shift(n), self._prec + n)
911
912    def __lshift__(self, n):
913        return self.parent()(self.polynomial() << n, self.prec())
914
915    def __rshift__(self, n):
916        return self.parent()(self.polynomial() >> n, self.prec())
917
918    def is_square(self):
919        """
920        Returns True if this function has a square root in this ring,
921        e.g. there is an element $y$ in \code{self.parent()} such that
922        $y^2 = \code{self}$.
923
924        ALGORITHM:
925            If the basering is a field, this is true whenver the power
926            series has even valuation and the leading coefficent is a
927            perfect square.
928
929            For an integral domain, it operates attempts the square root
930            in the fraction field and tests whether or not the result
931            lies in the original ring.
932
933        EXAMPLES:
934            sage: K.<t> = PowerSeriesRing(QQ, 't', 5)
935            sage: (1+t).is_square()
936            True
937            sage: (2+t).is_square()
938            False
939            sage: (2+t.change_ring(RR)).is_square()
940            True
941            sage: t.is_square()
942            False
943            sage: K.<t> = PowerSeriesRing(ZZ, 't', 5)
944            sage: (1+t).is_square()
945            False
946            sage: f = (1+t)^100
947            sage: f.is_square()
948            True
949
950        """
951        val = self.valuation()
952        if val is not infinity and val % 2 == 1:
953            return False
954        elif not self[val].is_square():
955            return False
956        elif is_Field(self.base_ring()):
957            return True
958        else:
959            try:
960                self.parent()(self.sqrt())
961                return True
962            except TypeError:
963                return False
964
965    def sqrt(self, prec=None, extend=False,
966             all=False, name=None):
967        r"""
968        The square root function.
969
970        INPUT:
971
972            prec -- integer (default: None): if not None and the series
973                 has infinite precision, truncates series at precision prec.
974            extend -- bool (default: False); if True, return a square
975                 root in an extension ring, if necessary. Otherwise,
976                 raise a ValueError if the square is not in the base
977                 power series ring.  For example, if extend is True
978                 the square root of a power series with odd degree
979                 leading coefficient is defined as an element of a formal
980                 extension ring.
981            name -- if extend is True, you must also specify the print name of
982                 formal square root.
983            all -- bool (default: False); if True, return all square
984                 roots of self, instead of just one.
985
986        ALGORITHM:
987            Newton's method
988            $$989 x_{i+1} = \frac{1}{2}( x_i + self/x_i ) 990$$
991
992        EXAMPLES:
993            sage: K.<t> = PowerSeriesRing(QQ, 't', 5)
994            sage: sqrt(t^2)
995            t
996            sage: sqrt(1+t)
997            1 + 1/2*t - 1/8*t^2 + 1/16*t^3 - 5/128*t^4 + O(t^5)
998            sage: sqrt(4+t)
999            2 + 1/4*t - 1/64*t^2 + 1/512*t^3 - 5/16384*t^4 + O(t^5)
1000            sage: u = sqrt(2+t, prec=2, extend=True, name = 'alpha'); u
1001            alpha
1002            sage: u^2
1003            2 + t
1004            sage: u.parent()
1005            Univariate Quotient Polynomial Ring in alpha over Power Series Ring in t over Rational Field with modulus x^2 + -2 - t
1006            sage: K.<t> = PowerSeriesRing(QQ, 't', 50)
1007            sage: sqrt(1+2*t+t^2)
1008            1 + t + O(t^50)
1009            sage: sqrt(t^2 +2*t^4 + t^6)
1010            t + t^3 + O(t^51)
1011            sage: sqrt(1 + t + t^2 + 7*t^3)^2
1012            1 + t + t^2 + 7*t^3 + O(t^50)
1013            sage: sqrt(K(0))
1014            0
1015            sage: sqrt(t^2)
1016            t
1017
1018            sage: K.<t> = PowerSeriesRing(CDF, 5)
1019            sage: v = sqrt(-1 + t + t^3, all=True); v
1020            [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),
1021             -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)]
1022            sage: [a^2 for a in v]
1023            [-1.0 + 1.0*t + 1.0*t^3 + O(t^5), -1.0 + 1.0*t + 1.0*t^3 + O(t^5)]
1024
1025        A formal square root:
1026            sage: K.<t> = PowerSeriesRing(QQ, 5)
1027            sage: f = 2*t + t^3 + O(t^4)
1028            sage: s = f.sqrt(extend=True, name='sqrtf'); s
1029            sqrtf
1030            sage: s^2
1031            2*t + t^3 + O(t^4)
1032            sage: parent(s)
1033            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)
1034
1035        AUTHORS:
1037            -- William Stein
1038        """
1039        if self.is_zero():
1040            ans = self._parent(0).O(self.prec()/2)
1041            if all:
1042                return [ans]
1043            else:
1044                return ans
1045
1046        if all and not self.base_ring().is_integral_domain():
1047            raise NotImplementedError, 'all roots not implemented over a non-integral domain'
1048
1049        formal_sqrt = False
1050        u = self.valuation_zero_part()
1051        # TODO, fix underlying element sqrt()
1052        try:
1053            try:
1054                s = u[0].sqrt(extend=False)
1055            except TypeError:
1056                s = u[0].sqrt()
1057        except ValueError:
1058            formal_sqrt = True
1059        if self.degree() == 0:
1060            if not formal_sqrt:
1061                a = self.parent()([s], self.prec())
1062                if all:
1063                    return [a, -a]
1064                else:
1065                    return a
1066
1067        val = self.valuation()
1068        if formal_sqrt or val % 2 == 1:
1069            if extend:
1070                if name is None:
1071                    raise ValueError, "the square root generates an extension, so you must specify the name of the square root"
1072                R = self._parent['x']
1073                S = R.quotient(R([-self,0,1]), names=name)
1074                a = S.gen()
1075                if all:
1076                    if not self.base_ring().is_integral_domain():
1077                        raise NotImplementedError, 'all roots not implemented over a non-integral domain'
1078                    return [a, -a]
1079                else:
1080                    return a
1081            else:
1082                raise ValueError, "power series does not have a square root since it has odd valuation."
1083
1084        pr = self.prec()
1085        if pr == infinity:
1086            if prec is None:
1087                pr = self._parent.default_prec()
1088            else:
1089                pr = prec
1090        prec = pr
1091
1092        R = s.parent()
1093        a = self.valuation_zero_part()
1094        P = self._parent
1095        if not R is P.base_ring():
1096            a = a.change_ring(R)
1097        half = ~R(2)
1098
1099        s = a.parent()([s])
1100        for cur_prec in sage.misc.misc.newton_method_sizes(prec)[1:]:
1101            (<PowerSeries>s)._prec = cur_prec
1102            s = half * (s + a/s)
1103
1104        ans = s
1105        if val != 0:
1106            ans *= P.gen(0)**(val/2)
1107
1108        if all:
1109            return [ans, -ans]  # since over an integral domain
1110        else:
1111            return ans
1112
1113    def square_root(self):
1114        """
1115        Return the square root of self in this ring. If this cannot be done
1116        then an error will be raised.
1117
1118        This function succeeds if and only if \code{self.is_square()}
1119
1120        EXAMPLES:
1121            sage: K.<t> = PowerSeriesRing(QQ, 't', 5)
1122            sage: (1+t).square_root()
1123            1 + 1/2*t - 1/8*t^2 + 1/16*t^3 - 5/128*t^4 + O(t^5)
1124            sage: (2+t).square_root()
1125            Traceback (most recent call last):
1126            ...
1127            ValueError: Square root does not live in this ring.
1128            sage: (2+t.change_ring(RR)).square_root()
1129            1.41421356237309 + 0.353553390593274*t - 0.0441941738241592*t^2 + 0.0110485434560399*t^3 - 0.00345266983001233*t^4 + O(t^5)
1130            sage: t.square_root()
1131            Traceback (most recent call last):
1132            ...
1133            ValueError: Square root not defined for power series of odd valuation.
1134            sage: K.<t> = PowerSeriesRing(ZZ, 't', 5)
1135            sage: f = (1+t)^20
1136            sage: f.square_root()
1137            1 + 10*t + 45*t^2 + 120*t^3 + 210*t^4 + O(t^5)
1138            sage: f = 1+t
1139            sage: f.square_root()
1140            Traceback (most recent call last):
1141            ...
1142            ValueError: Square root does not live in this ring.
1143
1144        AUTHOR:
1146        """
1147        val = self.valuation()
1148        if val is not infinity and val % 2 == 1:
1149            raise ValueError, "Square root not defined for power series of odd valuation."
1150        elif not self[val].is_square():
1151            raise ValueError, "Square root does not live in this ring."
1152        elif is_Field(self.base_ring()):
1153            return self.sqrt()
1154        else:
1155            try:
1156                return self.parent()(self.sqrt())
1157            except TypeError:
1158                raise ValueError, "Square root does not live in this ring."
1159
1160    def O(self, prec):
1161        r"""
1162        Return this series plus $O(x^\text{prec})$.  Does not change
1163        self.
1164        """
1165        if prec is infinity or prec >= self.prec():
1166            return self
1167        coeffs = self[:prec]
1168        return self._parent(coeffs, prec)
1169
1170
1171    def solve_linear_de(self, prec = infinity, b = None, f0 = None):
1172        r"""
1173        Obtains a power series solution to an inhomogeneous linear
1174        differential equation of the form:
1175           $$f'(t) = a(t) f(t) + b(t).$$
1176
1177        INPUT:
1178            self -- the power series $a(t)$
1179            b -- the power series $b(t)$ (default is zero)
1180            f0 -- the constant term of $f$ (initial condition'')
1181                 (default is 1)
1182            prec -- desired precision of result (this will be reduced if
1183                    either a or b have less precision available)
1184
1185        OUTPUT:
1186            the power series f, to indicated precision
1187
1188        ALGORITHM:
1189            A divide-and-conquer strategy; see the source code. Running time
1190            is approximately $M(n) \log n$, where $M(n)$ is the time required
1191            for a polynomial multiplication of length $n$ over the coefficient
1192            ring. (If you're working over something like RationalField(),
1193            running time analysis can be a little complicated because the
1194            coefficients tend to explode.)
1195
1196        NOTES:
1197            -- If the coefficient ring is a field of characteristic zero,
1198               then the solution will exist and is unique.
1199            -- For other coefficient rings, things are more complicated.
1200               A solution may not exist, and if it does it may not be unique.
1201               Generally, by the time the nth term has been computed, the
1202               algorithm will have attempted divisions by $n!$ in the
1203               coefficient ring. So if your coefficient ring has enough
1204               precision'', and if your coefficient ring can perform
1205               divisions even when the answer is not unique, and if you know
1206               in advance that a solution exists, then this function will find
1207               a solution (otherwise it will probably crash).
1208
1209        AUTHORS:
1210          -- David Harvey (2006-09-11): factored functionality out from
1211             exp() function, cleaned up precision tests a bit
1212
1213        EXAMPLES:
1214           sage: R.<t> = PowerSeriesRing(QQ, default_prec=10)
1215
1216           sage: a = 2 - 3*t + 4*t^2 + O(t^10)
1217           sage: b = 3 - 4*t^2 + O(t^7)
1218           sage: f = a.solve_linear_de(prec=5, b=b, f0=3/5)
1219           sage: f
1220            3/5 + 21/5*t + 33/10*t^2 - 38/15*t^3 + 11/24*t^4 + O(t^5)
1221           sage: f.derivative() - a*f - b
1222            O(t^4)
1223
1224           sage: a = 2 - 3*t + 4*t^2
1225           sage: b = b = 3 - 4*t^2
1226           sage: f = a.solve_linear_de(b=b, f0=3/5)
1227           Traceback (most recent call last):
1228           ...
1229           ValueError: cannot solve differential equation to infinite precision
1230
1231           sage: a.solve_linear_de(prec=5, b=b, f0=3/5)
1232            3/5 + 21/5*t + 33/10*t^2 - 38/15*t^3 + 11/24*t^4 + O(t^5)
1233        """
1234        if b is None:
1235            b = self._parent(0)
1236
1237        if f0 is None:
1238            f0 = 1
1239        f0 = self.base_ring()(f0)
1240
1241        # reduce precision to whatever is available from a and b
1242        prec = min(prec, self.prec() + 1, b.prec() + 1)
1243        if prec is infinity:
1244            raise ValueError, \
1245                 "cannot solve differential equation to infinite precision"
1246        prec = int(prec)
1247        if prec == 0:
1248            return self._parent(0, 0)
1249        if prec < 0:
1250            raise ValueError, \
1251                  "prec (=%s) must be a non-negative integer" % prec
1252
1253        base_ring = self._parent.base_ring()
1254        R = PolynomialRing(base_ring, self._parent.variable_name())
1255
1256        a_list = self.list()
1257        b_list = [base_ring(0)]
1258        b_list.extend(b.list())
1259        while len(b_list) < prec:
1260            b_list.append(base_ring(0))
1261        f = _solve_linear_de(R, 0, prec, a_list, b_list, f0)
1262        return self._parent(f, prec)
1263
1264
1265    def exp(self, prec=None):
1266        r"""
1267        Returns exp of this power series to the indicated precision.
1268
1269        INPUT:
1270            prec -- integer; default is self.parent().default_prec
1271
1272        ALGORITHM:
1273            See PowerSeries.solve_linear_de().
1274
1275        NOTES:
1276            -- Screwy things can happen if the coefficient ring is not
1277               a field of characteristic zero.
1278               See PowerSeries.solve_linear_de().
1279
1280        AUTHORS:
1281          -- David Harvey (2006-09-08): rewrote to use simplest possible
1282             lazy'' algorithm.
1283          -- David Harvey (2006-09-10): rewrote to use divide-and-conquer
1284             strategy.
1285          -- David Harvey (2006-09-11): factored functionality out to
1286             solve_linear_de().
1287
1288        EXAMPLES:
1289           sage: R.<t> = PowerSeriesRing(QQ, default_prec=10)
1290
1291         Check that $\exp(t)$ is, well, $\exp(t)$:
1292           sage: (t + O(t^10)).exp()
1293           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)
1294
1295         Check that $\exp(\log(1+t))$ is $1+t$:
1296           sage: (sum([-(-t)^n/n for n in range(1, 10)]) + O(t^10)).exp()
1297           1 + t + O(t^10)
1298
1299         Check that $\exp(2t + t^2 - t^5)$ is whatever it is:
1300           sage: (2*t + t^2 - t^5 + O(t^10)).exp()
1301           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)
1302
1303         Check requesting lower precision:
1304           sage: (t + t^2 - t^5 + O(t^10)).exp(5)
1305           1 + t + 3/2*t^2 + 7/6*t^3 + 25/24*t^4 + O(t^5)
1306
1307         Can't get more precision than the input:
1308           sage: (t + t^2 + O(t^3)).exp(10)
1309           1 + t + 3/2*t^2 + O(t^3)
1310
1311         Check some boundary cases:
1312           sage: (t + O(t^2)).exp(1)
1313           1 + O(t)
1314           sage: (t + O(t^2)).exp(0)
1315           O(t^0)
1316
1317        """
1318        if prec is None:
1319            prec = self._parent.default_prec()
1320        if not self[0].is_zero():
1321            raise ValueError, "constant term must to zero"
1322        return self.derivative().solve_linear_de(prec)
1323
1324
1325    def V(self, n):
1326        """
1327        If $f = \sum a_m x^m$, then this function returns $\sum a_m x^{nm}$.
1328        """
1329        v = self.list()
1330        m = 0
1331        w = []
1332        zero = self.base_ring()(0)
1333        for i in range(len(v)*n):
1334            if i%n != 0:
1335                w.append(zero)
1336            else:
1337                w.append(v[m])
1338                m += 1
1339        return self._parent(w, self.prec()*n)
1340
1341    def valuation(self):
1342        """
1343        Return the valuation of this power series.
1344
1345        This is equal to the valuation of the underlying polynomial.
1346
1347        EXAMPLES:
1348        Sparse examples:
1349            sage: R.<t> = PowerSeriesRing(QQ, sparse=True)
1350            sage: f = t^100000 + O(t^10000000)
1351            sage: f.valuation()
1352            100000
1353            sage: R(0).valuation()
1354            +Infinity
1355
1356        Dense examples:
1357            sage: R.<t> = PowerSeriesRing(ZZ)
1358            sage: f = 17*t^100 +O(t^110)
1359            sage: f.valuation()
1360            100
1361            sage: t.valuation()
1362            1
1363        """
1364        return self.polynomial().valuation()
1365
1366    def variable(self):
1367        """
1368        EXAMPLES:
1369            sage: R.<x> = PowerSeriesRing(Rationals())
1370            sage: f = x^2 + 3*x^4 + O(x^7)
1371            sage: f.variable()
1372            'x'
1373
1374        AUTHOR:
1375            -- David Harvey (2006-08-08)
1376        """
1377        return self._parent.variable_name()
1378
1379    def degree(self):
1380        """
1381        Return the degree of this power series, which is by definition
1382        the degree of the underlying polynomial.
1383
1384        EXAMPLES:
1385            sage: R.<t> = PowerSeriesRing(QQ, sparse=True)
1386            sage: f = t^100000 + O(t^10000000)
1387            sage: f.degree()
1388            100000
1389        """
1390        return self.polynomial().degree()
1391
1392    def derivative(self):
1393        raise NotImplementedError
1394
1395#def _normalize(v, prec):
1396#    n = len(v)-1
1397#    while n>=0 and (v[n] == 0 or (prec != infinity and n >= prec)):
1398#        del v[n]
1399#        n -= 1
1400
1401cdef class PowerSeries_poly(PowerSeries):
1402
1403    def __init__(self, parent, f=0, prec=infinity, int check=1, is_gen=0):
1404        """
1405        EXAMPLES:
1406            sage: R, q = PowerSeriesRing(CC, 'q').objgen()
1407            sage: R
1408            Power Series Ring in q over Complex Field with 53 bits of precision
1410            True
1411        """
1412        R = parent._poly_ring()
1413        if PY_TYPE_CHECK(f, Element):
1414            if (<Element>f)._parent is R:
1415                pass
1416            elif (<Element>f)._parent is R.base_ring():
1417                f = R([f])
1418            elif PY_TYPE_CHECK(f, PowerSeries_poly):
1419                prec = (<PowerSeries_poly>f)._prec
1420                f = R((<PowerSeries_poly>f).__f)
1421            else:
1422                f = R(f, check=check)
1423        else:
1424            f = R(f, check=check)
1425
1426        self.__f = f
1427        if check and not (prec is infinity):
1428            self.__f = self.__f.truncate(prec)
1429        PowerSeries.__init__(self, parent, prec, is_gen)
1430
1431    def __hash__(self):
1432        return hash(self.__f)
1433
1434    def __reduce__(self):
1435        # do *not* delete old versions.
1436        return make_powerseries_poly_v0, (self._parent, self.__f, self._prec, self.__is_gen)
1437
1438    def __richcmp__(left, right, int op):
1439        return (<Element>left)._richcmp(right, op)
1440
1441    def __pow__(self_t, r, dummy):
1442        cdef PowerSeries_poly self = self_t
1443        cdef int right = r
1444        if right != r:
1445            raise ValueError, "exponent must be an integer"
1446        if right < 0:
1447            return (~self)**(-right)
1448        if right == 0:
1449            return self._parent(1)
1450        if self.__is_gen:
1451            return PowerSeries_poly(self._parent, self.__f**right, check=False)
1452        if self.is_zero():
1453            return self
1454        return arith.generic_power(self, right, self._parent(1))
1455
1456    def polynomial(self):
1457        """
1458        EXAMPLE:
1459            sage: R.<t> = GF(7)[[]]
1460            sage: f = 3 - t^3 + O(t^5)
1461            sage: f.polynomial()
1462            6*t^3 + 3
1463        """
1464        return self.__f
1465
1466    def valuation(self):
1467        return self.__f.valuation()
1468
1469    def degree(self):
1470        return self.__f.degree()
1471
1472    def __call__(self, *xs):
1473        """
1474        EXAMPLE:
1475            sage: R.<t> = GF(7)[[]]
1476            sage: f = 3 - t^3 + O(t^5)
1477            sage: f(1)
1478            2
1479        """
1480        if isinstance(xs[0], tuple):
1481            xs = xs[0]
1482        x = xs[0]
1483        try:
1484            if x.parent() is self._parent:
1485                if not (self.prec() is infinity):
1487                    xs = list(xs); xs[0] = x; xs = tuple(xs) # tuples are immutable
1488        except AttributeError:
1489            pass
1490        return self.__f(xs)
1491
1492    def __setitem__(self, n, value):
1493        """
1494        EXAMPLES:
1495            sage: R.<t> = ZZ[[]]
1496            sage: f = 3 - t^3 + O(t^5)
1497            sage: f[1] = 5
1498            Traceback (most recent call last):
1499            ...
1500            IndexError: power series are immutable
1501        """
1502        raise IndexError, "power series are immutable"
1503
1504    def __getslice__(self, i, j):
1505        r"""
1506        Return slice of coefficient of this power series.
1507
1508        This calls slice on the underlying polynomial, and makes a power
1509        series out of the result, with precision the precision of self.
1510
1511        EXAMPLES:
1512            sage: R.<t> = ZZ[[]]
1513            sage: f = (2-t)^5 + O(t^7); f
1514            32 - 80*t + 80*t^2 - 40*t^3 + 10*t^4 - t^5 + O(t^7)
1515            sage: f[2:4]
1516            80*t^2 - 40*t^3 + O(t^7)
1517        """
1518        return PowerSeries_poly(self._parent, self.__f[i:j], prec=self.prec(), check=False)
1519
1520    def _unsafe_mutate(self, i, value):
1521        """
1522        SAGE assumes throughout that commutative ring elements are immutable.
1523        This is relevant for caching, etc.  But sometimes you need to change
1524        a power series and you really know what you're doing.  That's
1525        when this function is for you.
1526
1527        ** DO NOT USE THIS ** unless you know what you're doing.
1528
1529        EXAMPLES:
1530            sage: R.<t> = GF(7)[[]]
1531            sage: f = 3 + 6*t^3 + O(t^5)
1532            sage: f._unsafe_mutate(0, 5)
1533            sage: f
1534            5 + 6*t^3 + O(t^5)
1535
1536        Mutating can even bump up the precision.
1537            sage: f._unsafe_mutate(7,2)
1538            sage: f
1539            5 + 6*t^3 + 2*t^7 + O(t^8)
1540        """
1541        self.__f._unsafe_mutate(i, value)
1542        self._prec = max(self._prec, i+1)
1543
1544    def __getitem__(self, n):
1545        """
1546        Return the n-th coefficient.
1547
1548        Returns 0 for negative coefficients.  Raises an IndexError if
1549        try to access beyond known coefficients.
1550
1551        EXAMPLES:
1552            sage: R.<t> = QQ[[]]
1553            sage: f = 3/2 - 17/5*t^3 + O(t^5)
1554            sage: f[3]
1555            -17/5
1556            sage: f[-2]
1557            0
1558            sage: f[4]
1559            0
1560            sage: f[5]
1561            Traceback (most recent call last):
1562            ...
1563            IndexError: coefficient not known
1564            sage: f[1:4]
1565            -17/5*t^3 + O(t^5)
1566        """
1567        if n<0:
1568            return self.base_ring()(0)
1569        if n > self.__f.degree():
1570            if self._prec > n:
1571                return self.base_ring()(0)
1572            #elif isinstance(n, slice):
1573                # It makes no sense that this is needed and that
1574                # __getslice__ isn't just called by default...
1575            #    return self.__getslice__(slice[0],slice[1])
1576            else:
1577                raise IndexError, "coefficient not known"
1578        return self.__f[n]
1579
1580    def __iter__(self):
1581        """
1582        Return an interator over the coefficients of this power series.
1583
1584        EXAMPLES:
1585            sage: R.<t> = QQ[[]]
1586            sage: f = t + 17/5*t^3 + 2*t^4 + O(t^5)
1587            sage: for a in f: print a,
1588            0 1 0 17/5 2
1589        """
1590        return iter(self.__f)
1591
1592    def __neg__(self):
1593        """
1594        Return the negative of this power series.
1595
1596        EXAMPLES:
1597            sage: R.<t> = QQ[[]]
1598            sage: f = t + 17/5*t^3 + 2*t^4 + O(t^5)
1599            sage: -f
1600            -t - 17/5*t^3 - 2*t^4 + O(t^5)
1601        """
1602        return PowerSeries_poly(self._parent, -self.__f,
1603                                         self._prec, check=False)
1604
1605    cdef ModuleElement _add_c_impl(self, ModuleElement right_m):
1606        """
1607        EXAMPLES:
1608            sage: R.<x> = PowerSeriesRing(ZZ)
1609            sage: f = x^4 + O(x^5); f
1610            x^4 + O(x^5)
1611            sage: g = x^2 + O(x^3); g
1612            x^2 + O(x^3)
1613            sage: f+g
1614            x^2 + O(x^3)
1615        """
1616        cdef PowerSeries_poly right = <PowerSeries_poly>right_m
1617        return PowerSeries_poly(self._parent, self.__f + right.__f, \
1618                                         self.common_prec_c(right), check=True)
1619
1620    cdef ModuleElement _sub_c_impl(self, ModuleElement right_m):
1621        """
1622        Return difference of two power series.
1623
1624        EXAMPLES:
1625            sage: k.<w> = ZZ[]
1626            sage: R.<t> = k[[]]
1627            sage: w*t^2 -w*t +13 - (w*t^2 + w*t)
1628            13 + -2*w*t
1629        """
1630        cdef PowerSeries_poly right = <PowerSeries_poly>right_m
1631        return PowerSeries_poly(self._parent, self.__f - right.__f, \
1632                                         self.common_prec_c(right), check=True)
1633
1634    cdef RingElement _mul_c_impl(self, RingElement right_r):
1635        """
1636        Return the product of two power series.
1637
1638        EXAMPLES:
1639            sage: k.<w> = ZZ[[]]
1640            sage: (1+17*w+15*w^3+O(w^5))*(19*w^10+O(w^12))
1641            19*w^10 + 323*w^11 + O(w^12)
1642        """
1643        cdef PowerSeries_poly right = <PowerSeries_poly>right_r
1644        sp = self._prec
1645        rp = right._prec
1646        if is_Infinite(sp):
1647            if is_Infinite(rp):
1648                prec = infinity
1649            else:
1650                prec = rp + self.valuation()
1651        else:  # sp != infinity
1652            if is_Infinite(rp):
1653                prec = sp + right.valuation()
1654            else:
1655                prec = min(rp + self.valuation(), sp + right.valuation())
1656        return PowerSeries_poly(self._parent,
1657                                         self.__f * right.__f,
1658                                         prec,
1659                                         check=True)  # check, since truncation may be needed
1660
1661    def _rmul_(self, c):
1662        return PowerSeries_poly(self._parent, c * self.__f, self._prec, check=False)
1663
1664    def _lmul_(self, c):
1665        return PowerSeries_poly(self._parent, self.__f * c, self._prec, check=False)
1666
1667
1668    def __floordiv__(self, denom):
1669        try:
1670            return PowerSeries.__div__(self, denom)
1671        except (PariError, ZeroDivisionError), e: # PariError to general?
1672            if is_PowerSeries(denom) and denom.degree() == 0 and denom[0] in self._parent.base_ring():
1673                denom = denom[0]
1674            elif not denom in self._parent.base_ring():
1675                raise ZeroDivisionError, e
1676            return PowerSeries_poly(self._parent,
1677                                             self.__f // denom, self._prec)
1678
1679    def __lshift__(_self, n):
1680        cdef PowerSeries_poly self = _self
1681        return PowerSeries_poly(self._parent, self.__f << n, self._prec + n)
1682
1683    def __rshift__(_self, n):
1684        cdef PowerSeries_poly self = _self
1685        return PowerSeries_poly(self._parent, self.__f >> n, self._prec - n)
1686
1687    def copy(self):
1688        return PowerSeries_poly(self._parent, self.__f, self.prec(), check=False)
1689
1690    def list(self):
1691        return self.__f.list()
1692
1693    def dict(self):
1694        return self.__f.dict()
1695
1696    def derivative(self):
1697        """
1698        Return the derivative of this power series.
1699
1700        EXAMPLES:
1701            sage: R.<t> = PowerSeriesRing(QQ, sparse=True)
1702            sage: f = 2 + 3*t^2 + t^100000 + O(t^10000000); f
1703            2 + 3*t^2 + t^100000 + O(t^10000000)
1704            sage: f.derivative()
1705            6*t + 100000*t^99999 + O(t^9999999)
1706        """
1707        return PowerSeries_poly(self._parent, self.__f.derivative(),
1708                                         self.prec()-1, check=False)
1709
1710    def integral(self):
1711        """
1712        The integral of this power series with 0 constant term.
1713
1714        EXAMPLES:
1715            sage: k.<w> = QQ[[]]
1716            sage: (1+17*w+15*w^3+O(w^5)).integral()
1717            w + 17/2*w^2 + 15/4*w^4 + O(w^6)
1718        """
1719        return PowerSeries_poly(self._parent, self.__f.integral(),
1720                                         self.prec()+1, check=False)
1721
1722    def laurent_series(self):
1723        """
1724        Return the Laurent series associated to this power series, i.e., this
1725        series considered as a Laurent series.
1726
1727        EXAMPLES:
1728            sage: k.<w> = QQ[[]]
1729            sage: f = 1+17*w+15*w^3+O(w^5)
1730            sage: parent(f)
1731            Power Series Ring in w over Rational Field
1732            sage: g = f.laurent_series(); g
1733            1 + 17*w + 15*w^3 + O(w^5)
1734        """
1735        return self._parent.laurent_series_ring()(self)
1736
1737    def ogf(self):
1738        r"""
1739        Returns the ordinary generating function associated to self.
1740
1741        This function is known as serlaplace in GP/PARI.
1742
1743        EXAMPLES:
1744            sage: R.<t> = PowerSeriesRing(QQ)
1745            sage: f = t + t^2/factorial(2) + 2*t^3/factorial(3)
1746            sage: f.ogf()
1747            t + t^2 + 2*t^3
1748        """
1749        return self.parent()([self[i] * arith.factorial(i) for i in range(self.degree()+1)])
1750
1751    def egf(self):
1752        r"""
1753        Returns the exponential generating function associated to self.
1754
1755        This function is known as serlaplace in GP/PARI.
1756
1757        EXAMPLES:
1758            sage: R.<t> = PowerSeriesRing(QQ)
1759            sage: f = t + t^2 + 2*t^3
1760            sage: f.egf()
1761            t + 1/2*t^2 + 1/3*t^3
1762        """
1763        return self.parent()([self[i] / arith.factorial(i) for i in range(self.degree()+1)])
1764
1765    def reversion(self):
1766        """
1767        Return the reversion of f, i.e., the series g such that
1768        g(f(x)) = x.
1769
1770        EXAMPLES:
1771            sage: R.<x> = PowerSeriesRing(QQ)
1772            sage: f = 2*x + 3*x**2 - x**4 + O(x**5)
1773            sage: g = f.reversion()
1774            sage: g
1775            1/2*x - 3/8*x^2 + 9/16*x^3 - 131/128*x^4 + O(x^5)
1776            sage: f(g)
1777            x + O(x^5)
1778            sage: g(f)
1779            x + O(x^5)
1780        """
1781        if not isinstance(self.parent().base_ring(), rational_field.RationalField):
1782            raise NotImplementedError
1783        if self.prec() is infinity:
1784            raise RuntimeError, "series must have finite precision for reversion."
1785        f = self._pari_()
1786        g = f.serreverse()
1787        return PowerSeries_poly(self.parent(),g.Vecrev(),self.prec())
1788
1789    def _pari_(self):
1790        """
1791        Return PARI power series corresponding to this series.
1792
1793        This is currently only implemented over QQ and ZZ.
1794
1795        EXAMPLES:
1796            sage: k.<w> = QQ[[]]
1797            sage: f = 1+17*w+15*w^3+O(w^5)
1798            sage: pari(f)
1799            1 + 17*w + 15*w^3 + O(w^5)
1800            sage: pari(1 - 19*w + w^5)
1801            Traceback (most recent call last):
1802            ...
1803            RuntimeError: series precision must be finite for conversion to pari object.
1804        """
1805        if not isinstance(self.parent().base_ring(),
1806                          (rational_field.RationalField, integer_ring.IntegerRing)):
1807            raise NotImplementedError
1808        if self.prec() is infinity:
1809            raise RuntimeError, "series precision must be finite for conversion to pari object."
1810        return sage.libs.pari.all.pari(str(self))
1811
1812
1813
1814def _solve_linear_de(R, N, L, a, b, f0):
1815    r"""
1816    Internal function used by PowerSeries.solve_linear_de().
1817
1818    INPUT:
1819       R -- a PolynomialRing
1820       N -- integer >= 0
1821       L -- integer >= 1
1822       a -- list of coefficients of $a$, any length, all coefficients should
1823            belong to base ring of R.
1824       b -- list of coefficients of $b$, length at least $L$
1825            (only first $L$ coefficients are used), all coefficients
1826            should belong to base ring of R.
1827       f0 -- constant term of $f$ (only used if $N == 0$), should belong
1828            to base ring of R.
1829
1830    OUTPUT:
1831       List of coefficients of $f$ (length exactly $L$), where $f$ is a
1832       solution to the linear inhomogeneous differential equation:
1833          $$(t^N f)' = a t^N f + t^{N-1} b + O(t^{N+L-1}).$$
1834       The constant term of $f$ is taken to be f0 if $N == 0$, and otherwise
1835       is determined by $N f_0 = b_0$.
1836
1837    ALGORITHM:
1838        The algorithm implemented here is inspired by the fast lazy
1839        multiplication algorithm'' described in the paper Lazy multiplication
1840        of formal power series'' by J van der Hoeven (1997 ISSAC conference).
1841
1842        Our situation is a bit simpler than the one described there,
1843        because only one of the series being multiplied needs the lazy
1844        treatment; the other one is known fully in advance. I have
1845        reformulated the algorithm as an explicit divide-and-conquer
1846        recursion. Possibly it is slower than the one described by van der
1847        Hoeven, by a constant factor, but it seems easier to implement.
1848        Also, because we repeatedly split in half starting at the top level,
1849        instead of repeatedly doubling in size from the bottom level, it's
1850        easier to avoid problems with overshooting'' in the last iteration.
1851
1852        The idea is to split the problem into two instances with $L$
1853        about half the size. Take $L'$ to be the ceiling of
1854        $(L/2)$. First recursively find $g$ modulo $t^{L'}$ such that
1855          $$(t^N g)' = a t^N g + t^{N-1} b + O(t^{N+L'-1}).$$
1856
1857        Next we want to find $h$ modulo $t^{L-L'}$ such that
1858        $f = g + t^{L'} h$ is a solution of the original equation.
1859        We can find a suitable $h$ by recursively solving another
1860        differential equation of the same form, namely
1861          $$(t^{N+L'} h)' = a t^{N+L'} h + c t^{N+L'-1} + O(t^{N+L'-1}),$$
1862        where $c$ is determined by
1863          $$(t^N g)' - a t^N g - t^{N-1} b = -c t^{N+L'-1} + O(t^{N+L-1}).$$
1864        Once $g$ is known, $c$ can be recovered from this relation by computing
1865        the coefficients of $t^j$ of the product $a g$ for $j$ in the range
1866        $L'-1 \leq j < L-1$.
1867
1868        At the bottom of the recursion we have $L = 1$, in which case
1869        the solution is simply given by $f_0 = b_0/N$ (or by the supplied
1870        initial condition $f_0$ when $N == 0$).
1871
1872        The algorithm has to do one multiplication of length $N$, two of
1873        length $N/2$, four of length $N/4$, etc, hence giving runtime
1874        $O(M(N) \log N)$.
1875
1876    AUTHOR:
1877        -- David Harvey (2006-09-11): factored out of PowerSeries.exp().
1878
1879    """
1880    if L == 1:
1881        # base case
1882        if N == 0:
1883            return [f0]
1884        else:
1885            return [b[0] / N]
1886
1887    L2 = (L + 1) >> 1    # ceil(L/2)
1888
1889    g = _solve_linear_de(R, N, L2, a, b, f0)
1890
1891    term1 = R(g, check=False)
1892    term2 = R(a[:L], check=False)
1893    product = (term1 * term2).list()
1894
1895    # todo: perhaps next loop could be made more efficient
1896    c = b[L2 : L]
1897    for j in range(L2 - 1, min(L-1, len(product))):
1898        c[j - (L2-1)] = c[j - (L2-1)] + product[j]
1899
1900    h = _solve_linear_de(R, N + L2, L - L2, a, c, f0)
1901
1902    g.extend(h)
1903    return g
1904
1905
1906def make_powerseries_poly_v0(parent,  f, prec, is_gen):
1907    return PowerSeries_poly(parent, f, prec, 0, is_gen)
1908
1909def make_element_from_parent_v0(parent, *args):
1910    return parent(*args)
Note: See TracBrowser for help on using the repository browser.