GSOC 2016 Series Expansion - wrat/sympy GitHub Wiki

GSOC 2016- Abhishek Verma: Series Expansion

About Me

My name is Abhishek Verma and I am a second year Computer Science student at National Institute of Technology,JAIPUR.

Email : [email protected]

Github : wrat

Time-zone: UTC + 5:30

Programming And Mathematical Background

Presently I am pursuing computer science course and I have worked on various language (C and Java and python) with project in python related to Image processing(OpenCv). I have used Numpy in my project.My area of interest belongs to Artificial Intelligence(Machine Learning,NLP,Voice Recognition system).

Math has been a strong point for me. I have taken relevant courses in Math like Linear and Abstract algebra and I am quite familiar with the concepts.Presently I am taking Discrete Mathematics Course.

Contribution to sympy

Resently I have started contributing in sympy . pull requests are related to limits and Integration.

  • Merged:

  • #10792 integrals: Improved heurisch.py to get Unevaluated Integral Object of non-polynomial.

  • Unmerged:

    • #10801 limit working now with binomial

Issue Reported

  • #10680 Integral is not working non-polynomial( xlog(xlog(x**log(x))) ).

  • #10868 sympy is unable to solve limit for function log( x ) + asech( x ) at both point x = 0 and x = oo .

Project:

The Project is divided into sub-project based on Improvement(Method Implementation) and Extension.In Gsoc 2014 and Gsoc 2015 lot of work has been already done by Avichal Dayal and Sartaj singh on FormalPowerSeries and Asymtotic Expansion.SymPy now has support for Formal Power Series (series.formal).The algorithm[1] is more or less complete.According to me sympy need a good User-Interface for it's Series Representation and expansion for user.As like given in UD - Function expansions and series - user interface .

Some method has been implemented in FormalPowerSeries But not all.There lots of thing for extension and Improvement in module.The module should be made faster.There are also a lot of XFAIL tests that can be made to pass.I have devided it in form of Sub project.

Here ,Sub-Project.

  • Improvement

    • Work on issue
      • XFAIL TEST for Formal Series
      • Improve limit Algorithm (Comment - How do you plan to do this?)
      • Implement some method in present Formal Series
        • Inversion
        • composition
        • Convolution
        • Series Expansion With Assumptions of variable

  • Extension of Formal Power series

    • Implement Formal Power series for Special Function/Orthogonal Polynomial
    • Implementing Multi-Variate Formal Power series

  • Unifying series and expansion

Timeline:

* Week 1: Made Pass All @XFAIL TEST - CASES In Formal Series
* Week 2: Improve the limit or solve issue related to limit 
* Week 3-4: Implement Extra method in Formal Series(Inverse ,composition,convolution,series with assumption,etc)
* Week 5: Add test-cases for FormalSeries for new operation
* Week 6: Implement FormalSeries Special Funtion and Add Test - cases(fibbonaci,airy function,Bessel functions 
          etc.) 
* Week 7: Implement Multivariate Series for Normal function.
* Week 8: Add Test Case For Multivariate -Series.
* Week 9-10: Unifying all type of series(Import Whole Implementation of Formal Series to new series module 
             un_series() )
* Week 11: Work on Internal Representation of Series(Differentiate all types series)
* Week 12: Add test cases for new series module (un_series)
* Week 13: Documentation, docstrings and stylistic changes where necessary

Present Formal Series

>>> ps=fps(sin(x),x)
>>> ps
     3     5        
    x     x     βŽ› 6⎞
x - ── + ─── + O⎝x ⎠
    6    120        
>>>

Gives series of function with order term.

Main thing it can be represented like infinite series as mathmatica and maple compute for its series output.

>>> ps.infinite
       ∞                                                             
    ________                                                         
    β•²                                                                
     β•²       ⎧                k   1                                  
      β•²      βŽͺ                ─ - ─                                  
       β•²     βŽͺ                2   2  k                               
        β•²    βŽͺ            -1/4     β‹…x                                
         β•²   βŽͺ────────────────────────────────────  for Mod(k, 2) = 1
          β•²  ⎨               βŽ›     k   1⎞ βŽ›k   1⎞                    
