Coverage for /usr/lib/python3/dist-packages/sympy/concrete/gosper.py: 7%

1"""Gosper's algorithm for hypergeometric summation. """ 

2 

3from sympy.core import S, Dummy, symbols 

4from sympy.polys import Poly, parallel_poly_from_expr, factor 

5from sympy.utilities.iterables import is_sequence 

6 

7 

8def gosper_normal(f, g, n, polys=True): 

9 r""" 

10 Compute the Gosper's normal form of ``f`` and ``g``. 

11 

12 Explanation 

13 =========== 

14 

15 Given relatively prime univariate polynomials ``f`` and ``g``, 

16 rewrite their quotient to a normal form defined as follows: 

17 

18 .. math:: 

19 \frac{f(n)}{g(n)} = Z \cdot \frac{A(n) C(n+1)}{B(n) C(n)} 

20 

21 where ``Z`` is an arbitrary constant and ``A``, ``B``, ``C`` are 

22 monic polynomials in ``n`` with the following properties: 

23 

24 1. `\gcd(A(n), B(n+h)) = 1 \forall h \in \mathbb{N}` 

25 2. `\gcd(B(n), C(n+1)) = 1` 

26 3. `\gcd(A(n), C(n)) = 1` 

27 

28 This normal form, or rational factorization in other words, is a 

29 crucial step in Gosper's algorithm and in solving of difference 

30 equations. It can be also used to decide if two hypergeometric 

31 terms are similar or not. 

32 

33 This procedure will return a tuple containing elements of this 

34 factorization in the form ``(Z*A, B, C)``. 

35 

36 Examples 

37 ======== 

38 

39 >>> from sympy.concrete.gosper import gosper_normal 

40 >>> from sympy.abc import n 

41 

42 >>> gosper_normal(4*n+5, 2*(4*n+1)*(2*n+3), n, polys=False) 

43 (1/4, n + 3/2, n + 1/4) 

44 

45 """ 

46 (p, q), opt = parallel_poly_from_expr( 

47 (f, g), n, field=True, extension=True) 

48 

49 a, A = p.LC(), p.monic() 

50 b, B = q.LC(), q.monic() 

51 

52 C, Z = A.one, a/b 

53 h = Dummy('h') 

54 

55 D = Poly(n + h, n, h, domain=opt.domain) 

56 

57 R = A.resultant(B.compose(D)) 

58 roots = set(R.ground_roots().keys()) 

59 

60 for r in set(roots): 

61 if not r.is_Integer or r < 0: 

62 roots.remove(r) 

63 

64 for i in sorted(roots): 

65 d = A.gcd(B.shift(+i)) 

66 

67 A = A.quo(d) 

68 B = B.quo(d.shift(-i)) 

69 

70 for j in range(1, i + 1): 

71 C *= d.shift(-j) 

72 

73 A = A.mul_ground(Z) 

74 

75 if not polys: 

76 A = A.as_expr() 

77 B = B.as_expr() 

78 C = C.as_expr() 

79 

80 return A, B, C 

81 

82 

83def gosper_term(f, n): 

84 r""" 

85 Compute Gosper's hypergeometric term for ``f``. 

86 

87 Explanation 

88 =========== 

89 

90 Suppose ``f`` is a hypergeometric term such that: 

91 

92 .. math:: 

93 s_n = \sum_{k=0}^{n-1} f_k 

94 

95 and `f_k` does not depend on `n`. Returns a hypergeometric 

96 term `g_n` such that `g_{n+1} - g_n = f_n`. 

97 

98 Examples 

99 ======== 

100 

101 >>> from sympy.concrete.gosper import gosper_term 

102 >>> from sympy import factorial 

103 >>> from sympy.abc import n 

104 

105 >>> gosper_term((4*n + 1)*factorial(n)/factorial(2*n + 1), n) 

106 (-n - 1/2)/(n + 1/4) 

107 

108 """ 

109 from sympy.simplify import hypersimp 

110 r = hypersimp(f, n) 

111 

112 if r is None: 

