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