x +       β•±  βŽͺRisingFactorial⎜3/2, ─ - β”€βŽŸβ‹…βŽœβ”€ - β”€βŽŸ!                   
         β•±   βŽͺ               ⎝     2   2⎠ ⎝2   2⎠                    
        β•±    βŽͺ                                                       
       β•±     βŽͺ                 0                        otherwise    
      β•±      ⎩                                                       
     β•±                                                               
    β•±                                                                
    β€Ύβ€Ύβ€Ύβ€Ύβ€Ύβ€Ύβ€Ύβ€Ύ                                                         
     k = 2                                                           
>>>

Representation of infinite series.We can also play with infinite series. Some other operation like 1.convert to polynomial

>>> ps.polynomial()
  5    3    
 x    x     
─── - ── + x
120   6     
>>>

2.Truncate the series

>>> ps.truncate(10)
     3     5     7       9           
    x     x     x       x       βŽ› 10⎞
x - ── + ─── - ──── + ────── + O⎝x  ⎠
    6    120   5040   362880         
>>>

3.Function

>>> ps.function
sin(x)
>>>

4.We can differentiate and integrate

>>> ps.integrate()
      2    4        
     x    x     βŽ› 6⎞
-1 + ── - ── + O⎝x ⎠
     2    24        
>>>

Relevant Issue https://github.com/sympy/sympy/issues/9885

Improvement

Working On issue Related To series and limits -

Formal #XFAIL test case:

  • FormalPowerSeries:

     * Methods:
   1. simpleDE()

    ```
    Converts a function to a simple differential equation
    ```

  2. DE-to-RE()

    ```
     Converts a differential equation to recurrence equation
    ```

  3. solveRE()

   ```
    Solves recurrence equation to find a formula
   ```
  4. as_summation()

   ```
   Writes the infinite series in summation form
   ```

Most of test cases fails in these methods.I have observed that XFAIL Test cases most of time fail while Computing of DE(solve_DE() unable to solve the Differential Equation) . Here Example from XFAIL test case -

>>>f=x**n*sin(x**2)    
>>>fps(f, x).truncate(8)
>>>x**2*x**n - x**6*x**n/6 + O(x**(n + 8), x)

But it gives NotImplemented Error. Point in Regarding to solve the Issue: * Differential Equation from solve_de() method:

    ```
     DE : x**2*Derivative(g(x), x, x) + (-2*n*x - x)*Derivative(g(x), x) + (n**2 + 2*n + 4*x**4)*g(x)
     RE : (_k**2 - 2*_k*n + 6*_k + n**2 - 6*n + 8)*g(_k + 4) + 4*g(_k)
     sol = _solve_hyper_RE(f, x, RE, g, k)
     sol=None(which not be like this)
   ```

   So XFAIL Can be solve by Efficient Implementation of above method and If something to implement I am going to 
   implement it.

But If We observe the actual result then it will clear that how we have to solve this.

>>> x**n*fps(sin(x**2),x,0).truncate(18)
   βŽ›      6    10    14          ⎞
 n ⎜ 2   x    x     x       βŽ› 18⎞⎟
x β‹…βŽœx  - ── + ─── - ──── + O⎝x  ⎠⎟
   ⎝     6    120   5040         ⎠
>>>

here we can observe that expansion of fps(x**n*sin(x**2),x) is similar as if we separate the term x**n and expand sin(x**2) and multiply it x**n. So I am going to implement a module which is going to care about this.

>>>fps(f(x),x)

if f(x) has term like (x**n) then separate it from function and compute series for remaining function and then  
multiply with whole series with x**n.

Issue That I have found Regarding to Series Expansion:

1.Unable to compute series for function who has binomial Expansion ,like

>>> fps((1+x)**n,x)

Nothing is shown by present algorithm

2.Another terms like xlog(x) (in general xf(x) )

>>> fps(x**log(x),x)
 log(x)
x      (from present fps)

. But we have expansion for it

>>>fps(x**log(x),x,1)

1+ (x-1)**2 - (x-1)**3 + 17 * (x-1)**4 - 11  * (x-1)**5 + O((x-1)**6)
                         __              __
                         12              6

there are similar function whose expansion can not be computed by present series algorithm - x**sin(x) , exp(sin(x)), log(sin(x)) etc.

  • Except x**f(x) , other function can be solved by composition of series ("exp(sin(x)) composition of exp(x) and sin(x) ").

Improvement in limit algorithm:

Present Limit computation is based on gruntz algorithm and use series to get leading_term.Main step of this algorithm is to make set of Mrvterm(Most Rapid Varying term).Gruntz use series for limits calculation, and series algorithm use limits to calculate terms. gruntz working for exp-log function.and rewriting as tractable function enable most function to compute its limit.

Also gruntz can solve only fix Number of function.Most of function are not working Yet.Present algorithm can be extend further to work for more function.

Points Regarding to Improvement of limit

  • limit(3n*3(-n - 1)*(n + 1)2/n2, n, oo)

1

which is wrong ,here we required good checking in "ok() lamda" function function. Relevent issue #10610.

  • limit(4xhyper([S(1)/2],[S(3)/2, S(3)/2], -x**2), x, oo) which need expand(hyperexpand) and implementation of tractable function.

  • I can observe that , We need a Module which the keep track of function either it is compatible for gruntz algorithm or not.that module can check ,if function need simplification ,expand , or variable manipulation so that it become compatible for gruntz.ex -: hyper-function.

Here Example an

  • limit of fraction with oscillating term in the numerator:
>>> limit((n+cos(n))/n,n,oo)
>>>  n + cos(n)
lim ──────────
nβ”€β†’βˆž    n

present algorithm is unable to solve it . As far as I understand the algorithm, it is not able to handle oscillating stuff due to the properties of these function fields.we know the behaviour of cos at infinity. If we could find a way to express the fact that cos(x) is bounded oscillatory as x->oo, gruntz could probably be adapted to compute limits like this. There is some work done on this topic[5].This is based on so called transseries . I am going through book ,It will take some time to understand it.

Implement some method in present Formal Series

When we manipulate power series, it is sometimes convenient to think of the series as representing functions, which we can, for example, compose or invert, convolute.

  • Computing the Coefficient of power of series (f(x)**n where f(x) is formal power series)

Let A formal Power series given as w(x) = summation( a_n * x_n, (n, 0, oo) )

and a_0 != 0 then coefficient of Power series raised to powers can be computed by recurrence relation

w(x)**m = summation( c_n * x_n ( n*m, 0, oo))

c_0 = ( a_0 **m )
c_k = (1 /k * a_0) * summation( ( m*j - k + j) * a_j * c_(kβˆ’ j) ,( j, 1, k) )
follow this recurrence relation.


Sample run
>>>ps = fps( sin(x), x)
>>> ps.coefficient(n=3 ,pow=10) only one change in coefficient module, we can pass one extra parameter which is 
 (raised power).

>>> [1, - 5/3 , -59 / 72]
  • Convolution

the Cauchy product is the Discrete_convolution of two infinite series or can used for power series.

s1 = summation( a_n1 * x_n1, (n1, 0, oo) )
s2 = summation( a_n2 * x_n2, (n2, 0, oo) )
convolute(s1,s2) = summation( c_k * x_k, (k, 0, oo) )
c_k = summation( a_k * b_(k-l), (l, 0, k) )
  • Composition for example given two formal series:
f(x)=summation(a_n * x_n, (n, 0, oo))
g(x)=summation(b(n)*x**n, (n, 0, oo))

one can compute the composition like-

g(f(x))= summation(b_n * f(x)**n , (n, 0, oo))) =summation(c_n * X**n ,(n,0, oo))

where the coefficients c_n are determined by "expanding out" the powers of f(x): There exists a closed formula for composition of two power series and also to find the coefficients. This link mentions them clearly: http://en.wikipedia.org/wiki/Fa%C3%A0_di_Bruno%27s_formula

Example:

>>>g(x) = fps(exp(x),x)
         2    3    4     5        
        x    x    x     x     βŽ› 6⎞
1 + x + ── + ── + ── + ─── + O⎝x ⎠
        2    6    24   120        
>>> 
>>>f(x) = fps(sin(x),x)
     3     5        
    x     x     βŽ› 6⎞
x - ── + ─── + O⎝x ⎠
    6    120        

>>> fps(exp(x),x).ComposeSeries(sin(x)) "This will compose the exp(x)'s series with sin(x)'s series "
         2    4    5             
        x    x    x     βŽ› 6⎞
1 + x + ── - ── - ── + O⎝x ⎠
        2    8    15           
This replaces the variable in the power series for exp(x) by a power series for sin(x).
  • Inversion for A = a_0 + a_1 * x + .... Suppose that A has an inverse B = b_0 + b_1 * x + .... then the constant term a_0 b_0 of A * B is the constant term of the identity series, i.e. , it is 1.

