Ticket #6800: 12846.patch

File 12846.patch, 108.0 KB (added by Henryk.Trappmann, 11 years ago)

patch adds the file formal_powerseries.py

  • new file sage/rings/formal_powerseries.py

    # HG changeset patch
    # User Henryk Trappmann <bo198214@googlemail.com>
    # Date 1250925673 -7200
    # Node ID 172e984ca65257ee048e6d6489dcb2dfb5ccf5c6
    # Parent  684eea91ff224e5bc6259ca19f1576c4c082b9d3
    added formal_powerseries.py
    
    diff -r 684eea91ff22 -r 172e984ca652 sage/rings/formal_powerseries.py
    - +  
     1"""
     2Cached formal powerseries and formal Laurant series in one variable.
     3
     4Author: Henryk Trappmann
     5"""
     6
     7from sage.structure.sage_object import SageObject
     8from sage.rings.arith import factorial
     9from sage.rings.arith import binomial
     10from sage.rings.integer import Integer
     11from sage.calculus.calculus import var
     12from sage.symbolic.expression import Expression
     13from sage.calculus.functional import diff
     14from sage.rings.ring import Ring
     15from sage.rings.ring_element import RingElement
     16from sage.rings.rational_field import QQ, RationalField
     17from sage.rings.rational import Rational
     18from sage.rings.real_mpfr import RR, RealField
     19from sage.rings.complex_field import ComplexField_class
     20from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
     21from sage.rings.polynomial.polynomial_ring import PolynomialRing_field
     22from sage.misc.misc_c import prod
     23from sage.misc.functional import n
     24from sage.rings.infinity import Infinity
     25from sage.rings.power_series_ring_element import PowerSeries
     26from sage.rings.polynomial.polynomial_element import Polynomial
     27from sage.matrix.constructor import matrix
     28
     29def decidable0(K):
     30    """
     31    Returns true if K has a decidable 0.
     32    For example powerseries have no decidable 0.
     33
     34    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing,decidable0
     35    sage: decidable0(QQ)
     36    True
     37    sage: decidable0(FormalPowerSeriesRing(QQ))
     38    False
     39    """
     40    if K == Integer or K == int:
     41        return True
     42    if K == Rational:
     43        return True
     44    if isinstance(K,RationalField):
     45        return True
     46    if isinstance(K,RealField):
     47        return True
     48    if isinstance(K,PolynomialRing_field):
     49        return True
     50    #            if isinstance(K,RealField):
     51    #                return true
     52    #            if isinstance(K,ComplexField_class):
     53    #                return true
     54    return False
     55
     56def _isnat(n):
     57    """
     58    Returns True if n is a non-negative int or Integer.
     59
     60    sage: from sage.rings.formal_powerseries import _isnat
     61    sage: _isnat(5r)
     62    True
     63    sage: _isnat(5)
     64    True
     65    sage: _isnat(0)
     66    True
     67    sage: _isnat(-1r)
     68    False
     69    """
     70    if isinstance(n,int) or isinstance(n,Integer):
     71        if n >= 0:
     72            return True
     73    return False
     74
     75def _assert_nat(n):
     76    """
     77    Throws an exception if not _isnat(n).
     78    sage: from sage.rings.formal_powerseries import _assert_nat
     79    sage: #try:
     80    sage: #    _assert_nat(-1)                                             
     81    sage: #except AssertionError:
     82    sage: #    print 'bang'
     83    """
     84    assert _isnat(n), repr(n)+ " must be natural number."
     85
     86class FormalPowerSeriesRing(Ring):
     87    def by_lambda(self,f,min_index=0):
     88        """
     89        Returns the powerseries with coefficients f(n).
     90
     91        Alternative expression:
     92            self.by_lambda(f) == self(f)
     93            self.by_lambda(f,min_index) == self(f,min_index)
     94
     95        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     96        sage: P = FormalPowerSeriesRing(QQ)
     97        sage: P.by_lambda(lambda n: 1/factorial(n))
     98        [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...]
     99
     100        Initialization as formal Laurant series:
     101        sage: P.by_lambda(lambda n: n,-3)               
     102        [-3, -2, -1; 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, ...]
     103
     104        Corruptly setting min_index=3
     105        sage: P.by_lambda(lambda n: n,3)
     106        [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...]
     107
     108        Note that functions can not be serialized/pickled.
     109        If you want to have a serializable/picklable object, you can derive
     110        from FormalPowerSeries and define the method coeffs.
     111       
     112        sage: from sage.rings.formal_powerseries import FormalPowerSeries
     113        sage: #class F(FormalPowerSeries):^J    def coeffs(self,n): return n
     114        sage: #F(P)
     115        """
     116        return FormalPowerSeries(self,f,min_index)
     117
     118
     119    def by_iterator(self,g,min_index=0):
     120        """
     121        Returns a powerseries from the generator g.
     122        The powerseries coefficients start at min_index
     123        which is allowed be negative obtaining a formal Laurant series.
     124       
     125        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     126        sage: P = FormalPowerSeriesRing(QQ)
     127        sage: P.by_iterator(iter(ZZ),-2)
     128        [0, 1; -1, 2, -2, 3, -3, 4, -4, 5, -5, 6, -6, 7, -7, 8, -8, 9, -9, 10, -10, ...]
     129        """
     130
     131        return Iterated(self,g,min_index)
     132       
     133    def by_undefined(self,min_index=0):
     134        """
     135        Returns an undefined powerseries. This is useful to use method `define'.
     136
     137        Alternative expressions:
     138            self.by_undefined() == self()
     139            self.by_undefined(m) == self(None,m)
     140
     141        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     142        sage: P = FormalPowerSeriesRing(QQ)
     143        sage: a = P.by_undefined()         
     144        sage: a
     145        Undefined
     146        sage: P.by_undefined(2).min_index
     147        2
     148        """
     149        return FormalPowerSeries(self,min_index=min_index)
     150   
     151    def by_list(self,list,start=0):
     152        """
     153        Returns the powerseries with coefficients p[n] where
     154        p[n]==0 for 0<=n<start, p[m+start]==list[m] for all list indices m,
     155        and p[n]==0 for all later indices n.
     156
     157        Alternative expression:
     158            self.by_list(list) == self(list)
     159            self.by_list(list,start) == self(list,start)
     160
     161        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     162        sage: P = FormalPowerSeriesRing(QQ)                         
     163        sage: f = P.by_list([1,2,3,4,5])
     164        sage: f
     165        [1, 2, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     166        sage: P.by_list([1,2,3],5)
     167        [0, 0, 0, 0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     168        """
     169        return List(self,list,start).reclass()
     170
     171    def by_polynomial(self,p):
     172        """
     173        Returns the FormalPowerSeries from the given Polynomial.
     174        Alternative expression: self.by_polynomial(p) == self(p)
     175
     176        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     177        sage: PP = PolynomialRing(QQ,x)
     178        sage: P = FormalPowerSeriesRing(QQ)
     179        sage: pp = PP([1,2,3])
     180        sage: P.by_polynomial(pp)
     181        [1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     182        sage: x = PP(x)
     183        sage: P(2*x+x**2)
     184        [0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     185        """
     186        return self.by_list(p.padded_list())
     187
     188    def by_powerseries(self,p):
     189        """
     190        Returns the FormalPowerSeries from the given PowerSeries.
     191        Alternative expression: self.py_powerseries(ps)==self(ps)
     192
     193        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     194        sage: PS = PowerSeriesRing(QQ,x)
     195        sage: FPS = FormalPowerSeriesRing(QQ)
     196        sage: FPS.by_powerseries(PS([0,1,2,3]))
     197        [0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     198        """
     199        return self.by_polynomial(p.polynomial())
     200       
     201    def by_taylor(self,expr,v,at=0,**kwargs):
     202        """
     203        Returns the taylor series of `expr' with respect to `v' at `at'.
     204
     205        Alternative expressions:
     206            self.by_taylor(expr,v) == self(expr,v)
     207            self.by_taylor(expr,v,at) == self(expr,v,at)
     208
     209        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     210        sage: P = FormalPowerSeriesRing(QQ)
     211        sage: f = P.by_taylor(ln(x),x,1)   
     212        sage: f
     213        [0, 1, -1/2, 1/3, -1/4, 1/5, -1/6, 1/7, -1/8, 1/9, -1/10, 1/11, -1/12, ...]
     214        """
     215
     216        #print "after",min_index
     217        return Taylor(self,expr,v,at).reclass()
     218
     219    def by_constant(self,c,**kwargs):
     220        """
     221        Returns the powerseries with coefficients [c,0,0,...].
     222
     223        Alternative expression: self.by_constant(c) == self(c)
     224
     225        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     226        sage: P = FormalPowerSeriesRing(QQ)
     227        sage: f = P.by_constant(42)
     228        sage: f
     229        [42, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     230        """
     231       
     232        return Constant(self,c)
     233       
     234
     235    def __call__(self,p1=None,p2=None,p3=None,**kwargs):
     236        """
     237        Initialization by finite sequence of coefficients:
     238        Examples:
     239        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     240        sage: PQ = FormalPowerSeriesRing(QQ)
     241        sage: PQ([1,2,3])
     242        [1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     243        sage: PQ([])
     244        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     245       
     246        Initialization by coefficient function:
     247        Example:
     248        sage: PQ(lambda n: 1/factorial(n))
     249        [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...]
     250
     251        Initialization by expresion:
     252        Examples:
     253        sage: PQ(1+2*x+3*x^2,x)
     254        [1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     255        sage: PQ(exp(x),x)
     256        [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...]
     257        sage: PQ(ln(x),x,1)
     258        [0, 1, -1/2, 1/3, -1/4, 1/5, -1/6, 1/7, -1/8, 1/9, -1/10, 1/11, -1/12, ...]
     259
     260        Note: This is much slower than directly providing the coefficient function.
     261
     262        See also methods: by_const, by_undef, by_list, by_taylor, by_lambda
     263        """
     264
     265        if isinstance(p1,Integer) or isinstance(p1,int) or isinstance(p1,Rational):
     266            return self.by_constant(p1,**kwargs)
     267
     268        if isinstance(p1,list):
     269            if p2 == None:
     270                return self.by_list(p1,**kwargs)
     271            return self.by_list(p1,p2,**kwargs)
     272
     273        if isinstance(p1,Expression):
     274            if p3 == None:
     275                return self.by_taylor(p1,p2,**kwargs)
     276            return self.by_taylor(p1,p2,p3,**kwargs)
     277
     278        if isinstance(p1,Polynomial):
     279            return self.by_polynomial(p1)
     280
     281        if isinstance(p1,PowerSeries):
     282            return self.by_powerseries(p1)
     283
     284        #TODO generator if isinstance(p1,
     285
     286        if type(p1) is type(lambda n: 0):
     287            if p2 == None:
     288                return self.by_lambda(p1,**kwargs)
     289            return self.by_lambda(p1,p2,**kwargs)
     290
     291        if p1 == None:
     292            return self.by_undefined(p2)
     293
     294        raise TypeError, "unrecognized initialization input" + repr(type(p1))
     295
     296    def is_field(self):
     297        """
     298        Returns True if self is a field, i.e. if it can be used as
     299        formal laurant series.
     300
     301        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     302        sage: FormalPowerSeriesRing(IntegerRing()).is_field()
     303        False
     304        sage: FormalPowerSeriesRing(QQ).is_field()           
     305        True
     306        """
     307        return self.K.is_field()
     308
     309    def is_commutative(self):
     310        """
     311        The powerseries is commutative if the base ring is.
     312
     313        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     314        sage: FormalPowerSeriesRing(IntegerRing()).is_commutative()
     315        True
     316        """
     317        return self.K.is_commutative()
     318
     319    def is_exact(self):
     320        """
     321        The powerseries is exact if the base ring is.
     322
     323        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     324        sage: FormalPowerSeriesRing(RR).is_exact()
     325        False
     326        sage: FormalPowerSeriesRing(QQ).is_exact()           
     327        True
     328        """
     329        return self.K.is_exact()
     330
     331    def is_finite(self):
     332        """
     333        Powerseries ring is never finite (except the base ring has only 1 element).
     334        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     335        sage: FormalPowerSeriesRing(QQ).is_finite()     
     336        False
     337        """
     338        return False
     339
     340    def zero_element(self):
     341        """
     342        Returns the zero element of this power series ring.
     343
     344        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     345        sage: P = FormalPowerSeriesRing(QQ)
     346        sage: P.zero_element()             
     347        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     348        """
     349        return self.Zero
     350
     351    def one_element(self):
     352        """
     353        Returns the one element of this power series ring.
     354
     355        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     356        sage: P = FormalPowerSeriesRing(QQ)
     357        sage: P.one_element()             
     358        [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     359        """
     360        return self.One
     361       
     362    def __init__(self,base_ring):
     363        """
     364        Returns the powerseries ring over base_ring.
     365
     366        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     367        sage: FormalPowerSeriesRing(QQ)
     368        FormalPowerSeriesRing over Rational Field
     369        """
     370        if base_ring == None:
     371            return
     372
     373        self.K = base_ring
     374
     375        def PSS(seq):
     376            """ sage: None   # indirect doctest """
     377            return self.by_list(seq)
     378
     379        if self.K == int:
     380            self.K = Integer
     381
     382        K = self.K
     383        K1 = K.one_element()
     384        self.K1 = K1
     385
     386        self.Zero = Zero(self)
     387        self.One = One(self)
     388        self.Id = Id(self,min_index=1)
     389        self.Inc = Inc(self)
     390        self.Dec = Dec(self)
     391        self.Exp = Exp(self)
     392        self.Dec_exp = Dec_exp(self,min_index=1)
     393        self.Log_inc = Log_inc(self,min_index=1)
     394        self.Sin = Sin(self,min_index=1)
     395        self.Cos = Cos(self)
     396        self.Arcsin = Arcsin(self,min_index=1)
     397        self.Arctan = Arctan(self,min_index=1)
     398        self.Sinh = Sinh(self,min_index=1)
     399        self.Cosh = Cosh(self)
     400        self.Arcsinh = Arcsinh(self,min_index=1)
     401        self.Arctanh = Arctanh(self,min_index=1)
     402        self.Bernoulli = (self.Id / self.Exp.dec()).mul_fact()
     403        self.Bernoulli.__doc__ = """
     404        The n-th Bernoulli number is equal to
     405        the n-th derivative of 1/(exp(x)-1) at 0.
     406        """
     407        self.Tan = Tan(self,min_index=1)
     408        self.Tanh = Tanh(self,min_index=1)
     409        self.Xexp = Xexp(self,min_index=1)
     410        self.Lambert_w = Lambert_w(self,min_index=1)
     411        self.Sqrt_inc = Sqrt_inc(self)
     412
     413        #dont go into a recursion defining stirling1
     414        if isinstance(K,FormalPowerSeriesRing):
     415            return
     416
     417        self.Stirling1 = Stirling1(self)
     418
     419#         def lehmer_comtet(n,k): #A008296
     420#             """ sage: None   # indirect doctest """
     421#             r = 0
     422#             for l in range(k,n+1):
     423#                 r += binomial(l, k)*k**(l-k)*self.Stirling1[n][l]
     424#             return K(r)
     425
     426        self.Lehmer_comtet = Lehmer_comtet(self)
     427        self.A000248  = self.Lehmer_comtet
     428
     429        #self.selfpower_inc = PSF(lambda n: K(sum([ lehmer_comtet(n,k) for k in range(0,n+1))/factorial(n),K0))
     430        self.Selfpower_inc = Selfpower_inc(self)
     431
     432        self.Superroot_inc = Superroot_inc(self)
     433
     434        self.A003725 = A003725(self)
     435
     436        self.Selfroot_inc = Selfroot_inc(self)
     437       
     438        self.Inv_selfroot_inc = Inv_selfroot_inc(self)
     439
     440    def _repr_(self):
     441        """
     442        Description of this FormalPowerSeriesRing.
     443
     444        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     445        sage: FormalPowerSeriesRing(QQ)._repr_()           
     446        'FormalPowerSeriesRing over Rational Field'
     447        """
     448        return "FormalPowerSeriesRing over " + repr(self.K)
     449
     450    def _coerce_map_from_(self,T):
     451        """
     452        Returns true if type T can be coerced to a FormalPowerSeries
     453        with self as parent. This can be done always when
     454        T can be coerced to self.base_ring().
     455        This is used for operations like lmul and rmul.
     456
     457        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     458        sage: P = FormalPowerSeriesRing(RR)
     459        sage: P._coerce_map_from_(QQ)
     460        True
     461        sage: P._coerce_map_from_(int)
     462        True
     463        sage: P = FormalPowerSeriesRing(QQ)
     464        sage: P._coerce_map_from_(int)                     
     465        True
     466        sage: P._coerce_map_from_(RR)     
     467        False
     468        """
     469        #print self.K, T,not self.K.coerce_map_from(T) == None
     470        if not self.K.coerce_map_from(T) == None:
     471            return True
     472        if isinstance(T,FormalPowerSeriesRing):
     473            return not self.K.coerce_map_from(T.K) == None
     474        return False
     475           
     476
     477    def base_ring(self):
     478        """
     479        Returns the base ring of the powerseries ring.
     480
     481        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     482        sage: FormalPowerSeriesRing(QQ).base_ring() == QQ
     483        True
     484        """
     485        return self.K
     486
     487    def _test(self):
     488        """
     489        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     490        sage: P = FormalPowerSeriesRing(QQ)
     491       
     492        sage: P.Bernoulli
     493        [1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66, 0, -691/2730, 0, 7/6, ...]
     494
     495#         #takes too long:
     496#         sage: P(exp(x),x)-P.Exp
     497#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     498#         sage: P(log(x+1),x) - P.Log_inc
     499#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     500#         sage: P(cos(x),x)-P.Cos
     501#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     502#         sage: P(sin(x),x)-P.Sin
     503#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     504#         sage: P(arcsin(x),x) - P.Arcsin
     505#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     506#         sage: P(arctan(x),x) - P.Arctan
     507#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     508#         sage: P(sinh(x),x)-P.Sinh
     509#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     510#         sage: P(cosh(x),x)-P.Cosh
     511#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     512#         sage: P(arcsinh(x),x)-P.Arcsinh
     513#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     514#         sage: P(arctanh(x),x)-P.Arctanh
     515#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     516#         sage: P(sqrt(x+1),x)-P.Sqrt_inc
     517#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     518#         sage: P(x*exp(x),x)-P.Xexp
     519#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     520#         sage: P(exp(x)-1,x)-P.Exp.dec()
     521#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     522
     523        sage: P.Exp.dec().reclass() | P.Log_inc
     524        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     525        sage: P.Sin | P.Arcsin
     526        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     527        sage: P.Tan | P.Arctan
     528        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     529        sage: P.Tanh | P.Arctanh
     530        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     531        sage: P.Xexp | P.Lambert_w
     532        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     533        sage: P.Superroot_inc.dec().reclass() | P.Selfpower_inc
     534        [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     535        sage: P.Inv_selfroot_inc.dec().reclass() | P.Selfroot_inc
     536        [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     537
     538        sage: P.Superroot_inc ** P.Superroot_inc
     539        [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     540        sage: P.Tan - P.Sin / P.Cos
     541        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     542        sage: P.Sin*P.Sin + P.Cos*P.Cos
     543        [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     544        sage: p = P([3,2,1])
     545        sage: p.rcp()*p
     546        [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     547
     548
     549        sage: P._test()
     550        """
     551        pass
     552
     553# class UndefinedFormalPowerSeries(RingElement):
     554#     """
     555#     Undefined powerseries.
     556#     """
     557#     coeffs = None
     558#     def _repr_(a):
     559#         return "Undefined"
     560
     561class FormalPowerSeries(RingElement):
     562    """
     563    Formal power series:
     564
     565    A powerseries p is basically seen as an infinite sequence of coefficients
     566    The n-th coefficient is retrieved by p[n].
     567
     568    EXAMPLES:
     569        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     570        sage: PQ = FormalPowerSeriesRing(QQ)
     571        sage: #Predefined PowerSeries
     572        sage: expps = PQ.Exp
     573        sage: expps.polynomial(10,x)
     574        1/362880*x^9 + 1/40320*x^8 + 1/5040*x^7 + 1/720*x^6 + 1/120*x^5 + 1/24*x^4 + 1/6*x^3 + 1/2*x^2 + x + 1
     575        sage: expps
     576        [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...]
     577        sage: PQ.Zero
     578        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     579        sage: PQ.One
     580        [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     581        sage: PQ.Id
     582        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     583        sage: #finite powerseries
     584        sage: p = PQ([1,2,3,4])
     585        sage: p.polynomial(10,x)
     586        4*x^3 + 3*x^2 + 2*x + 1
     587        sage: p
     588        [1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     589        sage: one = PQ([1])
     590        sage: id = PQ([0,1])
     591
     592        sage: #power series are just functions that map the index to the coefficient         
     593        sage: expps[30]
     594        1/265252859812191058636308480000000
     595
     596        Formal Laurant Series can have negative minimum index
     597        sage: PQ(lambda n: n,-5)
     598        [-5, -4, -3, -2, -1; 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, ...]
     599
     600        Power series operations
     601        sage: p+p
     602        [2, 4, 6, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     603        sage: p-p
     604        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     605        sage: p*p
     606        [1, 4, 10, 20, 25, 24, 16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     607        sage: one/p
     608        [1, -2, 1, 0, 5, -14, 13, -4, 25, -90, 121, -72, 141, -550, 965, -844, 993, ...]
     609        sage: p.rcp()/p*p*p
     610        [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     611        sage: p**2
     612        [1, 4, 10, 20, 25, 24, 16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     613        sage: #composition only works for coefficient 0 being 0 in the second operand         
     614        sage: dexp = (expps - one).reclass()
     615        sage: expps(dexp)
     616        [1, 1, 1, 5/6, 5/8, 13/30, 203/720, 877/5040, 23/224, 1007/17280, ...]
     617
     618        sage: #we come into interesting regions ...
     619        sage: dexp(dexp)
     620        [0, 1, 1, 5/6, 5/8, 13/30, 203/720, 877/5040, 23/224, 1007/17280, ...]
     621        sage: dexp.nit(2)
     622        [0, 1, 1, 5/6, 5/8, 13/30, 203/720, 877/5040, 23/224, 1007/17280, ...]
     623        sage: dexp.it(-1)
     624        [0, 1, -1/2, 1/3, -1/4, 1/5, -1/6, 1/7, -1/8, 1/9, -1/10, 1/11, -1/12, ...]
     625        sage: dexp.it(-1)(dexp)
     626        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     627
     628        sage: hdexp = dexp.it(1/2)
     629        sage: hdexp
     630        [0, 1, 1/4, 1/48, 0, 1/3840, -7/92160, 1/645120, 53/3440640, -281/30965760, ...]
     631        sage: #verifying that shoudl be Zero                                                 
     632        sage: hdexp.it(2) - dexp
     633        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     634
     635        sage: #symbolic (parabolic) iteration                                                 
     636        sage: dexp.it(x)[5].coefficients()
     637        [[-1/180, 1], [1/24, 2], [-13/144, 3], [1/16, 4]]
     638        sage: q = dexp.it(1/x).it(x)
     639
     640        sage: expand(q[3])
     641        1/6
     642        sage: dexp[3]
     643        1/6
     644
     645        sage: #you can initialize power series with arbitrary functions on natural numbers   
     646        sage: #for example the power series of sqrt(2)^x can be given as
     647        sage: bsrt = FormalPowerSeriesRing(SR)(sqrt(2)^x,x)
     648
     649        sage: #making the 0-th coefficient 0 to get the decremented exponential
     650        sage: dbsrt = bsrt.set_item(0,0)
     651       
     652        sage: #and now starting hyperbolic iteration
     653        sage: dbsrt2 = dbsrt.it(x).it(1/x)
     654        sage: #Sage is not able to simplify
     655        sage: #simplify(dbsrt2[3])
     656        sage: #but numerically we can verify equality
     657        sage: RR(dbsrt2[3](x=0.73)-dbsrt[3]) < 10**(-18)
     658        True
     659    """
     660
     661    def __init__(self,parent,f=None,min_index=None,base_ring=None):
     662        """
     663        Returns the formal powerseries.
     664       
     665        This method should only be called from FormalPowerSeriesRing.
     666        See the description there how to obtain a FormalPowerSeries.
     667
     668        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing,FormalPowerSeries
     669        sage: FormalPowerSeries(FormalPowerSeriesRing(QQ))
     670        Undefined
     671        """
     672        if parent == None:
     673            if base_ring==None:
     674                parent = FormalPowerSeriesRing(QQ)
     675            else:
     676                parent = FormalPowerSeriesRing(base_ring)
     677
     678        RingElement.__init__(self, parent)
     679
     680        if not self.__class__.__dict__.has_key('coeffs'):
     681            self.coeffs = f
     682
     683        if min_index == None:
     684            min_index = 0
     685        self.min_index = min_index #the minimal non-zero index
     686        self._memo = {}
     687        self._powMemo = {}
     688        self._itMemo = {}
     689
     690        self.K = self.parent().K
     691        self.K1 = self.parent().K1
     692
     693        self.min_index = min_index
     694
     695        if self.min_index > 0:
     696            self._subclass(FormalPowerSeries0)
     697        #if not f==None:
     698        #    self.reclass()
     699           
     700#    def new(self,f=None,min_index=0,**kwargs):
     701#        return type(self)(self.parent(),f,**kwargs)
     702       
     703    def _subclass(self,T):
     704        """
     705        Makes the methods in T available to self.
     706
     707        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     708        sage: P = FormalPowerSeriesRing(QQ)
     709        sage: #class B:
     710        sage: #    def f(self): return 'B'
     711        sage: #b = P.Exp._subclass(B)                                         
     712        sage: #b.f()                                                         
     713        """
     714        if isinstance(self,T):
     715            return self
     716
     717#        if issubclass(T,self.__class__):
     718#            self.__class__ = T
     719#            return self
     720#
     721#        bs = ()
     722#        for C in self.__class__.__bases__:
     723#            if issubclass(T,C):
     724#                assert not T == C
     725#                bs += (T,)
     726#            else:
     727#                bs += (C,)
     728#        self.__class__.__bases__ = bs
     729
     730        class Sub(T,self.__class__): pass
     731
     732        self.__class__ = Sub
     733
     734        return self
     735       
     736    def new(self,f=None,min_index=None,**kwargs):
     737        """
     738        Returns a FormalPowerSeries from the same parent and class.
     739
     740        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     741        sage: p = FormalPowerSeriesRing(QQ)([1,2])
     742        sage: p.new(lambda n: n).parent() == p.parent()
     743        True
     744        sage: p = FormalPowerSeriesRing(QQ)([0,2])
     745        sage: type(p.new(lambda n: n)) == type(p)
     746        True
     747        sage: p = FormalPowerSeriesRing(QQ)([0,1])
     748        sage: type(p.new()) == type(p)
     749        True
     750        """
     751        res = FormalPowerSeries(self.parent(),f,min_index,**kwargs)
     752        if min_index == None:
     753            res.__class__ = self.__class__
     754            if isinstance(self,FormalPowerSeries0):
     755                res.min_index = 1
     756        return res
     757
     758#     def __nonzero__(self):
     759#         """
     760#         Returns always None,
     761#         as it is not decidable whether a powerseries is 0.
     762
     763#         sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     764#         sage: FormalPowerSeriesRing(QQ)(0).is_zero() == None
     765#         """
     766#         return None
     767
     768#     def is_one(self):
     769#         """
     770#         Returns always None,
     771#         as it is not decidable whether a powerseries is 1.
     772#         """
     773#         return None
     774
     775    def define(self,p):
     776        """
     777        Defines the power series self by another powerseries p
     778        which is allowed to refer to self.
     779
     780        For example one can define exp by taking the integral of its
     781        derivative which is again exp:
     782
     783        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     784        sage: P = FormalPowerSeriesRing(QQ)
     785        sage: f = P.by_undefined()         
     786        sage: f.define(integral(f,1))     
     787        sage: f                           
     788        [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...]
     789       
     790        Or generally one can define g = exp o f by taking the integral of
     791        g * f'. For example for f(x)=x**2, [0,0,1]:
     792        sage: g = P()
     793        sage: f = P([0,0,1])
     794        sage: fd = f.diff()
     795        sage: g.define(integral(g*fd,1))                 
     796        sage: g - (f | P.Exp)
     797        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     798        """
     799        self.coeffs = p.coeffs
     800       
     801    def __getitem__(self,n):
     802        """
     803        Returns the n-th coefficient of this powerseries: self[n].
     804        This is the coefficient of x^n.
     805
     806        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     807        sage: P = FormalPowerSeriesRing(QQ)
     808        sage: P([1,2,3])[1] == 2
     809        True
     810        """
     811
     812        if not self._memo.has_key(n):
     813            #self._memo[n] = simplify(expand(self.coeffs(n)))
     814            self._memo[n] = self.coeffs(n)
     815        return self._memo[n]
     816
     817    def __getslice__(self,i,j): # [i:j]
     818        """
     819        Returns the list of coefficients with indices in range(i,j).
     820
     821        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     822        sage: P = FormalPowerSeriesRing(QQ)
     823        sage: P(lambda n: n)[1:3]
     824        [1, 2]
     825        """
     826        return [self[k] for k in range(i,j)]
     827
     828    def set_item(a, index, value):
     829        """
     830        Returns the powerseries that has a[index] replaced by value.
     831
     832        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     833        sage: P = FormalPowerSeriesRing(QQ)
     834        sage: P([1,2,3]).set_item(1,42)
     835        [1, 42, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     836        sage: P([1,2]).set_item(0,0).min_index
     837        1
     838        """
     839        min_index = a.min_index
     840        if min_index == index and value == 0:
     841            min_index += 1
     842        return a.new(lambda n: value if n == index else a[n],min_index)
     843
     844    def __setitem__(a,index,value):
     845        """
     846        Replaces a[index] by value, returns None.
     847
     848        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     849        sage: P = FormalPowerSeriesRing(QQ)
     850        sage: p = P([1,2,3])
     851        sage: p[1] = 42
     852        sage: p
     853        [1, 42, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     854        sage: p = P([1,2])
     855        sage: p[0] = 0
     856        sage: p.min_index
     857        1
     858        """
     859        min_index = a.min_index
     860        if min_index == index and value == 0:
     861            min_index += 1
     862       
     863        a.min_index = min_index
     864        f = a.coeffs
     865        a.coeffs = lambda n: value if n == index else f(n)
     866
     867    def set_slice(a,i,j,seq):
     868        """
     869        Returns the powerseries that has a[i:j] replaced by seq, returns None.
     870
     871        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     872        sage: P = FormalPowerSeriesRing(QQ)
     873        sage: P(lambda n: n).set_slice(5,10,[42,43,44,45,46])
     874        [0, 1, 2, 3, 4, 42, 43, 44, 45, 46, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, ...]
     875        sage: P([1,2]).set_slice(0,1,[0]).min_index
     876        1
     877        """
     878
     879        min_index = a.min_index
     880        min_s=j-i
     881        for k in range(0,j-i):
     882            if not seq[k] == 0:
     883                min_s = k
     884                break
     885
     886        if i <= min_index and min_index <= i+min_s:
     887            min_index = i+min_s
     888        else:
     889            min_index = min(min_index,i+min_s)
     890        return a.new(lambda n: seq[n-i] if i<=n and n<j else a[n],min_index)
     891
     892    def __setslice__(a, i, j, seq):
     893        """
     894        Replaces a[i:j] by seq.
     895
     896        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     897        sage: P = FormalPowerSeriesRing(QQ)
     898        sage: p = P(lambda n: n)
     899        sage: p[5:10] = [42,43,44,45,46]
     900        sage: p
     901        [0, 1, 2, 3, 4, 42, 43, 44, 45, 46, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, ...]
     902        sage: p = P([1,2])
     903        sage: p[0:1] = [0]
     904        sage: p.min_index
     905        1
     906        """
     907
     908        min_index = a.min_index
     909        min_s=j-i
     910        for k in range(0,j-i):
     911            if not seq[k] == 0:
     912                min_s = k
     913                break
     914
     915        if i <= min_index and min_index <= i+min_s:
     916            min_index = i+min_s
     917        else:
     918            min_index = min(min_index,i+min_s)
     919
     920        a.min_index = min_index
     921        f = a.coeffs
     922        a.coeffs = lambda n: seq[n-i] if i<=n and n<j else f(n)
     923
     924    def extinct_before(a,min_index):
     925        """
     926        Returns the powerseries with elements before min_index replaced by 0.
     927
     928        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     929        sage: P = FormalPowerSeriesRing(QQ)
     930        sage: P(lambda n: n).extinct_before(5)
     931        [0, 0, 0, 0, 0, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...]
     932        """
     933        return ExtinctBefore(a,min_index)
     934
     935    def reclass(p):
     936        """
     937        Recalculates the class of this object which possibly changes
     938        to a subclass having more (and possibly different) operations
     939        available. Returns p.
     940
     941        Reclass queries p[0] and p[1], so for the sake of a lazy `define'
     942        this is not automatically done on creation.
     943
     944        Due to implementation effort a reclassed object can not be
     945        pickled/serialized with dumps.
     946
     947        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing, FormalPowerSeries, FormalPowerSeries0, FormalPowerSeries01
     948        sage: P = FormalPowerSeriesRing(QQ)
     949        sage: isinstance(P([0,1]).reclass(),FormalPowerSeries01)
     950        True
     951        sage: isinstance(P([0,2]).reclass(),FormalPowerSeries0)
     952        True
     953        sage: isinstance(P([1,1]).reclass(),FormalPowerSeries)
     954        True
     955        """
     956
     957#         if not decidable0(p.K):
     958#             if p.min_index > 0 and not isinstance(p,FormalPowerSeries0):
     959#                 p._subclass(FormalPowerSeries0)
     960#             return p
     961
     962        min_index = max(2,p.min_index)
     963        for n in range(p.min_index,2):
     964            if not p[n] == 0:
     965                min_index = n
     966                break
     967
     968        p.min_index = min_index
     969
     970        if p.min_index < 0:
     971            return p
     972
     973        if min_index > 0:
     974            if p[1] == 1:
     975                return p._subclass(FormalPowerSeries01)
     976            return p._subclass(FormalPowerSeries0)
     977        return p
     978           
     979    def mul_fact(a):
     980        """
     981        The sequence a[n]*n!
     982
     983        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     984        sage: P = FormalPowerSeriesRing(QQ)
     985        sage: P.Exp.mul_fact()                                   
     986        [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...]
     987        """
     988        return MulFact(a)
     989
     990    def div_fact(a):
     991        """
     992        Returns the sequence a[n]/n!.
     993
     994        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     995        sage: P = FormalPowerSeriesRing(QQ)
     996        sage: P(lambda n: 1).div_fact() - P.Exp
     997        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     998        """
     999        return DivFact(a)
     1000
     1001    def inc(a):
     1002        """
     1003        Increment: a + One.
     1004
     1005        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1006        sage: P = FormalPowerSeriesRing(QQ)
     1007        sage: P.Zero.inc()
     1008        [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1009        sage: P.Zero + P.One
     1010        [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1011        """
     1012        return IncMethod(a)
     1013
     1014    def dec(a):
     1015        """
     1016        Decrement: a - One.
     1017
     1018        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1019        sage: P = FormalPowerSeriesRing(QQ)
     1020
     1021        sage: P.Zero.dec()
     1022        [-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1023        """
     1024        return DecMethod(a)
     1025
     1026    def rmul(a,s):
     1027        """
     1028        Scalar multiplication from right with scalar s
     1029
     1030        Alternative expression: a.rmul(s) == a * s
     1031
     1032        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1033        sage: P = FormalPowerSeriesRing(QQ)
     1034        sage: p = P([1,2,3])
     1035        sage: p.rmul(2)
     1036        [2, 4, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1037        sage: p * (2/3)
     1038        [2/3, 4/3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1039        """
     1040        return RMul(a,s)
     1041
     1042    _rmul_ = rmul
     1043
     1044    def lmul(a,s):
     1045        """
     1046        Scalar multiplication from left with scalar s
     1047
     1048        Alternative expression: a.lmul(s) == s * a
     1049
     1050        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1051        sage: P = FormalPowerSeriesRing(QQ)
     1052        sage: p = P([1,2,3])
     1053        sage: p.lmul(2)
     1054        [2, 4, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1055        sage: (2/3) * p
     1056        [2/3, 4/3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1057        """
     1058        return LMul(a,s)
     1059
     1060    _lmul_ = lmul
     1061
     1062    def add(a,b):
     1063        """
     1064        Addition of two powerseries.
     1065
     1066        Alternative expression: a.add(b) == a+b
     1067
     1068        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1069        sage: P = FormalPowerSeriesRing(QQ)
     1070        sage: P([1,2,3]) + P([4,5,6])
     1071        [5, 7, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1072        """
     1073
     1074        return Add(a,b)
     1075
     1076    _add_ = add
     1077
     1078    def sub(a,b):
     1079        """
     1080        Subtraction of two powerseries: a-b.
     1081
     1082        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1083        sage: P = FormalPowerSeriesRing(QQ)
     1084        sage: P(lambda n: n) - P(lambda n: n)
     1085        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1086        sage: P([0,1]).sub(P([1,0]))
     1087        [-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1088        """
     1089        return Sub(a,b)
     1090
     1091    _sub_ = sub
     1092
     1093    def neg(a):
     1094        """
     1095        Negation of the powerseries a.
     1096
     1097        Alternative expression a.neg() == -a.
     1098
     1099        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1100        sage: P = FormalPowerSeriesRing(QQ)
     1101        sage: -P(lambda n: 1)
     1102        [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, ...]
     1103        """
     1104        return Neg(a)
     1105
     1106    _neg_ = neg
     1107
     1108    def mul(a,b):
     1109        """
     1110        Multiplication of two powerseries.
     1111        This a lazy algorithm: for initial 0 in a[k] and initial 0 in b[n-k]
     1112        the corresponding b[n-k] or a[k] is not evaluated.
     1113
     1114        Alternative expression: a.mul(b) == a*b
     1115
     1116        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1117        sage: P = FormalPowerSeriesRing(QQ)
     1118        sage: P(lambda n: 1) * P(lambda n:1)
     1119        [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, ...]
     1120        """
     1121        #multiplication of two powerseries
     1122
     1123        return Mul(a,b)
     1124
     1125    _mul_ = mul
     1126
     1127    def div(c,b):
     1128        """
     1129        Division.
     1130        Alternative expression: c.div(b) == c/b
     1131
     1132        Precondition: b != Zero
     1133        It satisfies: a/b*b=a, a*b/b=a
     1134
     1135        If c[0]==0 it returns a formal Laurant series, i.e. min_index < 0.
     1136
     1137        This operation is not lazy in b, it retrieves values starting from
     1138        i=min_index until it finds b[i]!=0
     1139
     1140        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1141        sage: P = FormalPowerSeriesRing(QQ)
     1142        sage: P([1,2,3])/P([2,3,4])
     1143        [1/2, 1/4, 1/8, -11/16, 25/32, 13/64, -239/128, 613/256, 73/512, ...]
     1144        sage: P.One/P(lambda n: (-1)**n)
     1145        [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1146        """
     1147        b.min_index = b.val()
     1148
     1149        return Div(c,b)
     1150
     1151    _div_ = div
     1152
     1153    def rcp(a):
     1154        """
     1155        Returns the reciprocal power series, i.e. One/a.
     1156        a must not be Zero.
     1157        If a[0] == 0 it returns a fromal Laurant series (min_index < 0).
     1158       
     1159        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1160        sage: P = FormalPowerSeriesRing(QQ)
     1161        sage: P(lambda n:1).rcp()
     1162        [1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1163        """
     1164        return a.parent().One/a
     1165
     1166    def npow_mult(a,n):
     1167        """
     1168        Power with natural exponent n computed in the most simple way by
     1169        multiplying the powerseries n times.
     1170        This function is cached, it remembers the power for each n.
     1171
     1172        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1173        sage: P = FormalPowerSeriesRing(QQ)
     1174        sage: P([1,2,3]).npow(2)/P([1,2,3])
     1175        [1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1176        """
     1177
     1178        _assert_nat(n)
     1179
     1180        if n==0:
     1181            return a.parent().One
     1182        if n==1:
     1183            return a
     1184        if not a._powMemo.has_key(n):
     1185            n1 = int(n) / 2
     1186            n2 = int(n) / 2 + n % 2
     1187            a._powMemo[n] = a.npow_mult(n1) * a.npow_mult(n2)
     1188        return a._powMemo[n]
     1189
     1190    def _s(f,k,m,n):
     1191        """
     1192        Returns the sum over
     1193        m_1+2*m_2+...+k*m_k = n, m_0+m_1+...+m_k = m
     1194        of
     1195        (m over m_1,m_2,...,m_k)* f[0]^{m_0} ... f[k]^{m_k}
     1196
     1197        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1198        sage: P = FormalPowerSeriesRing(QQ)
     1199        sage: P.Exp._s(3,5,2)
     1200        25/2
     1201        """
     1202#        print k,m,n
     1203        if f._powMemo.has_key((k,m,n)):
     1204            return f._powMemo[(k,m,n)]
     1205
     1206        if m == 0:
     1207            if n == 0:
     1208                return 1
     1209            return 0
     1210        if k < f.min_index:
     1211            return 0
     1212#        if n < f.min_index*m: #case only for speed optimizing
     1213#            return 0
     1214#         if k == f.min_index: #case only for speed optimizing
     1215#             if n == k*m:
     1216#                 return f[k]**m
     1217#             return 0
     1218#        print 'p'
     1219        res = 0
     1220        mk = 0
     1221        while min(f.min_index*m,0)+k*mk <= n and mk <= m:
     1222            v = f._s(k-1,m-mk,n-k*mk)
     1223            if not v == 0: #be lazy in evaluating f[k]
     1224                if not mk == 0: #be lazy in evaluating f[k]
     1225                    v *= f[k]**mk
     1226                res +=  v * binomial(m,mk)
     1227            mk+=1
     1228
     1229        f._powMemo[(k,m,n)] = res
     1230        return res
     1231       
     1232    def npow_combin(f,m):
     1233        """
     1234        Power with natural exponent m.
     1235        Computed via the cominatorial way similar to Faa di Bruno's formula.
     1236
     1237        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1238        sage: P = FormalPowerSeriesRing(QQ)
     1239        sage: P([1,2,3]).npow_combin(2)/P([1,2,3])
     1240        [1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1241        """
     1242
     1243        _assert_nat(m)
     1244        return Npow(f,m)
     1245
     1246    npow = npow_mult
     1247
     1248    def nipow(a,t):
     1249        """
     1250        Non-integer power of a (works also for integers but less efficiently).
     1251        a[0] must be nonzero and a[0]**t must be defined
     1252        as well multiplication of t with all coefficients.
     1253
     1254        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1255        sage: P = FormalPowerSeriesRing(QQ)
     1256        sage: P.Exp.nipow(1/2)
     1257        [1, 1/2, 1/8, 1/48, 1/384, 1/3840, 1/46080, 1/645120, 1/10321920, ...]
     1258
     1259        sage: PR = FormalPowerSeriesRing(RR)
     1260        sage: PR.Exp.nipow(0.5).n(20)                       
     1261        [1.0000, 0.50000, 0.12500, 0.020833, 0.0026042, 0.00026042, 0.000021701, ...]
     1262        """
     1263        return Nipow(a,t)
     1264
     1265    def pow(a,t):
     1266        """
     1267        The t-th (possibly non-integer) power of a.
     1268        a[0]**t has to be defined.
     1269
     1270        Alternative expression: a.pow(t) == a ** t
     1271
     1272        It satisfies:
     1273        a.nipow(0) = One
     1274        a.nipow(1) = a
     1275        a.nipow(s+t) == a.nipow(s)*a.nipow(t)
     1276
     1277        See also: npow, nipow
     1278       
     1279        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1280        sage: P = FormalPowerSeriesRing(QQ)
     1281        sage: P(lambda n: 1).pow(-1)
     1282        [1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1283        sage: P = FormalPowerSeriesRing(RealField(16))
     1284        sage: P([1,2,3])**(-0.37)
     1285        [1.000, -0.7400, -0.09619, 1.440, -2.228, 0.6642, 4.092, -9.079, 6.390, ...]
     1286        """
     1287
     1288        if isinstance(t,FormalPowerSeries):
     1289            P = a.parent()
     1290            return ( a.log() * t ) | P.Exp
     1291        if isinstance(t,Integer) or isinstance(t,int):
     1292            if t >= 0:
     1293                return a.npow(t)
     1294            return a.rcp().npow(-t)
     1295        return a.nipow(t)
     1296
     1297    __pow__ = pow
     1298
     1299    def sqrt(a):
     1300        """
     1301        Square root of a. a.sqrt() * a.sqrt() == a
     1302
     1303        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1304        sage: P = FormalPowerSeriesRing(QQ)
     1305        sage: P([1,2,1]).sqrt()
     1306        [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1307        sage: (P.One - P.Sin ** 2).sqrt() - P.Cos
     1308        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1309        """
     1310        return a.rt(2)
     1311
     1312    def rt(a,n):
     1313        """
     1314        n-th root of a. a.rt(n)**n == a
     1315
     1316        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1317        sage: P = FormalPowerSeriesRing(QQ)
     1318       
     1319        sage: P([1,3,3,1]).rt(3)
     1320        [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1321        """
     1322        return a.nipow(1/Integer(n))
     1323
     1324    def compose(b,a):
     1325        """
     1326        Composition (left after right), in mathematical notation b o a.
     1327        Alternative expressions: b.compose(a) == b(a) == a | b
     1328       
     1329        Precondition:
     1330            a[0] == 0.
     1331        Satisfies:
     1332            a(b).polynomial(m*n,x) == a.polynomial(m,b.polynomial(n,x))
     1333
     1334        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1335        sage: P = FormalPowerSeriesRing(QQ)
     1336        sage: P([1,2,3]).compose(P([0,1,2]))
     1337        [1, 2, 7, 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1338        sage: P([1,2,3])(P([0,1,2]))
     1339        [1, 2, 7, 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1340        """
     1341        a._assertp0()
     1342
     1343        return Compose(b,a)
     1344
     1345    __call__ = compose
     1346
     1347    def log(a):
     1348        """
     1349        Logarithm of powerseries a with a[0]==1.
     1350
     1351        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1352        sage: P = FormalPowerSeriesRing(QQ)
     1353        sage: P.Exp.log()           
     1354        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1355        """
     1356
     1357        P = a.parent()
     1358
     1359        dec_a = a.set_item(0,0)
     1360
     1361#        if decidable0(a.K):
     1362#            assert a[0] == 1
     1363        return dec_a | P.Log_inc
     1364   
     1365#     def __xor__(a,t): # ^
     1366#         #Not recognized as it seems to be mapped to ** in sage
     1367#         return NotImplemented
     1368
     1369    def bell_polynomials(a,k):
     1370        """
     1371        The sequence of Bell polynomials (partial/of the second kind)
     1372        [B_{0,k}(a[1],a[2],...),B_{1,k}(a[1],a[2],...),...]
     1373
     1374        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1375        sage: P = PolynomialRing(QQ,['c1','c2','c3','c4','c5'])
     1376        sage: c1 = P('c1'); c2= P('c2'); c3 = P('c3'); c4 = P('c4'); c5 = P('c5')
     1377        sage: c = FormalPowerSeriesRing(QQ)([0,c1,c2,c3,c4,c5])
     1378        sage: c.bell_polynomials(2)[6] == 6*c5*c1 + 15*c4*c2 + 10*c3**2
     1379        True
     1380        """
     1381        return (a.div_fact()**k).mul_fact().rmul(Integer(1)/factorial(k))
     1382
     1383    def bell_polynomial(a,n,k):
     1384        """
     1385        The Bell polynomial (partial/of the second kind).
     1386
     1387        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1388        sage: P = PolynomialRing(QQ,['c1','c2','c3','c4','c5'])
     1389        sage: c1 = P('c1'); c2= P('c2'); c3 = P('c3'); c4 = P('c4'); c5 = P('c5')
     1390        sage: c = FormalPowerSeriesRing(QQ)([0,c1,c2,c3,c4,c5])
     1391        sage: c.bell_polynomial(6,3) == 15*c4*c1**2 + 60*c3*c2*c1 + 15*c2**3
     1392        True
     1393        sage: c.bell_polynomial(6,2) == 6*c5*c1 + 15*c4*c2 + 10*c3**2
     1394        True
     1395        """
     1396        return a.bell_polynomials(k)[n]
     1397
     1398    def bell_complete(a,n):
     1399        """
     1400        The complete Bell polynomial.
     1401
     1402        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1403        sage: P = PolynomialRing(QQ,['c1','c2','c3','c4','c5'])
     1404        sage: c1 = P('c1');c2 = P('c2');c3 = P('c3');c4 = P('c4');c5 = P('c5')
     1405        sage: PS = FormalPowerSeriesRing(P)
     1406        sage: c = PS([0,c1,c2,c3,c4,c5])
     1407        sage: (PS.Exp(c.div_fact()) - c.bell_complete(5).div_fact())[1:6]
     1408        [0, 0, 0, 0, 0]
     1409        """
     1410        if n <= 0:
     1411            return a.parent().Zero
     1412       
     1413        res = a.bell_polynomials(1)
     1414        for k in range(2,n+1):
     1415            res += a.bell_polynomials(k)
     1416        return res
     1417       
     1418
     1419    def genfunc(a,n):
     1420        """
     1421        Returns the generating function of this powerseries up to term
     1422        a[n-1]*x**(n-1).
     1423
     1424        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1425        sage: P = FormalPowerSeriesRing(QQ)
     1426        sage: f = P.by_list([1,2,3],-2).polynomial(5)                 
     1427        sage: g = P.by_list([1,2,3],-2).genfunc(5)
     1428        sage: f(3.7)==g(3.7)
     1429        True
     1430        """
     1431        m = a.min_index
     1432        return lambda x: sum([a[k]*x**k for k in range(m,n)],a.K(0))
     1433       
     1434    def polynomial(a,n,x=var('x')):
     1435        """
     1436        Returns the associated polynomial for the first n coefficients.
     1437        f_0 + f_1*x + f_2*x^2 + ... + f_{n-1}*x^{n-1}
     1438
     1439        In case of a Laurant series with e.g. min_index=-2:
     1440        f_{-2} x^(-2) + f_{-1} x^(-1) + ... + f_{n-1}*x^{n-1}
     1441       
     1442        You can adjust the variable by setting x.
     1443
     1444        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1445        sage: P = FormalPowerSeriesRing(QQ)
     1446        sage: P.by_list([0,1,2]).polynomial(5).padded_list()
     1447        [0, 1, 2]
     1448        """
     1449
     1450#        return PolynomialRing(a.K,x)(sum([a[k]*x**k for k in range(n)],a.K(0)))
     1451        P = PolynomialRing(a.K,x)
     1452        m = a.min_index
     1453        if m >= 0:
     1454            return P(a[:n])
     1455        return P(a[m:n])/P(x**(-m))
     1456
     1457#     def subs(a,*args,**kwargs):
     1458#         def f(n):
     1459#             if n < a.min_index:
     1460#                 return 0
     1461#             if a[n] == 0:
     1462#                 return 0
     1463#             if a[n] == 1:
     1464#                 return 1
     1465#             return a[n].subs(*args,**kwargs)
     1466#         return FormalPowerSeries(f,a.min_index)
     1467
     1468    def _assertp0(a):
     1469        """
     1470        Asserts a[0]==0.
     1471
     1472        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1473        sage: P = FormalPowerSeriesRing(QQ)
     1474        sage: P([0,2])._assertp0()
     1475        """
     1476        assert a.min_index > 0, "min_index must be > 0, but is " + repr(a.min_index) + ". Use reclass() if necessary."
     1477
     1478    def _repr_(a):
     1479        """
     1480        Returns a string representation of a, as list of the first
     1481        coefficients.
     1482        If it is a formal Laurant series
     1483        the coefficients before index 0 are seperated by ";"
     1484
     1485        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1486        sage: FormalPowerSeriesRing(QQ).by_list([1,2,3],-2)
     1487        [1, 2; 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1488        sage: FormalPowerSeriesRing(QQ).by_list([1,2,3],2)
     1489        [0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1490        sage: FormalPowerSeriesRing(QQ).by_list([1,2,3],-2)._repr_()
     1491        '[1, 2; 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]'
     1492        """
     1493#         res = ""
     1494#         for n in range(80):
     1495#             coeff = a(n)
     1496#             s = ""
     1497#             if coeff != 0:
     1498#                 if coeff != 1:
     1499#                     s += repr(a(n)) + "*"
     1500#                 if n != 0:
     1501#                     s += "x"
     1502#                     if n != 1:
     1503#                         s += "^" + repr(n)
     1504#                 else:
     1505#                     s += repr(a(n))
     1506#                 s += " + "
     1507#             if len(res)+len(s) > 75: break
     1508#             else: res += s
     1509
     1510#         res += "O(x^" + repr(n) + ")"
     1511
     1512        if a.coeffs == None:
     1513            return "Undefined"
     1514
     1515        res = "["
     1516        if a.min_index < 0:
     1517            for k in range(a.min_index,0):
     1518                res += repr(a[k])
     1519                if k==-1:
     1520                    res += "; "
     1521                else:
     1522                    res += ", "
     1523        for n in range(80):
     1524            coeff = a[n]
     1525            s = repr(a[n]) + ", "
     1526            if len(res)+len(s) > 76: break
     1527            else: res += s
     1528
     1529        res += "...]";
     1530        #res = repr([ a(m) for m in range(10)])
     1531        return res
     1532
     1533#     def complex_contour(a,N,fp=0):
     1534#         """
     1535#         Returns a contour plot of this powerseries.
     1536#         Experimental yet.
     1537#         """
     1538#         r = abs(a[N])**(-1/Integer(N))
     1539#         l = r/sqrt(2.0)
     1540#         f = a.polynomial(N)
     1541#         x0=real(fp)
     1542#         y0=imag(fp)
     1543#         return contour_plot(lambda x,y: real(f(CC(x+i*y-fp))),(x0-l,x0+l),(y0-l,y0+l),fill=false) + contour_plot(lambda x,y: imag(f(CC(x+i*y-fp))),(x0-l,x0+l),(y0-l,y0+l),fill=false)       
     1544
     1545    def n(a,*args,**kwargs):
     1546        """
     1547        Applies n to the coefficients
     1548
     1549        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1550        sage: P = FormalPowerSeriesRing(QQ)
     1551        sage: P.Exp.n(digits=3)
     1552        [1.00, 1.00, 0.500, 0.167, 0.0417, 0.00833, 0.00139, 0.000198, 0.0000248, ...]
     1553        """
     1554        return N(a,*args,**kwargs)
     1555                   
     1556    def val(a):
     1557        """
     1558        Returns the first index i such that a[i] != 0
     1559        Does not terminate if a == 0
     1560
     1561        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1562        sage: P = FormalPowerSeriesRing(QQ)
     1563        sage: P([0,0,42]).val()
     1564        2
     1565        sage: P.by_list([1],-42).val()
     1566        -42
     1567        """
     1568        n = a.min_index
     1569        while a[n] == 0:
     1570            n += 1
     1571        return n
     1572
     1573    def __lshift__(a,m=1):
     1574        """
     1575        Returns the powerseries with coefficients shiftet to the left by m.
     1576        The coefficients a[0],...,a[m-1] get lost.
     1577
     1578        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1579        sage: P = FormalPowerSeriesRing(QQ)
     1580        sage: (P.Exp.mul_fact() << 1).div_fact() - P.Exp 
     1581        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1582        """
     1583        return Lshift(a,m)
     1584
     1585    def __rshift__(a,m=1):
     1586        """
     1587        Returns the powerseries with coefficients shifted to the right by m.
     1588        The resulting coefficients b[0],...,b[m-1] are zero.
     1589
     1590        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1591        sage: P = FormalPowerSeriesRing(QQ)
     1592        sage: (P.Exp.mul_fact() >> 1).div_fact() - P.Exp
     1593        [-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1594        """
     1595        return Rshift(a,m)
     1596
     1597    def diff(a,m=1):
     1598        """
     1599        Differentiates the powerseries m times.
     1600
     1601        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1602        sage: P = FormalPowerSeriesRing(QQ)
     1603        sage: P.Exp.diff(3) - P.Exp
     1604        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1605        sage: P.by_list([1],5).diff(5)[0] == factorial(5)
     1606        True
     1607        """
     1608        return Diff(a,m)
     1609
     1610#     def integral(a,c=0,m=1):
     1611#         """
     1612#         Integrates the powerseries m times with integration constant being 0
     1613#         """
     1614#         for i in range(-m,0):
     1615#             if a[i] != 0:
     1616#                 if m == 1:
     1617#                     raise ValueError, "Coefficient -1 must be 0, but it is: " + repr(a[-1])
     1618#                 else:
     1619#                     raise ValueError, "Coefficients from -"+repr(m)+" upto -1 must be 0, however a[" +repr(i)+"]=" + repr(a[i])
     1620
     1621#         def f(n):
     1622#             if 0 <= n and n < m:
     1623#                return c
     1624#             return a[n-m]/prod(Integer(k) for k in range(n-m+1,n+1))
     1625#         return a.new(f,a.min_index+m)
     1626
     1627    def integral(a,c=0):
     1628        """
     1629        Integrates the powerseries with integration constant c
     1630
     1631        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1632        sage: P = FormalPowerSeriesRing(QQ)
     1633        sage: P.Exp.integral(1) - P.Exp
     1634        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1635        sage: P.Exp.integral(0) + P.One - P.Exp 
     1636        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1637        """
     1638
     1639        return Integral(a,c)
     1640
     1641#     def indefinite_sum(f,c=0):
     1642#         def ids(n,m):
     1643#             N = m+1
     1644#             if n > N:
     1645#                 return 0
     1646#             if n < 0:
     1647#                 return 0
     1648#             if n == 0:
     1649#                 return c
     1650#             if n == N:
     1651#                 return 1/QQ(n)
     1652#             return - sum([ f[k]*binomial(k,n) for k in range(n+1,N+1)])/QQ(n)
     1653#         print ids(1,2), ids(2,2), ids(3,2)
     1654
     1655    ### finite approximate operations
     1656
     1657    def carleman_matrix(p,N,M=None):
     1658        """
     1659        The carleman_matrix has as nth row the coefficients of p^n.
     1660        It has N rows and N columns, except M specifies a different number of
     1661        rows.
     1662
     1663        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1664        sage: P = FormalPowerSeriesRing(QQ)
     1665        sage: P([1,1]).carleman_matrix(4) == matrix([[1,0,0,0],[1,1,0,0],[1,2,1,0],[1,3,3,1]])                     
     1666        True
     1667        """
     1668        if M == None:
     1669            M=N
     1670        return matrix([[p.npow(m)[n] for n in range(N) ] for m in range(M)])
     1671
     1672    def bell_matrix(a,N,M=None):
     1673        """
     1674        The Bell matrix with N (or M) rows and N columns of this power series.
     1675        The n-th column of the Bell matrix is the sequence of coefficients
     1676        of the n-th power of the power series.
     1677        The Bell matrix is the transpose of the Carleman matrix.
     1678
     1679        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1680        sage: P = FormalPowerSeriesRing(QQ)
     1681        sage: P([1,1]).bell_matrix(4) == matrix([[1,1,1,1],[0,1,2,3],[0,0,1,3],[0,0,0,1]])                     
     1682        True
     1683        """
     1684        if M == None:
     1685            M=N
     1686        return matrix([[a.npow(n)[m] for n in range(N)] for m in range(M)])
     1687
     1688class FormalPowerSeries0(FormalPowerSeries):
     1689    """
     1690    The formal powerseries f with f[0]==0.
     1691    """
     1692#     def __init__(self,f=None,min_index=1,**kwargs):
     1693#         """
     1694#         Initializes this FormalPowerSeries0.
     1695#         Should be called only from FormalPowerSeriesRing.
     1696
     1697#         sage: None
     1698#         """
     1699#         assert min_index >= 1
     1700#         super(FormalPowerSeries0,self).__init__(f,min_index,**kwargs)
     1701       
     1702    def __or__(a,b):
     1703        """
     1704        Composition (right after left): a|b.
     1705        For documentation see: `compose'.
     1706
     1707        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1708        sage: P = FormalPowerSeriesRing(QQ)
     1709        sage: P([0,1,2]) | P([1,2,3])
     1710        [1, 2, 7, 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1711        """
     1712        return b(a)
     1713
     1714    def nit(a,n):
     1715        """
     1716        n times iteration (n times composition), n is a natural number.
     1717
     1718        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1719        sage: P = FormalPowerSeriesRing(QQ)
     1720        sage: P(1/(x+1)-1,x).nit(2)
     1721        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1722        """
     1723        _assert_nat(n)
     1724
     1725        if n == 0:
     1726            return a.parent().Id
     1727        if n == 1:
     1728            return a
     1729        if not a._itMemo.has_key(n):
     1730            n1 = int(n) / 2
     1731            n2 = int(n) / 2 + n % 2
     1732            a._itMemo[n] = a.nit(n1).compose(a.nit(n2))
     1733        return a._itMemo[n]
     1734
     1735    def it(a,t):
     1736        """
     1737        The t-th (possibly non-integer) iteration.
     1738
     1739        It satisfies:
     1740        a.it(0) = Id
     1741        a.it(1) = a
     1742        a.it(s+t) == a.it(s)(a.it(t))
     1743
     1744        Alternative expression: a.it(t) == a & t
     1745
     1746        See also: nit, regit (for restrictions depending t's kind)
     1747
     1748        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1749        sage: P = FormalPowerSeriesRing(QQ)
     1750        sage: P([0,2]).it(-2)             
     1751        [0, 1/4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1752        sage: p = P([0,2]).regit_b(1/2)
     1753        sage: (p[0],p[1],p[2])
     1754        (0, sqrt(2), 0)
     1755        """
     1756        a._assertp0()
     1757
     1758        if isinstance(t,Integer) or isinstance(t,int):
     1759            if t >= 0:
     1760                return a.nit(t)
     1761            return a.inv().nit(-t)
     1762        return a.regit(t)
     1763
     1764    def regit(a,t):
     1765        """
     1766        Regular (non-integer) iteration (works also for integer t).
     1767        Preconditions:
     1768            a[1]**n != a[1] for all n.
     1769            Particularly a[1]!=0 and a[1]!=1.
     1770            a[1]**t must be defined.
     1771
     1772        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1773        sage: P = FormalPowerSeriesRing(QQ)
     1774        sage: P([0,2]).regit(1/2)[:3]
     1775        [0, sqrt(2), 0]
     1776        """
     1777
     1778        a._assertp0()
     1779
     1780        return Regit(a,t)
     1781
     1782    def regit_b(a,t):
     1783        """
     1784        Regular iteration via the schroeder function.
     1785
     1786        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1787        sage: P = FormalPowerSeriesRing(QQ)
     1788        sage: p = P([0,2])
     1789        sage: p.regit_b(1/2)[0:3]
     1790        [0, sqrt(2), 0]
     1791        """
     1792        s = a.schroeder()
     1793        return s.inv()(s.rmul(a[1]**t))
     1794
     1795
     1796
     1797    __and__ = it
     1798
     1799    def inv(a):
     1800        """
     1801        The inverse.
     1802
     1803        Precondition: a[1] != 0
     1804        Satisfies: a.inv()(a)=Id
     1805        Alternative expression: a.inv() == ~a
     1806
     1807        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1808        sage: P = FormalPowerSeriesRing(QQ)
     1809        sage: P.Xexp.inv() - P.Lambert_w   
     1810        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1811        sage: ~P.Inv_selfroot_inc.dec().reclass() - P.Selfroot_inc.dec()
     1812        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1813        """
     1814        a._assertp0()
     1815
     1816        return Inv(a)
     1817
     1818    __invert__ = inv
     1819
     1820#     def julia_b(a):
     1821#         """
     1822#         diff(it(a,t),t)(t=0) == ln(a[1])*julia(a)
     1823#         """
     1824#         a._assertp0()
     1825       
     1826#         Poly=PolynomialRing(a.K,'x')
     1827#         b = FormalPowerSeriesRing(Poly)()
     1828#         b.min_index = 1
     1829
     1830#         def f(n):
     1831#             """ sage: None   # indirect doctest """
     1832#             if decidable0(a.K):
     1833#                 assert a[1] != 0
     1834
     1835#             if n == 0:
     1836#                 return Poly([0])
     1837#             if n == 1:
     1838#                 return Poly([0,1])
     1839#             res = a[n]*(b[1]**n)-b[1]*a[n]
     1840#             res += sum([a[m]*b.npow(m)[n] - b[m]*a.npow(m)[n] for m in range(2,n)],a.K(0))
     1841#             res /= a[1]**n - a[1]
     1842#             return res
     1843#         b.coeffs = f
     1844
     1845#         def h(p):
     1846#             """ sage: None # indirect doctest """
     1847#             return sum([p.coeffs()[n]*n for n in range(p.degree()+1)],a.K(0))
     1848
     1849#         return a.new(lambda n: h(b[n]))
     1850
     1851    def julia(a):
     1852        """
     1853        Iterative logarithm or Julia function.
     1854        Has different equivalent definitions:
     1855        1. Solution j of: j o a = a' * j, j[1]=1
     1856        2. j = diff(f.it(t),t)(t=0)
     1857
     1858        Precondition: a[1]**n!=a[1] for all n>1
     1859
     1860        It has similar properties like the logarithm:
     1861        itlog(f^t) == t*itlog(f)
     1862
     1863        It can be used to define the regular Abel function abel(f) by
     1864        abel(f)' = 1/itlog(f)
     1865
     1866        Refs:
     1867        Eri Jabotinsky, Analytic iteration (1963), p 464
     1868        Jean Ecalle, Theorie des Invariants Holomorphes (1974), p 19
     1869
     1870        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1871        sage: P = FormalPowerSeriesRing(QQ)
     1872        sage: a = P([0,2,1])               
     1873        sage: j = a.julia()                   
     1874        sage: j(a) - a.diff()*j           
     1875        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1876        """
     1877       
     1878        return Julia(a)
     1879           
     1880       
     1881#     def itlog(a):
     1882#         """
     1883#         """
     1884
     1885#         #TODO this should be possible directly
     1886#         _t = var('_t')
     1887#         g = a.it(_t)
     1888#         def f(n):
     1889#             """ sage: None   # indirect doctest """
     1890#             return diff(g[n],_t)(_t=0)
     1891#         res = a.new(f)
     1892#         res.min_index = res.val()
     1893#         return res
     1894
     1895    def schroeder(a):
     1896        """
     1897        Returns the Schroeder powerseries s with s[1]=1
     1898        for a powerseries a with a[0]=0 and a[1]^n!=a[1] for all n
     1899
     1900        The Schroeder powerseries s is determined up to a multiplicative
     1901        constant by the functional equation:
     1902        s(a(z))=a[1]*s(z)
     1903
     1904        Let f(x) = 2*x + x**2, let s(x)=log(1+x), then s(f(x))=2*s(x):
     1905        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1906        sage: P = FormalPowerSeriesRing(QQ)
     1907        sage: P([0,2,1]).schroeder() - P.Log_inc
     1908        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1909        """
     1910        a._assertp0()
     1911
     1912        return Schroeder(a)
     1913
     1914    def inv_schroeder(a):
     1915        """
     1916        Returns the inverse Schroeder powerseries.
     1917
     1918        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1919        sage: P = FormalPowerSeriesRing(QQ)
     1920        sage: p = P([0,2,3])               
     1921        sage: p.inv_schroeder() | p.schroeder()
     1922        [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1923        """
     1924        return InvSchroeder(a)
     1925       
     1926    def abel(f):
     1927        """
     1928        The regular Abel function of a powerseries f (f[1]**n != f[1])
     1929        has the form a(x)=(ln(x)+ps(x))/ln(q)
     1930        where q=f[1]!=0,1 and ps is a powerseries
     1931       
     1932        This method returns ps.
     1933
     1934        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing,FormalPowerSeries0
     1935        sage: a = var('a')
     1936        sage: p = FormalPowerSeriesRing(PolynomialRing(QQ,a))(exp(a*x)-1,x)
     1937        sage: pah = p.abel()
     1938        sage: pac = p.abel2()
     1939        sage: pah
     1940        [0, 1/2*a/(-a + 1), (5/24*a^3 + 1/24*a^2)/(a^3 - a^2 - a + 1), ...]
     1941        sage: [pac[k] - pah[k]==0 for k in range(0,5)]
     1942        [True, True, True, True, True]
     1943        """
     1944        f._assertp0()
     1945
     1946        P = f.parent()
     1947        return (f.schroeder()<<1).dec().reclass() | P.Log_inc
     1948
     1949    def abel2(a):
     1950        """
     1951        A different implementation of the regular Abel function via
     1952        integration of 1/julia(a).
     1953
     1954        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing, FormalPowerSeries0
     1955        sage: a = var('a')
     1956        sage: p = FormalPowerSeriesRing(PolynomialRing(QQ,a))(exp(a*x)-1,x)
     1957        sage: p.abel2()
     1958        [0, -1/2*a/(a - 1), (5/12*a^3 + 1/12*a^2)/(2*a^3 - 2*a^2 - 2*a + 2), ...]
     1959       
     1960        """
     1961       
     1962        return a.julia().rcp().extinct_before(0).integral()
     1963
     1964
     1965class FormalPowerSeries01(FormalPowerSeries0):
     1966    """
     1967    The FormalPowerSeriess p with p[0]==0 and p[1]==1.
     1968    """
     1969
     1970    def valit(a):
     1971        """
     1972        Returns the first index i such that a[i] != Id[i]
     1973
     1974        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     1975        sage: P = FormalPowerSeriesRing(QQ)
     1976        sage: P([0,1,0,0,1]).valit()
     1977        4
     1978        """
     1979        if not a[0] == 0:
     1980            return 0
     1981        if not a[1] == 1:
     1982            return 1
     1983        n = 2
     1984        while a[n] == 0:
     1985            n+=1
     1986        return n
     1987
     1988#     def it_b(p,t):
     1989#         """
     1990#         A different direct implementation for `it'.
     1991
     1992#         sage: p = P.Dec_exp.it_b(1/2) 
     1993#         sage: (p | p) - P.Dec_exp   
     1994#         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     1995#         """
     1996#         N=p.valit()
     1997#         P = p.parent()
     1998#         q = P()
     1999#         def f(n):
     2000#             """ sage: None   # indirect doctest """
     2001#             if n < N:
     2002#                 return P.Id[n]
     2003#             if n == N:
     2004#                 return t * p[N]
     2005#             if n > N:
     2006#                 r=p[n]
     2007#                 r+=sum([p[m]*(q**m)[n] - q[m]*(p**m)[n] for m in range(N,n)])
     2008#                 return r
     2009           
     2010#         q.coeffs = f
     2011#         return q
     2012
     2013    def regit(a,t):
     2014        """
     2015        Regular iteration for powerseries with a[0]==0 and a[1]==1.
     2016        The iteration index t needs not to be an integer.
     2017        t should be a number in the base ring of this powerseries,
     2018        but at least have a defined multiplication with the elements
     2019        of the base ring.
     2020
     2021        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2022        sage: P = FormalPowerSeriesRing(QQ)
     2023        sage: p = P.Dec_exp.regit(1/2) 
     2024        sage: (p | p) - P.Dec_exp   
     2025        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     2026        """
     2027
     2028        return Regit01(a,t)
     2029
     2030    def julia(a):
     2031        """
     2032        Iterative logarithm or Julia function.
     2033        Has different equivalent definitions:
     2034        1. Solution j of: j o a = a' * j, j[1]=1
     2035        2. j = diff(f.it(t),t)(t=0)
     2036
     2037        It has similar properties like the logarithm:
     2038        itlog(f^t) == t*itlog(f)
     2039
     2040        It can be used to define the regular Abel function abel(f) by
     2041        abel(f)' = 1/itlog(f)
     2042
     2043        Refs:
     2044        Eri Jabotinsky, Analytic iteration (1963), p 464
     2045        Jean Ecalle, Theorie des Invariants Holomorphes (1974), p 19
     2046
     2047        sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2048        sage: P = FormalPowerSeriesRing(QQ)
     2049        sage: j = P.Dec_exp.julia()
     2050        sage: j(P.Dec_exp) - P.Dec_exp.diff() * j
     2051        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     2052        """
     2053        #diff(,t)(t=0) is the first coefficient of binomial(t,m)
     2054        #Stirling1(m)[k] is the kth coefficient of m!*binomial(t,m)
     2055       
     2056        res = Julia01(a)
     2057        #res.min_index = res.val()
     2058        return res
     2059       
     2060
     2061    def abel_coeffs(a):
     2062        """
     2063        The Abel function a of a power series f with f[0]=0 and f[1]=1
     2064        generally has the form
     2065        F(z) + p(z), where
     2066        F(z)=ji[-m]*z^{-m+1}/(-m+1)+...+ji[-2]*z^{-1}/(-1)+ji[-1]*log(z)
     2067        and p is a powerseries (with positive indices only)
     2068
     2069        ji[-1] is called the iterative residue (Ecalle)
     2070        ji is the reciprocal of the Julia function
     2071        (also called iterative logarithm) -which is meromorphic- of f
     2072
     2073        The method returns the sequence [ji[-1],[F[-m+1],...,F[-1],p[0],p[1],p[2],...]
     2074
     2075        These coefficients can be useful to compute the Abel function,
     2076        for which the powerseries is a non-converging asymptotic development.
     2077
     2078        The Abel function can then be gained by
     2079        lim_{n->oo} F(f&n(z))-n
     2080
     2081       
     2082        sage: from sage.rings.formal_powerseries import *
     2083        sage: P = FormalPowerSeriesRing(QQ)
     2084        sage: P([0,1,0,2]).abel_coeffs()
     2085        [3/2, [-1/4, 0; 0, 0, -5/4, 0, 21/8, 0, -35/4, 0, 2717/80, 0, -13429/100, 0, ...]]
     2086        sage: p = P([0,1,0,0,1,2,3])
     2087        sage: a = p.abel_coeffs()
     2088        sage: a
     2089        [6, [-1/3, 1, -1; 0, -10, 11/2, 17/9, -169/12, 349/30, 13/18, -544/21, 1727/24, ...]]
     2090        sage: ((p << 1).log().rmul(a[0]) + (p | a[1]) - a[1]).reclass()
     2091        [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     2092        """
     2093       
     2094        juli = a.julia().rcp()
     2095#         m = jul.val()
     2096#         juli = (jul << m).rcp()
     2097#         return [[ juli[m+i]/(i+1) for i in range(-m,-1) ],juli[m-1], (juli<<m).integral()]
     2098        resit = juli[-1]
     2099        #juli[-1]=0
     2100        return [resit,juli.set_item(-1,0).integral()]
     2101
     2102### Constructors ###
     2103
     2104class Constant(FormalPowerSeries):
     2105    def __init__(self,parent,c):
     2106        """
     2107        Description and tests at FormalPowerSeriesRing.by_constant
     2108        sage: None   # indirect doctest
     2109        """
     2110        FormalPowerSeries.__init__(self,parent)
     2111        self.c = parent.K(c)
     2112
     2113    def coeffs(self,n):
     2114        """ sage: None   # indirect doctest """
     2115        if n == 0:
     2116            return self.c
     2117        return 0
     2118
     2119
     2120class Iterated(FormalPowerSeries):
     2121    def __init__(self,parent,g,min_index):
     2122        """
     2123        Description and tests at FormalPowerSeriesRing.by_generator
     2124        sage: None   # indirect doctest
     2125        """
     2126        FormalPowerSeries.__init__(self,parent,min_index=min_index)
     2127        self.g = g
     2128       
     2129    def coeffs(res,n):
     2130        """ sage: None # indirect doctest """
     2131        g = res.g
     2132
     2133        if n<res.min_index:
     2134            return 0
     2135        if n==res.min_index:
     2136            return g.next()
     2137        x = res[n-1] #dummy to compute the prev value
     2138        return g.next()
     2139
     2140class List(FormalPowerSeries):
     2141    def __init__(self,parent,list,start):
     2142        """
     2143        Description and tests at FormalPowerSeriesRing.by_list
     2144        sage: None   # indirect doctest
     2145        """
     2146        FormalPowerSeries.__init__(self,parent)
     2147
     2148        l = len(list)
     2149        M=0
     2150        for M in range(l):
     2151            if list[M] != 0:
     2152                break
     2153        N=-1
     2154        for i in range(l):
     2155            N = l-i-1
     2156            if list[N] != 0:
     2157                break
     2158
     2159        self.min_index = start + M
     2160        self.max_index = start + N
     2161        self.start = start
     2162        self.list = list
     2163
     2164
     2165    def coeffs(self,k):
     2166        """ sage: None   # indirect doctest """
     2167        if k<self.min_index or k>self.max_index:
     2168            return 0
     2169        return self.list[k-self.start]
     2170
     2171class Taylor(FormalPowerSeries):
     2172    def __init__(self,parent,expr,v,at):
     2173        """
     2174        Description and tests at FormalPowerSeriesRing.by_taylor
     2175        sage: None   # indirect doctest
     2176        """
     2177        assert not v == None
     2178        assert isinstance(v,Expression) #this should be Variable
     2179
     2180        si = FormalPowerSeries.__init__
     2181        #coeffs always returns non-empty list, at least [0,0] is contained
     2182        min_index = expr.taylor(v,at,2).substitute({v:v+at}).coeffs(v)[0][1]
     2183        si(self,parent,min_index=min_index)
     2184
     2185        self.expr = expr
     2186        self.v = v
     2187        self.at = at
     2188
     2189    def coeffs(self,n):
     2190        """ sage: None   # indirect doctest """
     2191        #too slow
     2192        #if not at == 0 and n==0:
     2193        #    return expr({v:at})-at
     2194        #return simplify(diff(expr,v,n).substitute({v:at})/factorial(n))
     2195        expr = self.expr
     2196        v = self.v
     2197        at = self.at
     2198        return self.K(expr.taylor(v,at,n).substitute({v:v+at}).coeff(v,n))
     2199
     2200### Constants ###
     2201
     2202class Zero(FormalPowerSeries):
     2203    """
     2204    The zero element power series.
     2205
     2206    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2207    sage: P = FormalPowerSeriesRing(QQ)
     2208    sage: loads(dumps(P.Zero))
     2209    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     2210    """
     2211    def coeffs(self,n):
     2212        """ sage: None   # indirect doctest """
     2213        return 0
     2214
     2215class One(FormalPowerSeries):
     2216    """
     2217    The one element power series.
     2218
     2219    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2220    sage: P = FormalPowerSeriesRing(QQ)
     2221    sage: loads(dumps(P.One))
     2222    [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     2223    """
     2224    def coeffs(self,n):
     2225        """ sage: None   # indirect doctest """
     2226        if n == 0:
     2227            return self.K1
     2228        return 0
     2229
     2230class Id(FormalPowerSeries01):
     2231    """
     2232    The powerseries of the identity function id(x)=x.
     2233
     2234    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2235    sage: P = FormalPowerSeriesRing(QQ)
     2236    sage: loads(dumps(P.Id))
     2237    [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     2238    """
     2239    def coeffs(self,n):
     2240        """ sage: None   # indirect doctest """
     2241        if n == 1:
     2242            return self.K1
     2243        return 0
     2244
     2245class Inc(FormalPowerSeries):
     2246    """
     2247    The powerseries of the increment function inc(x)=x+1.
     2248
     2249    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2250    sage: P = FormalPowerSeriesRing(QQ)
     2251    sage: loads(dumps(P.Inc))
     2252    [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     2253    """
     2254    def coeffs(self,n):
     2255        """ sage: None   # indirect doctest """
     2256        if n == 0 or n == 1:
     2257            return self.K1
     2258        return 0
     2259
     2260class Dec(FormalPowerSeries):
     2261    """
     2262    The powerseries of the decrement function dec(x)=x-1.
     2263
     2264    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2265    sage: P = FormalPowerSeriesRing(QQ)
     2266    sage: loads(dumps(P.Dec))                                           
     2267    [-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
     2268    """
     2269    def coeffs(self,n):
     2270        """ sage: None   # indirect doctest """
     2271        if n == 0:
     2272            return -self.K1
     2273        if n == 1:
     2274            return self.K1
     2275        return 0
     2276           
     2277class Exp(FormalPowerSeries):
     2278    """
     2279    The powerseries of the exponential function.
     2280
     2281    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2282    sage: P = FormalPowerSeriesRing(QQ)
     2283    sage: loads(dumps(P.Exp))
     2284    [1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...]
     2285    """
     2286    def coeffs(self,n):
     2287        """ sage: None   # indirect doctest """
     2288        return self.K1/factorial(n)
     2289       
     2290class Dec_exp(FormalPowerSeries01):
     2291    """
     2292    The powerseries of the decremented exponential dec_exp(x)=exp(x)-1.
     2293
     2294    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2295    sage: P = FormalPowerSeriesRing(QQ)
     2296    sage: loads(dumps(P.Dec_exp))     
     2297    [0, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...]
     2298    """
     2299    def coeffs(self,n):
     2300        """ sage: None   # indirect doctest """
     2301        if n == 0: return 0
     2302        return self.K1/factorial(n)
     2303           
     2304class Log_inc(FormalPowerSeries01):
     2305    """
     2306    The powerseries of the logarithm at 1
     2307    or the powerseries of f(x)=log(x+1) at 0.
     2308
     2309    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2310    sage: P = FormalPowerSeriesRing(QQ)
     2311    sage: loads(dumps(P.Log_inc))
     2312    [0, 1, -1/2, 1/3, -1/4, 1/5, -1/6, 1/7, -1/8, 1/9, -1/10, 1/11, -1/12, ...]
     2313    """
     2314    def coeffs(self,n):
     2315        """ sage: None   # indirect doctest """
     2316        if n == 0: return 0
     2317        return self.K((-1)**(n+1)/Integer(n))
     2318
     2319class Sin(FormalPowerSeries01):
     2320    """
     2321    The powerseries of the sine.
     2322
     2323    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2324    sage: P = FormalPowerSeriesRing(QQ)
     2325    sage: loads(dumps(P.Sin))
     2326    [0, 1, 0, -1/6, 0, 1/120, 0, -1/5040, 0, 1/362880, 0, -1/39916800, 0, ...]
     2327    """
     2328    def coeffs(self,n):
     2329        """ sage: None   # indirect doctest """
     2330        if n % 2 == 0: return 0
     2331        return self.K((-1)**((n-1)/2)/factorial(n))
     2332
     2333class Cos(FormalPowerSeries):
     2334    """
     2335    The powerseries of the cosine.
     2336
     2337    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2338    sage: P = FormalPowerSeriesRing(QQ)
     2339    sage: loads(dumps(P.Cos))
     2340    [1, 0, -1/2, 0, 1/24, 0, -1/720, 0, 1/40320, 0, -1/3628800, 0, 1/479001600, ...]
     2341    """
     2342    def coeffs(self,n):
     2343        """ sage: None   # indirect doctest """
     2344        if n % 2 == 1: return 0
     2345        return self.K((-1)**(n/2)/factorial(n))
     2346
     2347class Arcsin(FormalPowerSeries01):
     2348    """
     2349    The powerseries of arcsin.
     2350
     2351    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2352    sage: P = FormalPowerSeriesRing(QQ)
     2353    sage: loads(dumps(P.Arcsin))
     2354    [0, 1, 0, 1/6, 0, 3/40, 0, 5/112, 0, 35/1152, 0, 63/2816, 0, 231/13312, 0, ...]
     2355    """
     2356    def coeffs(self,n):
     2357        """ sage: None   # indirect doctest """
     2358       
     2359        if n % 2 == 0:
     2360            return 0
     2361        evenprod = Integer(1)
     2362        oddprod = Integer(1)
     2363        for k in range(2,n):
     2364            if k % 2 == 0:
     2365                evenprod *= k
     2366            else:
     2367                oddprod *=k
     2368        return self.K(oddprod/evenprod/n)
     2369                   
     2370class Arctan(FormalPowerSeries01):
     2371    """
     2372    The powerseries of arctan.
     2373
     2374    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2375    sage: P = FormalPowerSeriesRing(QQ)
     2376    sage: loads(dumps(P.Arctan))
     2377    [0, 1, 0, -1/3, 0, 1/5, 0, -1/7, 0, 1/9, 0, -1/11, 0, 1/13, 0, -1/15, 0, ...]
     2378    """
     2379    def coeffs(self,n):
     2380        """ sage: None   # indirect doctest """
     2381        if n % 2 == 0: return 0
     2382        return self.K((-1)**(n/2)/Integer(n))
     2383
     2384class Sinh(FormalPowerSeries01):
     2385    """
     2386    The powerseries of sinh.
     2387
     2388    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2389    sage: P = FormalPowerSeriesRing(QQ)
     2390    sage: loads(dumps(P.Sinh))
     2391    [0, 1, 0, 1/6, 0, 1/120, 0, 1/5040, 0, 1/362880, 0, 1/39916800, 0, ...]
     2392    """
     2393    def coeffs(self,n):
     2394        """ sage: None   # indirect doctest """
     2395        if n % 2 == 0: return 0
     2396        return self.K(1/factorial(n))
     2397
     2398class Cosh(FormalPowerSeries):
     2399    """
     2400    The powerseries of cosh.
     2401
     2402    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2403    sage: P = FormalPowerSeriesRing(QQ)
     2404    sage: loads(dumps(P.Cosh))
     2405    [1, 0, 1/2, 0, 1/24, 0, 1/720, 0, 1/40320, 0, 1/3628800, 0, 1/479001600, 0, ...]
     2406    """
     2407    def coeffs(self,n):
     2408        """ sage: None   # indirect doctest """
     2409        if n % 2 == 1: return 0
     2410        return self.K(1/factorial(n))
     2411
     2412class Arcsinh(FormalPowerSeries01):
     2413    """
     2414    The powerseries of arcsinh.
     2415
     2416    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2417    sage: P = FormalPowerSeriesRing(QQ)
     2418    sage: loads(dumps(P.Arcsinh))
     2419    [0, 1, 0, -1/6, 0, 3/40, 0, -5/112, 0, 35/1152, 0, -63/2816, 0, 231/13312, ...]
     2420    """
     2421    def coeffs(self,n):
     2422        """ sage: None   # indirect doctest """
     2423        if n % 2 == 0:
     2424            return 0
     2425        evenprod = Integer(1)
     2426        oddprod = Integer(1)
     2427        for k in range(2,n):
     2428            if k % 2 == 0:
     2429                evenprod *= k
     2430            else:
     2431                oddprod *= k
     2432        return self.K((-1)**(n/2)*oddprod/evenprod/n)
     2433
     2434class Arctanh(FormalPowerSeries01):
     2435    """
     2436    The powerseries of arctanh.
     2437
     2438    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2439    sage: P = FormalPowerSeriesRing(QQ)
     2440    sage: loads(dumps(P.Arctanh))
     2441    [0, 1, 0, 1/3, 0, 1/5, 0, 1/7, 0, 1/9, 0, 1/11, 0, 1/13, 0, 1/15, 0, 1/17, ...]
     2442    """
     2443    def coeffs(self,n):
     2444        """ sage: None   # indirect doctest """
     2445        if n % 2 == 0: return 0
     2446        return self.K(1/Integer(n))
     2447
     2448class Tan(FormalPowerSeries01):
     2449    """
     2450    The powerseries of tan.
     2451
     2452    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2453    sage: P = FormalPowerSeriesRing(QQ)
     2454    sage: loads(dumps(P.Tan))
     2455    [0, 1, 0, 1/3, 0, 2/15, 0, 17/315, 0, 62/2835, 0, 1382/155925, 0, ...]
     2456    """
     2457    def coeffs(self,N):
     2458        """ sage: None   # indirect doctest """
     2459        if N % 2 == 0:
     2460            return 0
     2461        n = (N + 1) / 2
     2462        P = self.parent()
     2463        return P.K(P.Bernoulli[2*n] * (-4)**n * (1-4**n) / factorial(2*n))
     2464
     2465class Tanh(FormalPowerSeries01):
     2466    """
     2467    The powerseries of tanh.
     2468
     2469    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2470    sage: P = FormalPowerSeriesRing(QQ)
     2471    sage: loads(dumps(P.Tanh))
     2472    [0, 1, 0, -1/3, 0, 2/15, 0, -17/315, 0, 62/2835, 0, -1382/155925, 0, ...]
     2473    """
     2474    def coeffs(self,N):
     2475        """ sage: None   # indirect doctest """
     2476        if N % 2 == 0:
     2477            return 0
     2478        n = (N+1)/2
     2479        P = self.parent()
     2480        return P.K(P.Bernoulli[2*n] * (-1)**(2*n) * 4**n * (4**n-1) / factorial(2*n))
     2481
     2482class Xexp(FormalPowerSeries01):
     2483    """
     2484    The powerseries of x*exp(x)
     2485
     2486    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2487    sage: P = FormalPowerSeriesRing(QQ)
     2488    sage: loads(dumps(P.Xexp))   
     2489    [0, 1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, ...]
     2490    """
     2491    def coeffs(self,n):
     2492        """ sage: None   # indirect doctest """
     2493        if n==0: return 0
     2494        return self.K(1/factorial(n-1))
     2495
     2496class Lambert_w(FormalPowerSeries01):
     2497    """
     2498    The Lambert W function is the inverse of f(x)=x*e^x
     2499
     2500    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2501    sage: P = FormalPowerSeriesRing(QQ)
     2502    sage: loads(dumps(P.Lambert_w))
     2503    [0, 1, -1, 3/2, -8/3, 125/24, -54/5, 16807/720, -16384/315, 531441/4480, ...]
     2504    """
     2505    def coeffs(self,n):
     2506        """ sage: None   # indirect doctest """
     2507        if n==0: return 0
     2508        return self.K((-n)**(n-1)/factorial(n))
     2509
     2510class Sqrt_inc(FormalPowerSeries):
     2511    """
     2512    The powerseries of sqrt at 1, or of sqrt(x+1) at 0.
     2513
     2514    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2515    sage: P = FormalPowerSeriesRing(QQ)
     2516    sage: loads(dumps(P.Sqrt_inc))
     2517    [1, 1/2, -1/8, 1/16, -5/128, 7/256, -21/1024, 33/2048, -429/32768, ...]
     2518    """
     2519    def coeffs(self,n):
     2520        """ sage: None   # indirect doctest """
     2521        evenprod=Integer(1)
     2522        oddprod=Integer(1)
     2523        for k in range(2,2*n+1):
     2524            if k%2 == 0:
     2525                evenprod *= k
     2526            else:
     2527                oddprod *= k
     2528        return self.K((-1)**n *oddprod/evenprod/(1-2*n))
     2529
     2530# class BernoulliVar(FormalPowerSeries):
     2531#     """
     2532#     The variant of the Bernoulli numbers where B[1]=1/2
     2533#     instead of B[1]=-1/2
     2534#     """
     2535#     def coeffs(self,n):
     2536#         if n == 1:
     2537#             return self.K(1)/self.K(2)
     2538#         return self.Bernoulli[n]
     2539
     2540# class Faulhaber(FormalPowerSeries):
     2541#     """
     2542#     Faulhaber's formula.
     2543#     """
     2544#   def coeffs(self,n):
     2545       
     2546
     2547class Stirling1(FormalPowerSeries):
     2548    """
     2549    Returns the sequence of Stirling numbers of the first kind.
     2550    These are the coefficients of the polynomial x(x-1)(x-2)...(x-n+1).
     2551    stirling1[n][k] is the coefficient of x^k in the above polynomial.
     2552    """
     2553    def coeffs(self,n):
     2554        """ sage: None   # indirect doctest """
     2555        P = self.parent()
     2556        if n==0:
     2557            return P.One
     2558
     2559        res = Stirling1_Succ(P,self[n-1],n)
     2560        return res
     2561
     2562class Stirling1_Succ(FormalPowerSeries):
     2563    """
     2564    Constructs the sequence of Stirling numbers stirling1[n]
     2565    from g = stirling1[n-1]
     2566    """
     2567    def __init__(self,parent,g,n):
     2568        """ sage: None   # indirect doctest """
     2569        FormalPowerSeries.__init__(self,parent,min_index=1)
     2570        self.g = g
     2571        self.n = n
     2572       
     2573    def coeffs(self,k):
     2574        """ sage: None   # indirect doctest """
     2575        g = self.g
     2576        n = self.n
     2577        return g[k-1]-(n-1)*g[k]
     2578       
     2579
     2580class Lehmer_comtet(FormalPowerSeries):
     2581    """
     2582    The n-th Lehmer-Comtet number is the n-th derivative of x^x at 1.
     2583    See Sloane A008296.
     2584
     2585    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2586    sage: P = FormalPowerSeriesRing(QQ)
     2587    sage: loads(dumps(P.Lehmer_comtet))
     2588    [1, 1, 3, 10, 41, 196, 1057, 6322, 41393, 293608, 2237921, 18210094, ...]
     2589    """
     2590    def coeffs(self,n):
     2591        """ sage: None   # indirect doctest """
     2592        return self.K(sum([k**(n-k)*binomial(n,k) for k in range(n+1)]))
     2593
     2594class Selfpower_inc(FormalPowerSeries):
     2595    """
     2596    The powerseries of x^x at 1.
     2597
     2598    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2599    sage: P = FormalPowerSeriesRing(QQ)
     2600    sage: loads(dumps(P.Selfpower_inc))
     2601    [1, 1, 1, 1/2, 1/3, 1/12, 3/40, -1/120, 59/2520, -71/5040, 131/10080, ...]
     2602    """
     2603    def coeffs(self,n):
     2604        """ sage: None   # indirect doctest """
     2605        P = self.parent()
     2606        return P.K(sum([P.Stirling1[n][k]*P.A000248[k] for k in range(n+1)]))/factorial(n)
     2607
     2608class Superroot_inc(FormalPowerSeries):
     2609    """
     2610    The powerseries of the inverse of x^x developed at 1.
     2611
     2612    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2613    sage: P = FormalPowerSeriesRing(QQ)
     2614    sage: loads(dumps(P.Superroot_inc))
     2615    [1, 1, -1, 3/2, -17/6, 37/6, -1759/120, 13279/360, -97283/1008, ...]
     2616    """
     2617    def coeffs(self,n):
     2618        """ sage: None   # indirect doctest """
     2619        P = self.parent()
     2620        return P.K(sum([P.Stirling1[n][k]*Integer(1-k)**(k-1) for k in range(n+1)]))/factorial(n)
     2621
     2622class A003725(FormalPowerSeries):
     2623    """
     2624    The derivatives of exp(x*e^(-x)) at 0.
     2625
     2626    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2627    sage: P = FormalPowerSeriesRing(QQ)
     2628    sage: loads(dumps(P.A003725))
     2629    [1, 1, -1, -2, 9, -4, -95, 414, 49, -10088, 55521, -13870, -2024759, ...]
     2630    """
     2631    def coeffs(self,n):
     2632        """ sage: None   # indirect doctest """
     2633        return self.K(sum([ (-k)**(n-k)*binomial(n, k) for k in range(n+1)]))
     2634
     2635class Selfroot_inc(FormalPowerSeries):
     2636    """
     2637    The powerseries of x^(1/x) at 1.
     2638
     2639    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2640    sage: P = FormalPowerSeriesRing(QQ)
     2641    sage: loads(dumps(P.Selfroot_inc))
     2642    [1, 1, -1, 1/2, 1/6, -3/4, 131/120, -9/8, 1087/1260, -271/720, -2291/10080, ...]
     2643    """
     2644    def coeffs(self,n):
     2645        """ sage: None   # indirect doctest """
     2646        P = self.parent()
     2647        return P.K(sum([P.Stirling1[n][k]*P.A003725[k] for k in range(n+1)]))/factorial(n)
     2648
     2649class Inv_selfroot_inc(FormalPowerSeries):
     2650    """
     2651    The inverse of the self root x^(1/x) at 1.
     2652    The inverse of the self root for x=b is the fixed point of f(y)=b^y.
     2653
     2654    sage: from sage.rings.formal_powerseries import FormalPowerSeriesRing
     2655    sage: P = FormalPowerSeriesRing(QQ)
     2656    sage: loads(dumps(P.Inv_selfroot_inc))
     2657    [1, 1, 1, 3/2, 7/3, 4, 283/40, 4681/360, 123101/5040, 118001/2520, ...]
     2658    """
     2659    def coeffs(self,n):
     2660        """ sage: None   # indirect doctest """
     2661        P = self.parent()
     2662        if n<0:
     2663            return 0
     2664        if n==0:
     2665            return P.K1
     2666        r = 0
     2667        for k in range(1,n+1):
     2668            r += P.Stirling1[n][k]*P.K((k+1))**(k-1)
     2669
     2670        return P.K(r)/factorial(n)
     2671
     2672### Methods ###
     2673   
     2674class ExtinctBefore(FormalPowerSeries):
     2675    def __init__(self,a,min_index):
     2676        """
     2677        Description and tests at FormalPowerSeries.extinct_before
     2678        sage: None   # indirect doctest
     2679        """
     2680        FormalPowerSeries.extinct_before.__doc__
     2681
     2682        si = FormalPowerSeries.__init__
     2683        si(self,a.parent())
     2684        self.a=a
     2685        self.min_index=min_index
     2686
     2687    def coeffs(self,n):
     2688        """ sage: None   # indirect doctest """
     2689        if n < self.min_index:
     2690            return 0
     2691        return self.a[n]
     2692
     2693class MulFact(FormalPowerSeries):
     2694    def __init__(self,a):
     2695        """
     2696        Description and tests at FormalPowerSeries.mul_fact
     2697        sage: None   # indirect doctest
     2698        """
     2699        FormalPowerSeries.__init__(self,a.parent(),min_index=a.min_index)
     2700        self.a = a
     2701    def coeffs(self,n):
     2702        """ sage: None   # indirect doctest """
     2703        return self.a[n]*self.K(factorial(n))
     2704
     2705class DivFact(FormalPowerSeries):
     2706    def __init__(self,a):
     2707        """
     2708        Description and tests at FormalPowerSeries.div_fact
     2709        sage: None   # indirect doctest
     2710        """
     2711        FormalPowerSeries.__init__(self,a.parent(),min_index=a.min_index)
     2712        self.a = a
     2713    def coeffs(self,n):
     2714        """ sage: None   # indirect doctest """
     2715        return self.a[n]/self.K(factorial(n))
     2716       
     2717class IncMethod(FormalPowerSeries):
     2718    def __init__(self,a):
     2719        """
     2720        Description and tests at FormalPowerSeries.inc
     2721        sage: None   # indirect doctest
     2722        """
     2723        FormalPowerSeries.__init__(self,a.parent())
     2724        self.a = a
     2725    def coeffs(self,n):
     2726        """ sage: None   # indirect doctest """
     2727        if n == 0:
     2728            return self.a[0] + self.K1
     2729        return self.a[n]
     2730       
     2731class DecMethod(FormalPowerSeries):
     2732    def __init__(self,a):
     2733        """
     2734        Description and tests at FormalPowerSeries.dec
     2735        sage: None   # indirect doctest
     2736        """
     2737        FormalPowerSeries.__init__(self,a.parent())
     2738        self.a = a
     2739    def coeffs(self,n):
     2740        """ sage: None   # indirect doctest """
     2741        if n == 0:
     2742            return self.a[0] - self.K1
     2743        return self.a[n]
     2744       
     2745class RMul(FormalPowerSeries):
     2746    def __init__(self,a,s):
     2747        """
     2748        Description and tests at FormalPowerSeries.rmul
     2749        sage: None   # indirect doctest
     2750        """
     2751        FormalPowerSeries.__init__(self,a.parent(),min_index=a.min_index)
     2752        self.a = a
     2753        self.s = s
     2754
     2755    def coeffs(self,n):
     2756        """ sage: None   # indirect doctest """
     2757        return self.a[n] * self.s
     2758
     2759class LMul(FormalPowerSeries):
     2760    def __init__(self,a,s):
     2761        """
     2762        Description and tests at FormalPowerSeries.lmul
     2763        sage: None   # indirect doctest
     2764        """
     2765        FormalPowerSeries.__init__(self,a.parent(),min_index=a.min_index)
     2766        self.a = a
     2767        self.s = s
     2768
     2769    def coeffs(self,n):
     2770        """ sage: None   # indirect doctest """
     2771        return self.s * self.a[n]
     2772           
     2773class Add(FormalPowerSeries):
     2774    def __init__(self,a,b):
     2775        """
     2776        Description and tests at FormalPowerSeries.add
     2777        sage: None   # indirect doctest
     2778        """
     2779        si = FormalPowerSeries.__init__
     2780        si(self,a.parent(),min_index=min(a.min_index,b.min_index))
     2781        self.a = a; self.b = b
     2782
     2783    def coeffs(self,n):
     2784        """ sage: None   # indirect doctest """
     2785        a = self.a; b = self.b
     2786        if n < a.min_index:
     2787            if n < b.min_index:
     2788                return 0
     2789            return b[n]
     2790        if n < b.min_index:
     2791            return a[n]
     2792        return a[n]+b[n]
     2793
     2794class Sub(FormalPowerSeries):
     2795    def __init__(self,a,b):
     2796        """
     2797        Description and tests at FormalPowerSeries.sub
     2798        sage: None   # indirect doctest
     2799        """
     2800        si = FormalPowerSeries.__init__
     2801        si(self,a.parent(),min_index=min(a.min_index,b.min_index))
     2802        self.a = a; self.b = b
     2803
     2804    def coeffs(self,n):
     2805        """ sage: None   # indirect doctest """
     2806        a = self.a; b = self.b
     2807        if n < a.min_index:
     2808            if n < b.min_index:
     2809                return 0
     2810            #a[0]==0
     2811            return -b[n]
     2812        if n < b.min_index:
     2813            #b[0]==0
     2814            return a[n]
     2815        return a[n]-b[n]
     2816
     2817class Neg(FormalPowerSeries):
     2818    def __init__(self,a):
     2819        """
     2820        Description and tests at FormalPowerSeries.neg
     2821        sage: None   # indirect doctest
     2822        """
     2823        si = FormalPowerSeries.__init__
     2824        si(self,a.parent(),min_index=a.min_index)
     2825        self.a = a
     2826
     2827    def coeffs(self,n):
     2828        """ sage: None   # indirect doctest """
     2829        a = self.a
     2830        if n < a.min_index:
     2831            return 0
     2832        return -a[n]
     2833
     2834class Mul(FormalPowerSeries):
     2835    def __init__(self,a,b):
     2836        """
     2837        Description and tests at FormalPowerSeries.mul
     2838        sage: None   # indirect doctest
     2839        """
     2840        si = FormalPowerSeries.__init__
     2841        si(self,a.parent(),min_index=a.min_index+b.min_index)
     2842        self.a,self.b = a,b
     2843
     2844    def ab(self,m,n):
     2845        """
     2846        Lazy product of a[m] and b[n]
     2847        sage: None # indirect doctest
     2848        """
     2849        a = self.a
     2850        b = self.b
     2851        if a[m] == 0:
     2852            return 0
     2853        return a[m]*b[n]
     2854
     2855    def coeffs(self,n):
     2856        """ sage: None   # indirect doctest """
     2857        a,b = self.a,self.b
     2858        return sum([self.ab(k,n-k) for k in range(a.min_index,n+1-b.min_index)])
     2859
     2860class Div(FormalPowerSeries):
     2861    def __init__(self,c,b):
     2862        """
     2863        Description and tests at FormalPowerSeries.div
     2864        sage: None   # indirect doctest
     2865        """
     2866        si = FormalPowerSeries.__init__
     2867        si(self,c.parent(),min_index=c.min_index - b.min_index)
     2868        self.c=c
     2869        self.b=b
     2870
     2871    def _ab(a,m,n):
     2872        """
     2873        Lazy product of b[n] and a[m].
     2874        sage: None # indirect doctest
     2875        """
     2876        b = a.b
     2877        if b[n] == b.K(0):
     2878            return b.K(0)
     2879        return a[m]*b[n]
     2880
     2881    def coeffs(a,n):
     2882        """ sage: None   # indirect doctest """
     2883        c = a.c
     2884        b = a.b
     2885        if n<a.min_index:
     2886            return 0
     2887        r = c[n+b.min_index]
     2888        for k in range(a.min_index,n):
     2889            r -= a._ab(k,n+b.min_index-k)
     2890        return r/b[b.min_index]
     2891
     2892class Npow(FormalPowerSeries):
     2893    def __init__(self,f,m):
     2894        """
     2895        Description and tests at FormalPowerSeries.npow
     2896        sage: None   # indirect doctest
     2897        """
     2898        si = FormalPowerSeries.__init__
     2899        si(self,f.parent(),min_index=f.min_index*m)
     2900        self.f = f
     2901        self.m = m
     2902
     2903    def coeffs(self,n):
     2904        """ sage: None   # indirect doctest """
     2905        m = self.m
     2906        if n<self.min_index:
     2907            return 0
     2908        return self.f._s(n-self.f.min_index*(m-1),m,n)
     2909
     2910class Nipow(FormalPowerSeries):
     2911    def __init__(self,a,t):
     2912        """
     2913        Description and tests at FormalPowerSeries.nipow
     2914        sage: None   # indirect doctest
     2915        """
     2916        si = FormalPowerSeries.__init__
     2917        si(self,a.parent())
     2918        self.a = a
     2919        self.t = t
     2920
     2921    def coeffs(s,n):
     2922        """ sage: None   # indirect doctest """
     2923        a = s.a
     2924        t = s.t
     2925        if decidable0(a.K):
     2926            assert a[0] != 0, "0th coefficient is " + repr(a[0]) + ", but must be non-zero for non-integer powers"
     2927
     2928        da = a.set_item(0,0)
     2929       
     2930        if n>=0 and a[0] == 1:
     2931            #dont use symbolic arithmetic for ratonal powers of 1
     2932            return sum([binomial(t,k) * da.npow(k)[n] for k in range(n+1)],a.K(0))
     2933        return sum([binomial(t,k) * a[0]**t/a[0]**k * da.npow(k)[n] for k in range(n+1)],a.K(0))
     2934
     2935class Compose(FormalPowerSeries):
     2936    def __init__(self,b,a):
     2937        """
     2938        Description and tests at FormalPowerSeries.compose
     2939        sage: None   # indirect doctest
     2940        """
     2941        si = FormalPowerSeries.__init__
     2942        si(self,b.parent(),min_index=b.min_index*a.min_index)
     2943        self.b = b
     2944        self.a = a
     2945
     2946    def coeffs(self,n):
     2947        """ sage: None   # indirect doctest """
     2948        b = self.b
     2949        a = self.a
     2950        res = sum([b[k]*(a.npow(k)[n]) for k in range(n+1)])
     2951        if b.min_index < 0:
     2952            bi = a.rcp()
     2953            res += sum([b[k]*(bi.npow(-k)[n]) for k in range(b.min_index,0)],b.K(0))
     2954        return res
     2955class Lshift(FormalPowerSeries):
     2956    def __init__(self,a,m=1):
     2957        """
     2958        Description and tests at FormalPowerSeries.__lshift__
     2959        sage: None   # indirect doctest
     2960        """
     2961        si = FormalPowerSeries.__init__
     2962        si(self,a.parent())
     2963        self.a = a
     2964        self.m = m
     2965
     2966    def coeffs(self,n):
     2967        """ sage: None # indirect doctest """
     2968        return self.a[n+self.m]
     2969
     2970class Rshift(FormalPowerSeries):
     2971    def __init__(self,a,m=1):
     2972        """
     2973        Description and tests at FormalPowerSeries.__rshift__
     2974        sage: None   # indirect doctest
     2975        """
     2976        si = FormalPowerSeries.__init__
     2977        si(self,a.parent())
     2978        self.a = a
     2979        self.m = m
     2980
     2981    def coeffs(self,n):
     2982        """ sage: None # indirect doctest """
     2983        a = self.a
     2984        m = self.m
     2985        if n < m:
     2986            return 0
     2987        return a[n-m]
     2988
     2989class Diff(FormalPowerSeries):
     2990    def __init__(self,a,m=1):
     2991        """
     2992        Description and tests at FormalPowerSeries.diff
     2993        sage: None   # indirect doctest
     2994        """
     2995        si = FormalPowerSeries.__init__
     2996        si(self,a.parent(),min_index=self.deg(a.min_index,m))
     2997        self.a = a
     2998        self.m = m
     2999
     3000    def deg(self,v,m):
     3001        """ sage: None # indirect doctest """
     3002        if v >= 0:
     3003            return max(v-m,0)
     3004        return v-m
     3005
     3006    def coeffs(self,n):
     3007        """ sage: None   # indirect doctest """
     3008        a = self.a
     3009        m = self.m
     3010        if -m <= n and n < 0:
     3011            return 0
     3012        else:
     3013            return a[n+m]*prod(k for k in range(n+1,n+m+1))
     3014
     3015class Integral(FormalPowerSeries):
     3016    def __init__(self,a,c=0):
     3017        """
     3018        Description and tests at FormalPowerSeries.integral
     3019        sage: None   # indirect doctest
     3020        """
     3021        si = FormalPowerSeries.__init__
     3022        if c == 0:
     3023            si(self,a.parent(),min_index=a.min_index+1)
     3024        else:
     3025            si(self,a.parent())
     3026
     3027        self.a = a
     3028        self.c = c
     3029
     3030    def coeffs(self,n):
     3031        """ sage: None   # indirect doctest """
     3032        a = self.a
     3033        c = self.c
     3034
     3035        if n == 0:
     3036            return c
     3037        if n < a.min_index+1:
     3038            return 0
     3039
     3040        #This test can not be moved to the beginning because
     3041        #it would break lazyness
     3042        if a.min_index < 0:
     3043            assert a[-1] == 0, "Coefficient at -1 must be 0, but it is: " + repr(a[-1])
     3044        return a[n-1]/Integer(n)
     3045   
     3046class N(FormalPowerSeries):
     3047    def __init__(self,a,*args,**kwargs):
     3048        """
     3049        Description and tests at FormalPowerSeries.n
     3050        sage: None # indirect doctest
     3051        """
     3052        si = FormalPowerSeries.__init__
     3053        si(self,FormalPowerSeriesRing(n(0,*args,**kwargs).parent()))
     3054
     3055        self.a = a
     3056        self.args = args
     3057        self.kwargs = kwargs
     3058
     3059    def coeffs(self,k):
     3060        """ sage: None # indirect doctest """
     3061        return n(self.a[k],*self.args,**self.kwargs)
     3062       
     3063class Regit(FormalPowerSeries0):
     3064    def __init__(self,a,t):
     3065        """
     3066        Description and tests at FormalPowerSeries.regit
     3067        sage: None   # indirect doctest
     3068        """
     3069        si = FormalPowerSeries.__init__
     3070        si(self,a.parent(),min_index=1)
     3071        self.a = a
     3072        self.t = t
     3073
     3074    def coeffs(b,n):
     3075        """ sage: None   # indirect doctest """
     3076        a = b.a
     3077        t = b.t
     3078        if decidable0(a.K):
     3079            assert a[1]!=1, a[1]
     3080            assert a[1]!=0, a[1]
     3081
     3082        #print "(" + repr(n)
     3083        if n == 0:
     3084            #print ")"
     3085            return 0
     3086        if n == 1:
     3087            #print ")"
     3088            return a[1]**t
     3089        res = a[n]*(b[1]**n)-b[1]*a[n]
     3090
     3091        for m in range(2,n):
     3092            res += a[m]*b.npow(m)[n] - b[m]*a.npow(m)[n]
     3093
     3094        res /= a[1]**n - a[1]
     3095        #print ")"
     3096        return res
     3097
     3098class Inv(FormalPowerSeries0):
     3099    def __init__(self,a):
     3100        """
     3101        Description and tests at FormalPowerSeries.inv
     3102        sage: None   # indirect doctest
     3103        """
     3104        si = FormalPowerSeries.__init__
     3105        si(self,a.parent(),min_index=1)
     3106        self.a=a
     3107
     3108    def coeffs(b,n):
     3109        """ sage: None # indirect doctest """
     3110        a = b.a
     3111        if n==0:
     3112            return 0
     3113        if n==1:
     3114            return 1/a[1]
     3115        if n>1:
     3116            return - sum([ b[k]*a.npow(k)[n] for k in range(1,n)])/a[1]**n
     3117
     3118class Julia(FormalPowerSeries01):
     3119    def __init__(self,a):
     3120        """
     3121        Description and tests at FormalPowerSeries.julia
     3122        sage: None   # indirect doctest
     3123        """
     3124        si = FormalPowerSeries.__init__
     3125        si(self,a.parent(),min_index=1)
     3126        self.a=a
     3127
     3128    def coeffs(j,n):
     3129        """ sage: None #indirect doctest """
     3130
     3131        a = j.a
     3132        if n < 0:
     3133            return 0
     3134
     3135        if n == 1:
     3136            return 1
     3137       
     3138        ap = a.diff()
     3139
     3140        r = a.K(0)
     3141
     3142        for k in range(1,n):
     3143            r+=ap[n-k]*j[k]
     3144
     3145        for k in range(1,n):
     3146            r-=j[k]*a.npow(k)[n]
     3147
     3148        return r/(a[1]**n-a[1])
     3149
     3150class Schroeder(FormalPowerSeries01):
     3151    def __init__(self,a):
     3152        """
     3153        Description and tests at FormalPowerSeries.schroeder
     3154        sage: None   # indirect doctest
     3155        """
     3156        si = FormalPowerSeries.__init__
     3157        si(self,a.parent(),min_index=1)
     3158        self.a=a
     3159
     3160    def coeffs(s,n):
     3161        """ sage: None   # indirect doctest """
     3162        a = s.a
     3163        if decidable0(a.K):
     3164            assert a[1] != 0, a[1]
     3165            assert a[1] != 1, a[1]
     3166       
     3167        q = a[1]
     3168
     3169        if n <= 0:
     3170            return 0
     3171        if n == 1:
     3172            return 1
     3173        return sum([s[m]*a.npow(m)[n] for m in range(1,n)])/(q - q**n)
     3174
     3175class InvSchroeder(FormalPowerSeries01):
     3176    def __init__(self,a):
     3177        """
     3178        Description and tests at FormalPowerSeries.inv_schroeder
     3179        sage: None   # indirect doctest
     3180        """
     3181        si = FormalPowerSeries.__init__
     3182        si(self,a.parent(),min_index=1)
     3183        self.a=a
     3184
     3185    def coeffs(s,n):
     3186        """sage: None #indirect doctest"""
     3187        a = s.a
     3188        q = a[1]
     3189
     3190        if n <= 0:
     3191            return 0
     3192        if n == 1:
     3193            return 1
     3194        return sum([a[m]*s.npow(m)[n] for m in range(2,n+1)])/(q**n-q)
     3195   
     3196class Regit01(FormalPowerSeries01):
     3197    def __init__(self,a,t):
     3198        """
     3199        Description and tests at FormalPowerSeries01.regit
     3200        sage: None   # indirect doctest
     3201        """
     3202        si = FormalPowerSeries.__init__
     3203        si(self,a.parent(),min_index=1)
     3204        self.a = a
     3205        self.t = t
     3206
     3207    def coeffs(self,n):
     3208        """ sage: None   # indirect doctest """
     3209        a = self.a
     3210        t = self.t
     3211        if decidable0(a.K):
     3212            assert a[0] == 0, "The index of the lowest non-zero coefficient must be 1, but is " + repr(a.min_index)
     3213            assert a[1] == 1, "The first coefficient must be 1, but is " + repr(a[1])
     3214
     3215        if n == 1: return 1
     3216        #def c(m):
     3217        #    return (-1)**(n-1-m)*binomial(t,m)*binomial(t-1-m,n-1-m)
     3218        #res = sum([c(m)*a.nit(m)[n] for m in range(n)],a.K(0))
     3219        #return res
     3220
     3221        r = a.K(0)
     3222        for m in range(n):
     3223            s = a.K(0)
     3224            for k in range(m+1):
     3225                s += binomial(m,k)*(-1)**(m-k)*a.nit(k)[n]
     3226            s *= binomial(t,m)
     3227            r += s
     3228        return r
     3229
     3230class Julia01(FormalPowerSeries01):
     3231    def __init__(self,a):
     3232        """
     3233        Description and tests at FormalPowerSeries01.julia
     3234        sage: None   # indirect doctest
     3235        """
     3236        P = a.parent()
     3237        si = FormalPowerSeries.__init__
     3238        si(self,P,min_index=1)
     3239        self.a=a
     3240
     3241    def coeffs(self,n):
     3242        """ sage: None # indirect doctest """
     3243        a = self.a
     3244        r = a.K(0)
     3245        P = a.parent()
     3246
     3247        for m in range(n):
     3248            s = a.K(0)
     3249            for k in range(m+1):
     3250                s += binomial(m,k)*(-1)**(m-k)*a.nit(k)[n]
     3251            s *= P.Stirling1[m][1]/factorial(m)
     3252            r += s
     3253
     3254        return r
     3255