1from sympy import Integer
2from sympy.core import Symbol
3from sympy.utilities import public
4
5@public
6def approximants(l, X=Symbol('x'), simplify=False):
7    """
8    Return a generator for consecutive Pade approximants for a series.
9    It can also be used for computing the rational generating function of a
10    series when possible, since the last approximant returned by the generator
11    will be the generating function (if any).
12
13    Explanation
14    ===========
15
16    The input list can contain more complex expressions than integer or rational
17    numbers; symbols may also be involved in the computation. An example below
18    show how to compute the generating function of the whole Pascal triangle.
19
20    The generator can be asked to apply the sympy.simplify function on each
21    generated term, which will make the computation slower; however it may be
22    useful when symbols are involved in the expressions.
23
24    Examples
25    ========
26
27    >>> from sympy.series import approximants
28    >>> from sympy import lucas, fibonacci, symbols, binomial
29    >>> g = [lucas(k) for k in range(16)]
30    >>> [e for e in approximants(g)]
31    [2, -4/(x - 2), (5*x - 2)/(3*x - 1), (x - 2)/(x**2 + x - 1)]
32
33    >>> h = [fibonacci(k) for k in range(16)]
34    >>> [e for e in approximants(h)]
35    [x, -x/(x - 1), (x**2 - x)/(2*x - 1), -x/(x**2 + x - 1)]
36
37    >>> x, t = symbols("x,t")
38    >>> p=[sum(binomial(k,i)*x**i for i in range(k+1)) for k in range(16)]
39    >>> y = approximants(p, t)
40    >>> for k in range(3): print(next(y))
41    1
42    (x + 1)/((-x - 1)*(t*(x + 1) + (x + 1)/(-x - 1)))
43    nan
44
45    >>> y = approximants(p, t, simplify=True)
46    >>> for k in range(3): print(next(y))
47    1
48    -1/(t*(x + 1) - 1)
49    nan
50
51    See Also
52    ========
53
54    See function sympy.concrete.guess.guess_generating_function_rational and
55    function mpmath.pade
56
57    """
58    p1, q1 = [Integer(1)], [Integer(0)]
59    p2, q2 = [Integer(0)], [Integer(1)]
60    while len(l):
61        b = 0
62        while l[b]==0:
63            b += 1
64            if b == len(l):
65                return
66        m = [Integer(1)/l[b]]
67        for k in range(b+1, len(l)):
68            s = 0
69            for j in range(b, k):
70                s -= l[j+1] * m[b-j-1]
71            m.append(s/l[b])
72        l = m
73        a, l[0] = l[0], 0
74        p = [0] * max(len(p2), b+len(p1))
75        q = [0] * max(len(q2), b+len(q1))
76        for k in range(len(p2)):
77            p[k] = a*p2[k]
78        for k in range(b, b+len(p1)):
79            p[k] += p1[k-b]
80        for k in range(len(q2)):
81            q[k] = a*q2[k]
82        for k in range(b, b+len(q1)):
83            q[k] += q1[k-b]
84        while p[-1]==0: p.pop()
85        while q[-1]==0: q.pop()
86        p1, p2 = p2, p
87        q1, q2 = q2, q
88
89        # yield result
90        from sympy import denom, lcm, simplify as simp
91        c = 1
92        for x in p:
93            c = lcm(c, denom(x))
94        for x in q:
95            c = lcm(c, denom(x))
96        out = ( sum(c*e*X**k for k, e in enumerate(p))
97              / sum(c*e*X**k for k, e in enumerate(q)) )
98        if simplify: yield(simp(out))
99        else: yield out
100    return
101