113 return None # 'f' is *not* a hypergeometric term 

114 

115 p, q = r.as_numer_denom() 

116 

117 A, B, C = gosper_normal(p, q, n) 

118 B = B.shift(-1) 

119 

120 N = S(A.degree()) 

121 M = S(B.degree()) 

122 K = S(C.degree()) 

123 

124 if (N != M) or (A.LC() != B.LC()): 

125 D = {K - max(N, M)} 

126 elif not N: 

127 D = {K - N + 1, S.Zero} 

128 else: 

129 D = {K - N + 1, (B.nth(N - 1) - A.nth(N - 1))/A.LC()} 

130 

131 for d in set(D): 

132 if not d.is_Integer or d < 0: 

133 D.remove(d) 

134 

135 if not D: 

136 return None # 'f(n)' is *not* Gosper-summable 

137 

138 d = max(D) 

139 

140 coeffs = symbols('c:%s' % (d + 1), cls=Dummy) 

141 domain = A.get_domain().inject(*coeffs) 

142 

143 x = Poly(coeffs, n, domain=domain) 

144 H = A*x.shift(1) - B*x - C 

145 

146 from sympy.solvers.solvers import solve 

147 solution = solve(H.coeffs(), coeffs) 

148 

149 if solution is None: 

150 return None # 'f(n)' is *not* Gosper-summable 

151 

152 x = x.as_expr().subs(solution) 

153 

154 for coeff in coeffs: 

155 if coeff not in solution: 

156 x = x.subs(coeff, 0) 

157 

158 if x.is_zero: 

159 return None # 'f(n)' is *not* Gosper-summable 

160 else: 

161 return B.as_expr()*x/C.as_expr() 

162 

163 

164def gosper_sum(f, k): 

165 r""" 

166 Gosper's hypergeometric summation algorithm. 

167 

168 Explanation 

169 =========== 

170 

171 Given a hypergeometric term ``f`` such that: 

172 

173 .. math :: 

174 s_n = \sum_{k=0}^{n-1} f_k 

175 

176 and `f(n)` does not depend on `n`, returns `g_{n} - g(0)` where 

177 `g_{n+1} - g_n = f_n`, or ``None`` if `s_n` cannot be expressed 

178 in closed form as a sum of hypergeometric terms. 

179 

180 Examples 

181 ======== 

182 

183 >>> from sympy.concrete.gosper import gosper_sum 

184 >>> from sympy import factorial 

185 >>> from sympy.abc import n, k 

186 

187 >>> f = (4*k + 1)*factorial(k)/factorial(2*k + 1) 

188 >>> gosper_sum(f, (k, 0, n)) 

189 (-factorial(n) + 2*factorial(2*n + 1))/factorial(2*n + 1) 

190 >>> _.subs(n, 2) == sum(f.subs(k, i) for i in [0, 1, 2]) 

191 True 

192 >>> gosper_sum(f, (k, 3, n)) 

193 (-60*factorial(n) + factorial(2*n + 1))/(60*factorial(2*n + 1)) 

194 >>> _.subs(n, 5) == sum(f.subs(k, i) for i in [3, 4, 5]) 

195 True 

196 

197 References 

198 ========== 

199 

200 .. [1] Marko Petkovsek, Herbert S. Wilf, Doron Zeilberger, A = B, 

201 AK Peters, Ltd., Wellesley, MA, USA, 1997, pp. 73--100 

202 

203 """ 

204 indefinite = False 

205 

206 if is_sequence(k): 

207 k, a, b = k 

208 else: 

209 indefinite = True 

210 

211 g = gosper_term(f, k) 

212 

213 if g is None: 

214 return None 

215 

216 if indefinite: 

217 result = f*g 

218 else: 

219 result = (f*(g + 1)).subs(k, b) - (f*g).subs(k, a) 

220 

221 if result is S.NaN: 

222 try: 

223 result = (f*(g + 1)).limit(k, b) - (f*g).limit(k, a) 

224 except NotImplementedError: 

225 result = None 

226 

227 return factor(result)