b_0 = 1 / a_0 (provided a_0 is invertible)

 >>> ps1 = Summation(a(n)*x**n, (n, 0, oo))
 >>> ps1.invert()
     Summation(b(n)*x**n, (n, 0, oo))
     where b(n) = -Summation(a(i) * b(n-i), (i, 0, oo))/a(0) (Explicit Recursive Formula)

Example:


>>> fps(sin(x),x)
     3     5        
    x     x     βŽ› 6⎞
x - ── + ─── + O⎝x ⎠
    6    120        

>>> fps(sin(x),x).inverse()
     3      5        
    x    3β‹…x     βŽ› 6⎞
x + ── + ──── + O⎝x ⎠
    6     40

  • Series Expansion With Assumptions of variable In present formal series if we want compute the series of piecewise function , Lets take example:

     >>>fps(acos(x),x,1).truncate(2)

            _______   √2β‹…β…ˆβ‹…βˆšxβ‹…(x - 1)   βŽ›       2       ⎞
    √2β‹…β…ˆβ‹…β•²β•± x - 1  - ─────────────── + O⎝(x - 1) ; x β†’ 1⎠
                     12                            

  >>> p = Piecewise((acos(x), x>1))
  >>> fps(p,x,1)
    {acos(x)  for x > 1
  >>>

    But it can be solvable by mathmatica series methods: Example from mathmatica:

    >>> Series(p,{x,1,1})
        _____    βŽ›                    ⎞
 2β‹…β…ˆβ‹…β•²β•± x - 1 + O⎝(x - 1)**3/2 ; x β†’ 1⎠

Extension:

Special functions series implementation:

Also I was thinking of having formalSeries made separately for special functions similar to what eval_aseries and eval_nseries do. Although this algorithm can find series for error functions, exponential
integral functions etc, as well, it will be better to implement for them separately. here list of these special function-

The algorithm has been extended to handle many special functions(Required steps already implemented), in particular orthogonal poly-nomials.list of function for which we will compute FormalSeries-

  1. Fibonacci polynomials Fibonacci(n,x)
  2. the Bessel functions BesselJ(n,x), BesselY(n,x), BesselI(n,x), and BesselK(n,x)
  3. Hankel functions Hankel1(n,x), and Hankel2(n,x)
  4. the Kummer functions KummerM(a,b,x), and KummerU(a,b,x)

other function also there in the list. .Compute Series -

  Example From Reference[1]:

  ```
   we consider the Fibonacci polynomials
      RE : Fn (x) = x F nβˆ’1 (x) + F nβˆ’2 (x)
      F0 (x) = 0
      F1 (x) = 1.
  ```

As above given four step for implementing the Formal Series ,same Approach can be also implemented for these
type function. But here several things that We should keep in our mind that:

  1. the only precondition which must be met by a function to be covered by the algorithm is that its derivative is defined in terms of functions which also may be handled by our algorithm.

  2. The general derivative of an orthogonal polynomial of degree n can be defined in terms of the orthogonal polynomials of degree n,and n βˆ’ 1.

  3. We can teach our procedure to use this recurrence equation by assigning the table FPSRecursion.The second index specifies the number of arguments of the function family.

  4. Derivative Rule -

        ```
     Fn (x).diff(x,1) = (n βˆ’ 1) * x * Fn (x) + 2 * n * Fnβˆ’1 (x) / (x**2+4)
    ```

Step By step Solution for Fibbonacci:

  1. solve_DE(Fibbonaci(n,x) , x)

     ```
      if(f(x) is fibbonaci type):
       make_de(use fibbinaci derivative rule):
         DE : (x**2+4) * F(x).diff(x,2) + (n-1) *(n+1) *F(x) + 3*x * F(x).diff(x,1)        
     ```
  2. SolveRE():

   ```
   FPSRecursion() ### this table to solve RE
    Recurrence Rule :
    x = 0: Fn (0) = F nβˆ’2 (0)
  
    x=0 as initial point 
   ```

  3. _compute_fps(Fibonacci(n,x), x):

    ```
 
 Series: (sin(1 / 2*pi*n)**2 * (n-1)*(n-2) * Summation((-1)**k (- 1/2 n + 1/2 + k)! (1/2 n + 1/2 + k)! * x**2k,{ 
                k, 0 , oo})) / ((n -2*k-1)*(n + 2 k + 1) * (2 k)!) 
               
    + n*cos(1/2 n Pi)**2 * Summation( (-1)**k (-1/2*n + k)! * (1/2 n + k)! x**(2 k + 1) / (2 k + 1)!,k, 0 , oo}) 
                                                     /  (- 1/2 n)! * (1/2 n)! 

 ```

Implementation of Multivariate Series:

Single Variable FormalPowerSeries can be easily computed from present algorithm but ,it is not Not Implemented for Multivariate function.But it can be extended for multivariable.Here What Mathmatica Does with these type series

Series[f,{x,x0,nx},{y,y0,ny},…]
successively finds series expansions with respect to x, then y, etc.

Example from mathmatica-

Power series in two variables:
 Series[sin[x+y],{x,0,3},{y,0,3}]
 >>>( (y - y**3/6 + O[y]**4) + (1- y**2/2 +O[y]**4)*x +(-y/2 +y**3/12 +O[y]**4) * x**2 +(-1/6 + y**2 /12 +O[y]**4)+  
      O[x]**4)

So power series with multivariable can be computed.and Present algorithm is in state of Extension. We need some method which can combine two(number of variable) individual series.first it will calculate series w.r.t to one variable then w.r.t another variable . But this approach only good ,till there are two variable . what if more than two variable are there.

step by solution of above example

>>> s = fps(sin(x+y),x,0)
>>> s1 = fps(s1,y)
>>> s1
( (y - y**3/6 + O[y]**4) + (1- y**2/2 +O[y]**4)*x + ( -y/2 + y**3/12 + O[y]**4) * x**2 +(-1/6 + y**2 /12 +O[y]**4)+  
      O[x]**4)
  • catch the exception of NotImplementedError("multivariate formal power series") in actual code and call the module
    Multi_fps(),where we can handle the multivariate series.In Multi_fps() we will separate the variable and calculate series as I explain above.

So it is time to implement it. and it is easy to implement. [5] explained about multi-series of Inverse Function. presently I am reading this book, and try to understand algorithm to computing these types of series.If it will compatible to implement and not time consuming, then I will implement it but I am sure about it.

Unifying the Series and Expansion

I think that sympy series needs a good user_interface for it's series.User should acknowledged by more information about series.so I plan to implement a new type of series.After Improvement and Implementation , Same FormalPowerseries Algorithm[1] Can used to implement this type of series.new class can be construct to work with all types of series and it's operation . What exactly idea is FormalPowerSeries ---> series(FormalPowerSeries, TaylorSeries,LaurentSeries, Expansion(truncated series with order term O() )).In this Way we will get unified version of series.

Basic class structure: <--un_series-->

  • Points Regarding to implement it-
    • I am not going to implement a whole new series ,just Older formalPowerSeries file can be use for this.I am going
      create a new class named is un_series().where I just handle the properties of different - different type of
      series.so that while computing the series user will get more information about function and it's corresponding series.
    • Methods Which are given below just require some checks while computing series.And It will be easy to implement
      all these.
    • New Series Will be more compatible for user.

For Example: sample Run-


>>> s=sin(x).un_series(x,0,Kind="formal")
>>> s
  x + x**3/6 + x**5/120 + ...
  or it can interpreted as 
>>> summation(a(n)*x**n, (n, 0, oo))
>>> s.coefficients(n=5)
[1, 0, 1/6, 0, 1/120]

Sometime more convenient to obtain asymptotic expansion for some expression or function:

>>> s=sin(x).un_series(x,0,kind="formal")
>>> s.as_expansion(x,6)
x + x**3/6 + x**5/120 + O(x**6)

Operation

Result from series , is infinite series and expansion is till Order term,But it is unreasonably to operate for series only with a few first terms " x + x3/6 + x5/120 " as ordinary polynomial while it is a infinite series. More sense to operate with expansion with Big-O notation:

>>> s.as_expansion(x,6)
x + x**3/6 + x**5/120 + O(x**6)
>>>s1=cos(x).un_series(x,0)
>>> s.as_expansion(x,0) + s1.as_expansion(x,0)
1 + x - x**2/2 - x**3/6 + x**4/24 + x**5/120 + O(x**6)
>>> s1.as_expansion(x,6)**2 
1 - x**2 + x**4/3 + O(x**6)

Big-O signals to user what precision is and to do not forget terms of orders. Also it help for operations: eats high order terms which is not needed for user:

>>> cos(x).as_expansion(x)**2 + 1 + x + x**6 + 777*x**7 + x**8
2 + x - x**2 + x**4/3 + O(x**6)

Note:

All Operation that was part of FormalPowerSeries can also be performed for new un_series( ) . Difference only that while performing operation, un_series( ) should have track about property of function(what is type of series it represent).

Internal and screen representations

User in above examples can do not suspect of internal representation of series. But really .un_series() method yield the expression of special type.

One difference is internal representation: that PowerSeries track coefficients near ( x^k ) terms, but TaylorSeries track coefficients near ( \frac{x^k}{n!} ) ones.

>>> s = sin(x).un_series(x,0)
>>> type(s)
<class 'sympy.un_series.FormalPowerSeries'>

>>>a=a.as_expansion(x,0)
>>> type(a)
sin(x).un_series(x,0)
<class 'sympy.un_series.SeriesExpansion'>


>>> s=sin(x).un_series(x,0,Kind="formal")
>>> s.coefficients(n=5)
[1, 0, 1/6, 0, 1/120]
>>> s=sin(x).un_series(x,0,Kind="taylor")
>>> s
x + x**3/6 + x**5/120 + ...(same as it for formal)
>>> s.coefficients(n=5)
[1, 0, 1, 0, 1]

Operation with un_series same as implemented for formalpowerseries.Although I have to add some operation. These classes track some info about:

>>> s.expression
sin(x)
>>> s.variable
x
>>> s.coefficients(n=3)
[0, 1, 0, 1/6]

>>> a.precision
6

Some properties:

>>> s.is_periodical()
False

>>> s.is_finite()
False

>>> s.to_taylor().coefficients(n=8)
[0, 1, 0, -1, 0, 1, 0, -1]

>>> s.to_taylor().is_periodical()
True

>>> s.to_taylor().period()
4

>>> s.to_taylor().sequence
Sequence(periodical=True, [0, 1, 0, -1])

Some analytic: .isAnalytical .hasLaurentSeries .hasTaylorSeries

Laurent


>>> s = (sin(x)/x**4).series(x)
>>> s
x**(-3) - 1/(6*x) + x/120 - x**3/5040 + x**5/362880 + ...
>>> type(s)
<class 'sympy.series.FormalLaurentSeries'>

>>> s.principalPart()
x**(-3) - 1/(6*x)

>>> (sin(x)/x**4).series(x, kind="taylor")
SeriesError:

ClassicalLaurent (infinite only negative powers)

>>> s = (exp(1/x)).un_series(x, oo)
>>> s
1 + 1/x +  1/(2*x**2) + 1/(6*x**3) + 1/(24*x**4) + 1/(120*x**5) + ...(if power of x is ClassicalLaurent , we can detect while handling it)
>>> type(s)
<class 'sympy.un_series.ClassicalLaurentSeries'>
>>> s.x0 | s.point
oo
>>> s.principalPart().is_finite
False
>>> s.positivePart().is_finite
True

ClassicalLaurent infinite negative end positive powers: >>> s = (exp(x) + exp(1/x)).series(x) ... 1/(6x3) + 1/(2x2) +

1/x + 1 + x + x2/2 + x**3/6 + ...

Multiplication is not allowed in this case.

>>> s*s
SeriesError:  " Multiplication is not allowed "

In this case may be an annulus must be token into consideration. An *annulus* it is the area between two concentric circles

We can distinguish different type of series internally so that user would not worry about , what is result ,and what kind of series any function is going to represent.

Notes:

I have no other plans in summer and I would dedicate 40 hours a week to complete my work according to my timeline.

Most of the references have pseudo-code given in their paper and I would refer them for my help. I plan to be an active SymPy contributor after GSoC period and continue work on the series module in particular.

References:

[1] On computing limits in a symbolic manipulation world – Dominik Gruntz

[2] "Formal Power Series" by Dominik Gruntz and Wolfram Koepf

[3] Power Series in Computer Algebra - Wolfram Koepf

[4] "A New Algorithm Computing for Asymptotic Series" by Dominik Gruntz

[5] ""On the computation of limsups" by Joris van der Hoeven."

[6] http://mathworld.wolfram.com/

[7] Wikipedia - Formal Power Series

Relevant group discussions