Coverage for /usr/lib/python3/dist-packages/sympy/simplify/hyperexpand.py: 11%
1225 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2Expand Hypergeometric (and Meijer G) functions into named
3special functions.
5The algorithm for doing this uses a collection of lookup tables of
6hypergeometric functions, and various of their properties, to expand
7many hypergeometric functions in terms of special functions.
9It is based on the following paper:
10 Kelly B. Roach. Meijer G Function Representations.
11 In: Proceedings of the 1997 International Symposium on Symbolic and
12 Algebraic Computation, pages 205-211, New York, 1997. ACM.
14It is described in great(er) detail in the Sphinx documentation.
15"""
16# SUMMARY OF EXTENSIONS FOR MEIJER G FUNCTIONS
17#
18# o z**rho G(ap, bq; z) = G(ap + rho, bq + rho; z)
19#
20# o denote z*d/dz by D
21#
22# o It is helpful to keep in mind that ap and bq play essentially symmetric
23# roles: G(1/z) has slightly altered parameters, with ap and bq interchanged.
24#
25# o There are four shift operators:
26# A_J = b_J - D, J = 1, ..., n
27# B_J = 1 - a_j + D, J = 1, ..., m
28# C_J = -b_J + D, J = m+1, ..., q
29# D_J = a_J - 1 - D, J = n+1, ..., p
30#
31# A_J, C_J increment b_J
32# B_J, D_J decrement a_J
33#
34# o The corresponding four inverse-shift operators are defined if there
35# is no cancellation. Thus e.g. an index a_J (upper or lower) can be
36# incremented if a_J != b_i for i = 1, ..., q.
37#
38# o Order reduction: if b_j - a_i is a non-negative integer, where
39# j <= m and i > n, the corresponding quotient of gamma functions reduces
40# to a polynomial. Hence the G function can be expressed using a G-function
41# of lower order.
42# Similarly if j > m and i <= n.
43#
44# Secondly, there are paired index theorems [Adamchik, The evaluation of
45# integrals of Bessel functions via G-function identities]. Suppose there
46# are three parameters a, b, c, where a is an a_i, i <= n, b is a b_j,
47# j <= m and c is a denominator parameter (i.e. a_i, i > n or b_j, j > m).
48# Suppose further all three differ by integers.
49# Then the order can be reduced.
50# TODO work this out in detail.
51#
52# o An index quadruple is called suitable if its order cannot be reduced.
53# If there exists a sequence of shift operators transforming one index
54# quadruple into another, we say one is reachable from the other.
55#
56# o Deciding if one index quadruple is reachable from another is tricky. For
57# this reason, we use hand-built routines to match and instantiate formulas.
58#
59from collections import defaultdict
60from itertools import product
61from functools import reduce
62from math import prod
64from sympy import SYMPY_DEBUG
65from sympy.core import (S, Dummy, symbols, sympify, Tuple, expand, I, pi, Mul,
66 EulerGamma, oo, zoo, expand_func, Add, nan, Expr, Rational)
67from sympy.core.mod import Mod
68from sympy.core.sorting import default_sort_key
69from sympy.functions import (exp, sqrt, root, log, lowergamma, cos,
70 besseli, gamma, uppergamma, expint, erf, sin, besselj, Ei, Ci, Si, Shi,
71 sinh, cosh, Chi, fresnels, fresnelc, polar_lift, exp_polar, floor, ceiling,
72 rf, factorial, lerchphi, Piecewise, re, elliptic_k, elliptic_e)
73from sympy.functions.elementary.complexes import polarify, unpolarify
74from sympy.functions.special.hyper import (hyper, HyperRep_atanh,
75 HyperRep_power1, HyperRep_power2, HyperRep_log1, HyperRep_asin1,
76 HyperRep_asin2, HyperRep_sqrts1, HyperRep_sqrts2, HyperRep_log2,
77 HyperRep_cosasin, HyperRep_sinasin, meijerg)
78from sympy.matrices import Matrix, eye, zeros
79from sympy.polys import apart, poly, Poly
80from sympy.series import residue
81from sympy.simplify.powsimp import powdenest
82from sympy.utilities.iterables import sift
84# function to define "buckets"
85def _mod1(x):
86 # TODO see if this can work as Mod(x, 1); this will require
87 # different handling of the "buckets" since these need to
88 # be sorted and that fails when there is a mixture of
89 # integers and expressions with parameters. With the current
90 # Mod behavior, Mod(k, 1) == Mod(1, 1) == 0 if k is an integer.
91 # Although the sorting can be done with Basic.compare, this may
92 # still require different handling of the sorted buckets.
93 if x.is_Number:
94 return Mod(x, 1)
95 c, x = x.as_coeff_Add()
96 return Mod(c, 1) + x
99# leave add formulae at the top for easy reference
100def add_formulae(formulae):
101 """ Create our knowledge base. """
102 a, b, c, z = symbols('a b c, z', cls=Dummy)
104 def add(ap, bq, res):
105 func = Hyper_Function(ap, bq)
106 formulae.append(Formula(func, z, res, (a, b, c)))
108 def addb(ap, bq, B, C, M):
109 func = Hyper_Function(ap, bq)
110 formulae.append(Formula(func, z, None, (a, b, c), B, C, M))
112 # Luke, Y. L. (1969), The Special Functions and Their Approximations,
113 # Volume 1, section 6.2
115 # 0F0
116 add((), (), exp(z))
118 # 1F0
119 add((a, ), (), HyperRep_power1(-a, z))
121 # 2F1
122 addb((a, a - S.Half), (2*a, ),
123 Matrix([HyperRep_power2(a, z),
124 HyperRep_power2(a + S.Half, z)/2]),
125 Matrix([[1, 0]]),
126 Matrix([[(a - S.Half)*z/(1 - z), (S.Half - a)*z/(1 - z)],
127 [a/(1 - z), a*(z - 2)/(1 - z)]]))
128 addb((1, 1), (2, ),
129 Matrix([HyperRep_log1(z), 1]), Matrix([[-1/z, 0]]),
130 Matrix([[0, z/(z - 1)], [0, 0]]))
131 addb((S.Half, 1), (S('3/2'), ),
132 Matrix([HyperRep_atanh(z), 1]),
133 Matrix([[1, 0]]),
134 Matrix([[Rational(-1, 2), 1/(1 - z)/2], [0, 0]]))
135 addb((S.Half, S.Half), (S('3/2'), ),
136 Matrix([HyperRep_asin1(z), HyperRep_power1(Rational(-1, 2), z)]),
137 Matrix([[1, 0]]),
138 Matrix([[Rational(-1, 2), S.Half], [0, z/(1 - z)/2]]))
139 addb((a, S.Half + a), (S.Half, ),
140 Matrix([HyperRep_sqrts1(-a, z), -HyperRep_sqrts2(-a - S.Half, z)]),
141 Matrix([[1, 0]]),
142 Matrix([[0, -a],
143 [z*(-2*a - 1)/2/(1 - z), S.Half - z*(-2*a - 1)/(1 - z)]]))
145 # A. P. Prudnikov, Yu. A. Brychkov and O. I. Marichev (1990).
146 # Integrals and Series: More Special Functions, Vol. 3,.
147 # Gordon and Breach Science Publisher
148 addb([a, -a], [S.Half],
149 Matrix([HyperRep_cosasin(a, z), HyperRep_sinasin(a, z)]),
150 Matrix([[1, 0]]),
151 Matrix([[0, -a], [a*z/(1 - z), 1/(1 - z)/2]]))
152 addb([1, 1], [3*S.Half],
153 Matrix([HyperRep_asin2(z), 1]), Matrix([[1, 0]]),
154 Matrix([[(z - S.Half)/(1 - z), 1/(1 - z)/2], [0, 0]]))
156 # Complete elliptic integrals K(z) and E(z), both a 2F1 function
157 addb([S.Half, S.Half], [S.One],
158 Matrix([elliptic_k(z), elliptic_e(z)]),
159 Matrix([[2/pi, 0]]),
160 Matrix([[Rational(-1, 2), -1/(2*z-2)],
161 [Rational(-1, 2), S.Half]]))
162 addb([Rational(-1, 2), S.Half], [S.One],
163 Matrix([elliptic_k(z), elliptic_e(z)]),
164 Matrix([[0, 2/pi]]),
165 Matrix([[Rational(-1, 2), -1/(2*z-2)],
166 [Rational(-1, 2), S.Half]]))
168 # 3F2
169 addb([Rational(-1, 2), 1, 1], [S.Half, 2],
170 Matrix([z*HyperRep_atanh(z), HyperRep_log1(z), 1]),
171 Matrix([[Rational(-2, 3), -S.One/(3*z), Rational(2, 3)]]),
172 Matrix([[S.Half, 0, z/(1 - z)/2],
173 [0, 0, z/(z - 1)],
174 [0, 0, 0]]))
175 # actually the formula for 3/2 is much nicer ...
176 addb([Rational(-1, 2), 1, 1], [2, 2],
177 Matrix([HyperRep_power1(S.Half, z), HyperRep_log2(z), 1]),
178 Matrix([[Rational(4, 9) - 16/(9*z), 4/(3*z), 16/(9*z)]]),
179 Matrix([[z/2/(z - 1), 0, 0], [1/(2*(z - 1)), 0, S.Half], [0, 0, 0]]))
181 # 1F1
182 addb([1], [b], Matrix([z**(1 - b) * exp(z) * lowergamma(b - 1, z), 1]),
183 Matrix([[b - 1, 0]]), Matrix([[1 - b + z, 1], [0, 0]]))
184 addb([a], [2*a],
185 Matrix([z**(S.Half - a)*exp(z/2)*besseli(a - S.Half, z/2)
186 * gamma(a + S.Half)/4**(S.Half - a),
187 z**(S.Half - a)*exp(z/2)*besseli(a + S.Half, z/2)
188 * gamma(a + S.Half)/4**(S.Half - a)]),
189 Matrix([[1, 0]]),
190 Matrix([[z/2, z/2], [z/2, (z/2 - 2*a)]]))
191 mz = polar_lift(-1)*z
192 addb([a], [a + 1],
193 Matrix([mz**(-a)*a*lowergamma(a, mz), a*exp(z)]),
194 Matrix([[1, 0]]),
195 Matrix([[-a, 1], [0, z]]))
196 # This one is redundant.
197 add([Rational(-1, 2)], [S.Half], exp(z) - sqrt(pi*z)*(-I)*erf(I*sqrt(z)))
199 # Added to get nice results for Laplace transform of Fresnel functions
200 # https://functions.wolfram.com/07.22.03.6437.01
201 # Basic rule
202 #add([1], [Rational(3, 4), Rational(5, 4)],
203 # sqrt(pi) * (cos(2*sqrt(polar_lift(-1)*z))*fresnelc(2*root(polar_lift(-1)*z,4)/sqrt(pi)) +
204 # sin(2*sqrt(polar_lift(-1)*z))*fresnels(2*root(polar_lift(-1)*z,4)/sqrt(pi)))
205 # / (2*root(polar_lift(-1)*z,4)))
206 # Manually tuned rule
207 addb([1], [Rational(3, 4), Rational(5, 4)],
208 Matrix([ sqrt(pi)*(I*sinh(2*sqrt(z))*fresnels(2*root(z, 4)*exp(I*pi/4)/sqrt(pi))
209 + cosh(2*sqrt(z))*fresnelc(2*root(z, 4)*exp(I*pi/4)/sqrt(pi)))
210 * exp(-I*pi/4)/(2*root(z, 4)),
211 sqrt(pi)*root(z, 4)*(sinh(2*sqrt(z))*fresnelc(2*root(z, 4)*exp(I*pi/4)/sqrt(pi))
212 + I*cosh(2*sqrt(z))*fresnels(2*root(z, 4)*exp(I*pi/4)/sqrt(pi)))
213 *exp(-I*pi/4)/2,
214 1 ]),
215 Matrix([[1, 0, 0]]),
216 Matrix([[Rational(-1, 4), 1, Rational(1, 4)],
217 [ z, Rational(1, 4), 0],
218 [ 0, 0, 0]]))
220 # 2F2
221 addb([S.Half, a], [Rational(3, 2), a + 1],
222 Matrix([a/(2*a - 1)*(-I)*sqrt(pi/z)*erf(I*sqrt(z)),
223 a/(2*a - 1)*(polar_lift(-1)*z)**(-a)*
224 lowergamma(a, polar_lift(-1)*z),
225 a/(2*a - 1)*exp(z)]),
226 Matrix([[1, -1, 0]]),
227 Matrix([[Rational(-1, 2), 0, 1], [0, -a, 1], [0, 0, z]]))
228 # We make a "basis" of four functions instead of three, and give EulerGamma
229 # an extra slot (it could just be a coefficient to 1). The advantage is
230 # that this way Polys will not see multivariate polynomials (it treats
231 # EulerGamma as an indeterminate), which is *way* faster.
232 addb([1, 1], [2, 2],
233 Matrix([Ei(z) - log(z), exp(z), 1, EulerGamma]),
234 Matrix([[1/z, 0, 0, -1/z]]),
235 Matrix([[0, 1, -1, 0], [0, z, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]))
237 # 0F1
238 add((), (S.Half, ), cosh(2*sqrt(z)))
239 addb([], [b],
240 Matrix([gamma(b)*z**((1 - b)/2)*besseli(b - 1, 2*sqrt(z)),
241 gamma(b)*z**(1 - b/2)*besseli(b, 2*sqrt(z))]),
242 Matrix([[1, 0]]), Matrix([[0, 1], [z, (1 - b)]]))
244 # 0F3
245 x = 4*z**Rational(1, 4)
247 def fp(a, z):
248 return besseli(a, x) + besselj(a, x)
250 def fm(a, z):
251 return besseli(a, x) - besselj(a, x)
253 # TODO branching
254 addb([], [S.Half, a, a + S.Half],
255 Matrix([fp(2*a - 1, z), fm(2*a, z)*z**Rational(1, 4),
256 fm(2*a - 1, z)*sqrt(z), fp(2*a, z)*z**Rational(3, 4)])
257 * 2**(-2*a)*gamma(2*a)*z**((1 - 2*a)/4),
258 Matrix([[1, 0, 0, 0]]),
259 Matrix([[0, 1, 0, 0],
260 [0, S.Half - a, 1, 0],
261 [0, 0, S.Half, 1],
262 [z, 0, 0, 1 - a]]))
263 x = 2*(4*z)**Rational(1, 4)*exp_polar(I*pi/4)
264 addb([], [a, a + S.Half, 2*a],
265 (2*sqrt(polar_lift(-1)*z))**(1 - 2*a)*gamma(2*a)**2 *
266 Matrix([besselj(2*a - 1, x)*besseli(2*a - 1, x),
267 x*(besseli(2*a, x)*besselj(2*a - 1, x)
268 - besseli(2*a - 1, x)*besselj(2*a, x)),
269 x**2*besseli(2*a, x)*besselj(2*a, x),
270 x**3*(besseli(2*a, x)*besselj(2*a - 1, x)
271 + besseli(2*a - 1, x)*besselj(2*a, x))]),
272 Matrix([[1, 0, 0, 0]]),
273 Matrix([[0, Rational(1, 4), 0, 0],
274 [0, (1 - 2*a)/2, Rational(-1, 2), 0],
275 [0, 0, 1 - 2*a, Rational(1, 4)],
276 [-32*z, 0, 0, 1 - a]]))
278 # 1F2
279 addb([a], [a - S.Half, 2*a],
280 Matrix([z**(S.Half - a)*besseli(a - S.Half, sqrt(z))**2,
281 z**(1 - a)*besseli(a - S.Half, sqrt(z))
282 *besseli(a - Rational(3, 2), sqrt(z)),
283 z**(Rational(3, 2) - a)*besseli(a - Rational(3, 2), sqrt(z))**2]),
284 Matrix([[-gamma(a + S.Half)**2/4**(S.Half - a),
285 2*gamma(a - S.Half)*gamma(a + S.Half)/4**(1 - a),
286 0]]),
287 Matrix([[1 - 2*a, 1, 0], [z/2, S.Half - a, S.Half], [0, z, 0]]))
288 addb([S.Half], [b, 2 - b],
289 pi*(1 - b)/sin(pi*b)*
290 Matrix([besseli(1 - b, sqrt(z))*besseli(b - 1, sqrt(z)),
291 sqrt(z)*(besseli(-b, sqrt(z))*besseli(b - 1, sqrt(z))
292 + besseli(1 - b, sqrt(z))*besseli(b, sqrt(z))),
293 besseli(-b, sqrt(z))*besseli(b, sqrt(z))]),
294 Matrix([[1, 0, 0]]),
295 Matrix([[b - 1, S.Half, 0],
296 [z, 0, z],
297 [0, S.Half, -b]]))
298 addb([S.Half], [Rational(3, 2), Rational(3, 2)],
299 Matrix([Shi(2*sqrt(z))/2/sqrt(z), sinh(2*sqrt(z))/2/sqrt(z),
300 cosh(2*sqrt(z))]),
301 Matrix([[1, 0, 0]]),
302 Matrix([[Rational(-1, 2), S.Half, 0], [0, Rational(-1, 2), S.Half], [0, 2*z, 0]]))
304 # FresnelS
305 # Basic rule
306 #add([Rational(3, 4)], [Rational(3, 2),Rational(7, 4)], 6*fresnels( exp(pi*I/4)*root(z,4)*2/sqrt(pi) ) / ( pi * (exp(pi*I/4)*root(z,4)*2/sqrt(pi))**3 ) )
307 # Manually tuned rule
308 addb([Rational(3, 4)], [Rational(3, 2), Rational(7, 4)],
309 Matrix(
310 [ fresnels(
311 exp(
312 pi*I/4)*root(
313 z, 4)*2/sqrt(
314 pi) ) / (
315 pi * (exp(pi*I/4)*root(z, 4)*2/sqrt(pi))**3 ),
316 sinh(2*sqrt(z))/sqrt(z),
317 cosh(2*sqrt(z)) ]),
318 Matrix([[6, 0, 0]]),
319 Matrix([[Rational(-3, 4), Rational(1, 16), 0],
320 [ 0, Rational(-1, 2), 1],
321 [ 0, z, 0]]))
323 # FresnelC
324 # Basic rule
325 #add([Rational(1, 4)], [S.Half,Rational(5, 4)], fresnelc( exp(pi*I/4)*root(z,4)*2/sqrt(pi) ) / ( exp(pi*I/4)*root(z,4)*2/sqrt(pi) ) )
326 # Manually tuned rule
327 addb([Rational(1, 4)], [S.Half, Rational(5, 4)],
328 Matrix(
329 [ sqrt(
330 pi)*exp(
331 -I*pi/4)*fresnelc(
332 2*root(z, 4)*exp(I*pi/4)/sqrt(pi))/(2*root(z, 4)),
333 cosh(2*sqrt(z)),
334 sinh(2*sqrt(z))*sqrt(z) ]),
335 Matrix([[1, 0, 0]]),
336 Matrix([[Rational(-1, 4), Rational(1, 4), 0 ],
337 [ 0, 0, 1 ],
338 [ 0, z, S.Half]]))
340 # 2F3
341 # XXX with this five-parameter formula is pretty slow with the current
342 # Formula.find_instantiations (creates 2!*3!*3**(2+3) ~ 3000
343 # instantiations ... But it's not too bad.
344 addb([a, a + S.Half], [2*a, b, 2*a - b + 1],
345 gamma(b)*gamma(2*a - b + 1) * (sqrt(z)/2)**(1 - 2*a) *
346 Matrix([besseli(b - 1, sqrt(z))*besseli(2*a - b, sqrt(z)),
347 sqrt(z)*besseli(b, sqrt(z))*besseli(2*a - b, sqrt(z)),
348 sqrt(z)*besseli(b - 1, sqrt(z))*besseli(2*a - b + 1, sqrt(z)),
349 besseli(b, sqrt(z))*besseli(2*a - b + 1, sqrt(z))]),
350 Matrix([[1, 0, 0, 0]]),
351 Matrix([[0, S.Half, S.Half, 0],
352 [z/2, 1 - b, 0, z/2],
353 [z/2, 0, b - 2*a, z/2],
354 [0, S.Half, S.Half, -2*a]]))
355 # (C/f above comment about eulergamma in the basis).
356 addb([1, 1], [2, 2, Rational(3, 2)],
357 Matrix([Chi(2*sqrt(z)) - log(2*sqrt(z)),
358 cosh(2*sqrt(z)), sqrt(z)*sinh(2*sqrt(z)), 1, EulerGamma]),
359 Matrix([[1/z, 0, 0, 0, -1/z]]),
360 Matrix([[0, S.Half, 0, Rational(-1, 2), 0],
361 [0, 0, 1, 0, 0],
362 [0, z, S.Half, 0, 0],
363 [0, 0, 0, 0, 0],
364 [0, 0, 0, 0, 0]]))
366 # 3F3
367 # This is rule: https://functions.wolfram.com/07.31.03.0134.01
368 # Initial reason to add it was a nice solution for
369 # integrate(erf(a*z)/z**2, z) and same for erfc and erfi.
370 # Basic rule
371 # add([1, 1, a], [2, 2, a+1], (a/(z*(a-1)**2)) *
372 # (1 - (-z)**(1-a) * (gamma(a) - uppergamma(a,-z))
373 # - (a-1) * (EulerGamma + uppergamma(0,-z) + log(-z))
374 # - exp(z)))
375 # Manually tuned rule
376 addb([1, 1, a], [2, 2, a+1],
377 Matrix([a*(log(-z) + expint(1, -z) + EulerGamma)/(z*(a**2 - 2*a + 1)),
378 a*(-z)**(-a)*(gamma(a) - uppergamma(a, -z))/(a - 1)**2,
379 a*exp(z)/(a**2 - 2*a + 1),
380 a/(z*(a**2 - 2*a + 1))]),
381 Matrix([[1-a, 1, -1/z, 1]]),
382 Matrix([[-1,0,-1/z,1],
383 [0,-a,1,0],
384 [0,0,z,0],
385 [0,0,0,-1]]))
388def add_meijerg_formulae(formulae):
389 a, b, c, z = list(map(Dummy, 'abcz'))
390 rho = Dummy('rho')
392 def add(an, ap, bm, bq, B, C, M, matcher):
393 formulae.append(MeijerFormula(an, ap, bm, bq, z, [a, b, c, rho],
394 B, C, M, matcher))
396 def detect_uppergamma(func):
397 x = func.an[0]
398 y, z = func.bm
399 swapped = False
400 if not _mod1((x - y).simplify()):
401 swapped = True
402 (y, z) = (z, y)
403 if _mod1((x - z).simplify()) or x - z > 0:
404 return None
405 l = [y, x]
406 if swapped:
407 l = [x, y]
408 return {rho: y, a: x - y}, G_Function([x], [], l, [])
410 add([a + rho], [], [rho, a + rho], [],
411 Matrix([gamma(1 - a)*z**rho*exp(z)*uppergamma(a, z),
412 gamma(1 - a)*z**(a + rho)]),
413 Matrix([[1, 0]]),
414 Matrix([[rho + z, -1], [0, a + rho]]),
415 detect_uppergamma)
417 def detect_3113(func):
418 """https://functions.wolfram.com/07.34.03.0984.01"""
419 x = func.an[0]
420 u, v, w = func.bm
421 if _mod1((u - v).simplify()) == 0:
422 if _mod1((v - w).simplify()) == 0:
423 return
424 sig = (S.Half, S.Half, S.Zero)
425 x1, x2, y = u, v, w
426 else:
427 if _mod1((x - u).simplify()) == 0:
428 sig = (S.Half, S.Zero, S.Half)
429 x1, y, x2 = u, v, w
430 else:
431 sig = (S.Zero, S.Half, S.Half)
432 y, x1, x2 = u, v, w
434 if (_mod1((x - x1).simplify()) != 0 or
435 _mod1((x - x2).simplify()) != 0 or
436 _mod1((x - y).simplify()) != S.Half or
437 x - x1 > 0 or x - x2 > 0):
438 return
440 return {a: x}, G_Function([x], [], [x - S.Half + t for t in sig], [])
442 s = sin(2*sqrt(z))
443 c_ = cos(2*sqrt(z))
444 S_ = Si(2*sqrt(z)) - pi/2
445 C = Ci(2*sqrt(z))
446 add([a], [], [a, a, a - S.Half], [],
447 Matrix([sqrt(pi)*z**(a - S.Half)*(c_*S_ - s*C),
448 sqrt(pi)*z**a*(s*S_ + c_*C),
449 sqrt(pi)*z**a]),
450 Matrix([[-2, 0, 0]]),
451 Matrix([[a - S.Half, -1, 0], [z, a, S.Half], [0, 0, a]]),
452 detect_3113)
455def make_simp(z):
456 """ Create a function that simplifies rational functions in ``z``. """
458 def simp(expr):
459 """ Efficiently simplify the rational function ``expr``. """
460 numer, denom = expr.as_numer_denom()
461 numer = numer.expand()
462 # denom = denom.expand() # is this needed?
463 c, numer, denom = poly(numer, z).cancel(poly(denom, z))
464 return c * numer.as_expr() / denom.as_expr()
466 return simp
469def debug(*args):
470 if SYMPY_DEBUG:
471 for a in args:
472 print(a, end="")
473 print()
476class Hyper_Function(Expr):
477 """ A generalized hypergeometric function. """
479 def __new__(cls, ap, bq):
480 obj = super().__new__(cls)
481 obj.ap = Tuple(*list(map(expand, ap)))
482 obj.bq = Tuple(*list(map(expand, bq)))
483 return obj
485 @property
486 def args(self):
487 return (self.ap, self.bq)
489 @property
490 def sizes(self):
491 return (len(self.ap), len(self.bq))
493 @property
494 def gamma(self):
495 """
496 Number of upper parameters that are negative integers
498 This is a transformation invariant.
499 """
500 return sum(bool(x.is_integer and x.is_negative) for x in self.ap)
502 def _hashable_content(self):
503 return super()._hashable_content() + (self.ap,
504 self.bq)
506 def __call__(self, arg):
507 return hyper(self.ap, self.bq, arg)
509 def build_invariants(self):
510 """
511 Compute the invariant vector.
513 Explanation
514 ===========
516 The invariant vector is:
517 (gamma, ((s1, n1), ..., (sk, nk)), ((t1, m1), ..., (tr, mr)))
518 where gamma is the number of integer a < 0,
519 s1 < ... < sk
520 nl is the number of parameters a_i congruent to sl mod 1
521 t1 < ... < tr
522 ml is the number of parameters b_i congruent to tl mod 1
524 If the index pair contains parameters, then this is not truly an
525 invariant, since the parameters cannot be sorted uniquely mod1.
527 Examples
528 ========
530 >>> from sympy.simplify.hyperexpand import Hyper_Function
531 >>> from sympy import S
532 >>> ap = (S.Half, S.One/3, S(-1)/2, -2)
533 >>> bq = (1, 2)
535 Here gamma = 1,
536 k = 3, s1 = 0, s2 = 1/3, s3 = 1/2
537 n1 = 1, n2 = 1, n2 = 2
538 r = 1, t1 = 0
539 m1 = 2:
541 >>> Hyper_Function(ap, bq).build_invariants()
542 (1, ((0, 1), (1/3, 1), (1/2, 2)), ((0, 2),))
543 """
544 abuckets, bbuckets = sift(self.ap, _mod1), sift(self.bq, _mod1)
546 def tr(bucket):
547 bucket = list(bucket.items())
548 if not any(isinstance(x[0], Mod) for x in bucket):
549 bucket.sort(key=lambda x: default_sort_key(x[0]))
550 bucket = tuple([(mod, len(values)) for mod, values in bucket if
551 values])
552 return bucket
554 return (self.gamma, tr(abuckets), tr(bbuckets))
556 def difficulty(self, func):
557 """ Estimate how many steps it takes to reach ``func`` from self.
558 Return -1 if impossible. """
559 if self.gamma != func.gamma:
560 return -1
561 oabuckets, obbuckets, abuckets, bbuckets = [sift(params, _mod1) for
562 params in (self.ap, self.bq, func.ap, func.bq)]
564 diff = 0
565 for bucket, obucket in [(abuckets, oabuckets), (bbuckets, obbuckets)]:
566 for mod in set(list(bucket.keys()) + list(obucket.keys())):
567 if (mod not in bucket) or (mod not in obucket) \
568 or len(bucket[mod]) != len(obucket[mod]):
569 return -1
570 l1 = list(bucket[mod])
571 l2 = list(obucket[mod])
572 l1.sort()
573 l2.sort()
574 for i, j in zip(l1, l2):
575 diff += abs(i - j)
577 return diff
579 def _is_suitable_origin(self):
580 """
581 Decide if ``self`` is a suitable origin.
583 Explanation
584 ===========
586 A function is a suitable origin iff:
587 * none of the ai equals bj + n, with n a non-negative integer
588 * none of the ai is zero
589 * none of the bj is a non-positive integer
591 Note that this gives meaningful results only when none of the indices
592 are symbolic.
594 """
595 for a in self.ap:
596 for b in self.bq:
597 if (a - b).is_integer and (a - b).is_negative is False:
598 return False
599 for a in self.ap:
600 if a == 0:
601 return False
602 for b in self.bq:
603 if b.is_integer and b.is_nonpositive:
604 return False
605 return True
608class G_Function(Expr):
609 """ A Meijer G-function. """
611 def __new__(cls, an, ap, bm, bq):
612 obj = super().__new__(cls)
613 obj.an = Tuple(*list(map(expand, an)))
614 obj.ap = Tuple(*list(map(expand, ap)))
615 obj.bm = Tuple(*list(map(expand, bm)))
616 obj.bq = Tuple(*list(map(expand, bq)))
617 return obj
619 @property
620 def args(self):
621 return (self.an, self.ap, self.bm, self.bq)
623 def _hashable_content(self):
624 return super()._hashable_content() + self.args
626 def __call__(self, z):
627 return meijerg(self.an, self.ap, self.bm, self.bq, z)
629 def compute_buckets(self):
630 """
631 Compute buckets for the fours sets of parameters.
633 Explanation
634 ===========
636 We guarantee that any two equal Mod objects returned are actually the
637 same, and that the buckets are sorted by real part (an and bq
638 descendending, bm and ap ascending).
640 Examples
641 ========
643 >>> from sympy.simplify.hyperexpand import G_Function
644 >>> from sympy.abc import y
645 >>> from sympy import S
647 >>> a, b = [1, 3, 2, S(3)/2], [1 + y, y, 2, y + 3]
648 >>> G_Function(a, b, [2], [y]).compute_buckets()
649 ({0: [3, 2, 1], 1/2: [3/2]},
650 {0: [2], y: [y, y + 1, y + 3]}, {0: [2]}, {y: [y]})
652 """
653 dicts = pan, pap, pbm, pbq = [defaultdict(list) for i in range(4)]
654 for dic, lis in zip(dicts, (self.an, self.ap, self.bm, self.bq)):
655 for x in lis:
656 dic[_mod1(x)].append(x)
658 for dic, flip in zip(dicts, (True, False, False, True)):
659 for m, items in dic.items():
660 x0 = items[0]
661 items.sort(key=lambda x: x - x0, reverse=flip)
662 dic[m] = items
664 return tuple([dict(w) for w in dicts])
666 @property
667 def signature(self):
668 return (len(self.an), len(self.ap), len(self.bm), len(self.bq))
671# Dummy variable.
672_x = Dummy('x')
674class Formula:
675 """
676 This class represents hypergeometric formulae.
678 Explanation
679 ===========
681 Its data members are:
682 - z, the argument
683 - closed_form, the closed form expression
684 - symbols, the free symbols (parameters) in the formula
685 - func, the function
686 - B, C, M (see _compute_basis)
688 Examples
689 ========
691 >>> from sympy.abc import a, b, z
692 >>> from sympy.simplify.hyperexpand import Formula, Hyper_Function
693 >>> func = Hyper_Function((a/2, a/3 + b, (1+a)/2), (a, b, (a+b)/7))
694 >>> f = Formula(func, z, None, [a, b])
696 """
698 def _compute_basis(self, closed_form):
699 """
700 Compute a set of functions B=(f1, ..., fn), a nxn matrix M
701 and a 1xn matrix C such that:
702 closed_form = C B
703 z d/dz B = M B.
704 """
705 afactors = [_x + a for a in self.func.ap]
706 bfactors = [_x + b - 1 for b in self.func.bq]
707 expr = _x*Mul(*bfactors) - self.z*Mul(*afactors)
708 poly = Poly(expr, _x)
710 n = poly.degree() - 1
711 b = [closed_form]
712 for _ in range(n):
713 b.append(self.z*b[-1].diff(self.z))
715 self.B = Matrix(b)
716 self.C = Matrix([[1] + [0]*n])
718 m = eye(n)
719 m = m.col_insert(0, zeros(n, 1))
720 l = poly.all_coeffs()[1:]
721 l.reverse()
722 self.M = m.row_insert(n, -Matrix([l])/poly.all_coeffs()[0])
724 def __init__(self, func, z, res, symbols, B=None, C=None, M=None):
725 z = sympify(z)
726 res = sympify(res)
727 symbols = [x for x in sympify(symbols) if func.has(x)]
729 self.z = z
730 self.symbols = symbols
731 self.B = B
732 self.C = C
733 self.M = M
734 self.func = func
736 # TODO with symbolic parameters, it could be advantageous
737 # (for prettier answers) to compute a basis only *after*
738 # instantiation
739 if res is not None:
740 self._compute_basis(res)
742 @property
743 def closed_form(self):
744 return reduce(lambda s,m: s+m[0]*m[1], zip(self.C, self.B), S.Zero)
746 def find_instantiations(self, func):
747 """
748 Find substitutions of the free symbols that match ``func``.
750 Return the substitution dictionaries as a list. Note that the returned
751 instantiations need not actually match, or be valid!
753 """
754 from sympy.solvers import solve
755 ap = func.ap
756 bq = func.bq
757 if len(ap) != len(self.func.ap) or len(bq) != len(self.func.bq):
758 raise TypeError('Cannot instantiate other number of parameters')
759 symbol_values = []
760 for a in self.symbols:
761 if a in self.func.ap.args:
762 symbol_values.append(ap)
763 elif a in self.func.bq.args:
764 symbol_values.append(bq)
765 else:
766 raise ValueError("At least one of the parameters of the "
767 "formula must be equal to %s" % (a,))
768 base_repl = [dict(list(zip(self.symbols, values)))
769 for values in product(*symbol_values)]
770 abuckets, bbuckets = [sift(params, _mod1) for params in [ap, bq]]
771 a_inv, b_inv = [{a: len(vals) for a, vals in bucket.items()}
772 for bucket in [abuckets, bbuckets]]
773 critical_values = [[0] for _ in self.symbols]
774 result = []
775 _n = Dummy()
776 for repl in base_repl:
777 symb_a, symb_b = [sift(params, lambda x: _mod1(x.xreplace(repl)))
778 for params in [self.func.ap, self.func.bq]]
779 for bucket, obucket in [(abuckets, symb_a), (bbuckets, symb_b)]:
780 for mod in set(list(bucket.keys()) + list(obucket.keys())):
781 if (mod not in bucket) or (mod not in obucket) \
782 or len(bucket[mod]) != len(obucket[mod]):
783 break
784 for a, vals in zip(self.symbols, critical_values):
785 if repl[a].free_symbols:
786 continue
787 exprs = [expr for expr in obucket[mod] if expr.has(a)]
788 repl0 = repl.copy()
789 repl0[a] += _n
790 for expr in exprs:
791 for target in bucket[mod]:
792 n0, = solve(expr.xreplace(repl0) - target, _n)
793 if n0.free_symbols:
794 raise ValueError("Value should not be true")
795 vals.append(n0)
796 else:
797 values = []
798 for a, vals in zip(self.symbols, critical_values):
799 a0 = repl[a]
800 min_ = floor(min(vals))
801 max_ = ceiling(max(vals))
802 values.append([a0 + n for n in range(min_, max_ + 1)])
803 result.extend(dict(list(zip(self.symbols, l))) for l in product(*values))
804 return result
809class FormulaCollection:
810 """ A collection of formulae to use as origins. """
812 def __init__(self):
813 """ Doing this globally at module init time is a pain ... """
814 self.symbolic_formulae = {}
815 self.concrete_formulae = {}
816 self.formulae = []
818 add_formulae(self.formulae)
820 # Now process the formulae into a helpful form.
821 # These dicts are indexed by (p, q).
823 for f in self.formulae:
824 sizes = f.func.sizes
825 if len(f.symbols) > 0:
826 self.symbolic_formulae.setdefault(sizes, []).append(f)
827 else:
828 inv = f.func.build_invariants()
829 self.concrete_formulae.setdefault(sizes, {})[inv] = f
831 def lookup_origin(self, func):
832 """
833 Given the suitable target ``func``, try to find an origin in our
834 knowledge base.
836 Examples
837 ========
839 >>> from sympy.simplify.hyperexpand import (FormulaCollection,
840 ... Hyper_Function)
841 >>> f = FormulaCollection()
842 >>> f.lookup_origin(Hyper_Function((), ())).closed_form
843 exp(_z)
844 >>> f.lookup_origin(Hyper_Function([1], ())).closed_form
845 HyperRep_power1(-1, _z)
847 >>> from sympy import S
848 >>> i = Hyper_Function([S('1/4'), S('3/4 + 4')], [S.Half])
849 >>> f.lookup_origin(i).closed_form
850 HyperRep_sqrts1(-1/4, _z)
851 """
852 inv = func.build_invariants()
853 sizes = func.sizes
854 if sizes in self.concrete_formulae and \
855 inv in self.concrete_formulae[sizes]:
856 return self.concrete_formulae[sizes][inv]
858 # We don't have a concrete formula. Try to instantiate.
859 if sizes not in self.symbolic_formulae:
860 return None # Too bad...
862 possible = []
863 for f in self.symbolic_formulae[sizes]:
864 repls = f.find_instantiations(func)
865 for repl in repls:
866 func2 = f.func.xreplace(repl)
867 if not func2._is_suitable_origin():
868 continue
869 diff = func2.difficulty(func)
870 if diff == -1:
871 continue
872 possible.append((diff, repl, f, func2))
874 # find the nearest origin
875 possible.sort(key=lambda x: x[0])
876 for _, repl, f, func2 in possible:
877 f2 = Formula(func2, f.z, None, [], f.B.subs(repl),
878 f.C.subs(repl), f.M.subs(repl))
879 if not any(e.has(S.NaN, oo, -oo, zoo) for e in [f2.B, f2.M, f2.C]):
880 return f2
882 return None
885class MeijerFormula:
886 """
887 This class represents a Meijer G-function formula.
889 Its data members are:
890 - z, the argument
891 - symbols, the free symbols (parameters) in the formula
892 - func, the function
893 - B, C, M (c/f ordinary Formula)
894 """
896 def __init__(self, an, ap, bm, bq, z, symbols, B, C, M, matcher):
897 an, ap, bm, bq = [Tuple(*list(map(expand, w))) for w in [an, ap, bm, bq]]
898 self.func = G_Function(an, ap, bm, bq)
899 self.z = z
900 self.symbols = symbols
901 self._matcher = matcher
902 self.B = B
903 self.C = C
904 self.M = M
906 @property
907 def closed_form(self):
908 return reduce(lambda s,m: s+m[0]*m[1], zip(self.C, self.B), S.Zero)
910 def try_instantiate(self, func):
911 """
912 Try to instantiate the current formula to (almost) match func.
913 This uses the _matcher passed on init.
914 """
915 if func.signature != self.func.signature:
916 return None
917 res = self._matcher(func)
918 if res is not None:
919 subs, newfunc = res
920 return MeijerFormula(newfunc.an, newfunc.ap, newfunc.bm, newfunc.bq,
921 self.z, [],
922 self.B.subs(subs), self.C.subs(subs),
923 self.M.subs(subs), None)
926class MeijerFormulaCollection:
927 """
928 This class holds a collection of meijer g formulae.
929 """
931 def __init__(self):
932 formulae = []
933 add_meijerg_formulae(formulae)
934 self.formulae = defaultdict(list)
935 for formula in formulae:
936 self.formulae[formula.func.signature].append(formula)
937 self.formulae = dict(self.formulae)
939 def lookup_origin(self, func):
940 """ Try to find a formula that matches func. """
941 if func.signature not in self.formulae:
942 return None
943 for formula in self.formulae[func.signature]:
944 res = formula.try_instantiate(func)
945 if res is not None:
946 return res
949class Operator:
950 """
951 Base class for operators to be applied to our functions.
953 Explanation
954 ===========
956 These operators are differential operators. They are by convention
957 expressed in the variable D = z*d/dz (although this base class does
958 not actually care).
959 Note that when the operator is applied to an object, we typically do
960 *not* blindly differentiate but instead use a different representation
961 of the z*d/dz operator (see make_derivative_operator).
963 To subclass from this, define a __init__ method that initializes a
964 self._poly variable. This variable stores a polynomial. By convention
965 the generator is z*d/dz, and acts to the right of all coefficients.
967 Thus this poly
968 x**2 + 2*z*x + 1
969 represents the differential operator
970 (z*d/dz)**2 + 2*z**2*d/dz.
972 This class is used only in the implementation of the hypergeometric
973 function expansion algorithm.
974 """
976 def apply(self, obj, op):
977 """
978 Apply ``self`` to the object ``obj``, where the generator is ``op``.
980 Examples
981 ========
983 >>> from sympy.simplify.hyperexpand import Operator
984 >>> from sympy.polys.polytools import Poly
985 >>> from sympy.abc import x, y, z
986 >>> op = Operator()
987 >>> op._poly = Poly(x**2 + z*x + y, x)
988 >>> op.apply(z**7, lambda f: f.diff(z))
989 y*z**7 + 7*z**7 + 42*z**5
990 """
991 coeffs = self._poly.all_coeffs()
992 coeffs.reverse()
993 diffs = [obj]
994 for c in coeffs[1:]:
995 diffs.append(op(diffs[-1]))
996 r = coeffs[0]*diffs[0]
997 for c, d in zip(coeffs[1:], diffs[1:]):
998 r += c*d
999 return r
1002class MultOperator(Operator):
1003 """ Simply multiply by a "constant" """
1005 def __init__(self, p):
1006 self._poly = Poly(p, _x)
1009class ShiftA(Operator):
1010 """ Increment an upper index. """
1012 def __init__(self, ai):
1013 ai = sympify(ai)
1014 if ai == 0:
1015 raise ValueError('Cannot increment zero upper index.')
1016 self._poly = Poly(_x/ai + 1, _x)
1018 def __str__(self):
1019 return '<Increment upper %s.>' % (1/self._poly.all_coeffs()[0])
1022class ShiftB(Operator):
1023 """ Decrement a lower index. """
1025 def __init__(self, bi):
1026 bi = sympify(bi)
1027 if bi == 1:
1028 raise ValueError('Cannot decrement unit lower index.')
1029 self._poly = Poly(_x/(bi - 1) + 1, _x)
1031 def __str__(self):
1032 return '<Decrement lower %s.>' % (1/self._poly.all_coeffs()[0] + 1)
1035class UnShiftA(Operator):
1036 """ Decrement an upper index. """
1038 def __init__(self, ap, bq, i, z):
1039 """ Note: i counts from zero! """
1040 ap, bq, i = list(map(sympify, [ap, bq, i]))
1042 self._ap = ap
1043 self._bq = bq
1044 self._i = i
1046 ap = list(ap)
1047 bq = list(bq)
1048 ai = ap.pop(i) - 1
1050 if ai == 0:
1051 raise ValueError('Cannot decrement unit upper index.')
1053 m = Poly(z*ai, _x)
1054 for a in ap:
1055 m *= Poly(_x + a, _x)
1057 A = Dummy('A')
1058 n = D = Poly(ai*A - ai, A)
1059 for b in bq:
1060 n *= D + (b - 1).as_poly(A)
1062 b0 = -n.nth(0)
1063 if b0 == 0:
1064 raise ValueError('Cannot decrement upper index: '
1065 'cancels with lower')
1067 n = Poly(Poly(n.all_coeffs()[:-1], A).as_expr().subs(A, _x/ai + 1), _x)
1069 self._poly = Poly((n - m)/b0, _x)
1071 def __str__(self):
1072 return '<Decrement upper index #%s of %s, %s.>' % (self._i,
1073 self._ap, self._bq)
1076class UnShiftB(Operator):
1077 """ Increment a lower index. """
1079 def __init__(self, ap, bq, i, z):
1080 """ Note: i counts from zero! """
1081 ap, bq, i = list(map(sympify, [ap, bq, i]))
1083 self._ap = ap
1084 self._bq = bq
1085 self._i = i
1087 ap = list(ap)
1088 bq = list(bq)
1089 bi = bq.pop(i) + 1
1091 if bi == 0:
1092 raise ValueError('Cannot increment -1 lower index.')
1094 m = Poly(_x*(bi - 1), _x)
1095 for b in bq:
1096 m *= Poly(_x + b - 1, _x)
1098 B = Dummy('B')
1099 D = Poly((bi - 1)*B - bi + 1, B)
1100 n = Poly(z, B)
1101 for a in ap:
1102 n *= (D + a.as_poly(B))
1104 b0 = n.nth(0)
1105 if b0 == 0:
1106 raise ValueError('Cannot increment index: cancels with upper')
1108 n = Poly(Poly(n.all_coeffs()[:-1], B).as_expr().subs(
1109 B, _x/(bi - 1) + 1), _x)
1111 self._poly = Poly((m - n)/b0, _x)
1113 def __str__(self):
1114 return '<Increment lower index #%s of %s, %s.>' % (self._i,
1115 self._ap, self._bq)
1118class MeijerShiftA(Operator):
1119 """ Increment an upper b index. """
1121 def __init__(self, bi):
1122 bi = sympify(bi)
1123 self._poly = Poly(bi - _x, _x)
1125 def __str__(self):
1126 return '<Increment upper b=%s.>' % (self._poly.all_coeffs()[1])
1129class MeijerShiftB(Operator):
1130 """ Decrement an upper a index. """
1132 def __init__(self, bi):
1133 bi = sympify(bi)
1134 self._poly = Poly(1 - bi + _x, _x)
1136 def __str__(self):
1137 return '<Decrement upper a=%s.>' % (1 - self._poly.all_coeffs()[1])
1140class MeijerShiftC(Operator):
1141 """ Increment a lower b index. """
1143 def __init__(self, bi):
1144 bi = sympify(bi)
1145 self._poly = Poly(-bi + _x, _x)
1147 def __str__(self):
1148 return '<Increment lower b=%s.>' % (-self._poly.all_coeffs()[1])
1151class MeijerShiftD(Operator):
1152 """ Decrement a lower a index. """
1154 def __init__(self, bi):
1155 bi = sympify(bi)
1156 self._poly = Poly(bi - 1 - _x, _x)
1158 def __str__(self):
1159 return '<Decrement lower a=%s.>' % (self._poly.all_coeffs()[1] + 1)
1162class MeijerUnShiftA(Operator):
1163 """ Decrement an upper b index. """
1165 def __init__(self, an, ap, bm, bq, i, z):
1166 """ Note: i counts from zero! """
1167 an, ap, bm, bq, i = list(map(sympify, [an, ap, bm, bq, i]))
1169 self._an = an
1170 self._ap = ap
1171 self._bm = bm
1172 self._bq = bq
1173 self._i = i
1175 an = list(an)
1176 ap = list(ap)
1177 bm = list(bm)
1178 bq = list(bq)
1179 bi = bm.pop(i) - 1
1181 m = Poly(1, _x) * prod(Poly(b - _x, _x) for b in bm) * prod(Poly(_x - b, _x) for b in bq)
1183 A = Dummy('A')
1184 D = Poly(bi - A, A)
1185 n = Poly(z, A) * prod((D + 1 - a) for a in an) * prod((-D + a - 1) for a in ap)
1187 b0 = n.nth(0)
1188 if b0 == 0:
1189 raise ValueError('Cannot decrement upper b index (cancels)')
1191 n = Poly(Poly(n.all_coeffs()[:-1], A).as_expr().subs(A, bi - _x), _x)
1193 self._poly = Poly((m - n)/b0, _x)
1195 def __str__(self):
1196 return '<Decrement upper b index #%s of %s, %s, %s, %s.>' % (self._i,
1197 self._an, self._ap, self._bm, self._bq)
1200class MeijerUnShiftB(Operator):
1201 """ Increment an upper a index. """
1203 def __init__(self, an, ap, bm, bq, i, z):
1204 """ Note: i counts from zero! """
1205 an, ap, bm, bq, i = list(map(sympify, [an, ap, bm, bq, i]))
1207 self._an = an
1208 self._ap = ap
1209 self._bm = bm
1210 self._bq = bq
1211 self._i = i
1213 an = list(an)
1214 ap = list(ap)
1215 bm = list(bm)
1216 bq = list(bq)
1217 ai = an.pop(i) + 1
1219 m = Poly(z, _x)
1220 for a in an:
1221 m *= Poly(1 - a + _x, _x)
1222 for a in ap:
1223 m *= Poly(a - 1 - _x, _x)
1225 B = Dummy('B')
1226 D = Poly(B + ai - 1, B)
1227 n = Poly(1, B)
1228 for b in bm:
1229 n *= (-D + b)
1230 for b in bq:
1231 n *= (D - b)
1233 b0 = n.nth(0)
1234 if b0 == 0:
1235 raise ValueError('Cannot increment upper a index (cancels)')
1237 n = Poly(Poly(n.all_coeffs()[:-1], B).as_expr().subs(
1238 B, 1 - ai + _x), _x)
1240 self._poly = Poly((m - n)/b0, _x)
1242 def __str__(self):
1243 return '<Increment upper a index #%s of %s, %s, %s, %s.>' % (self._i,
1244 self._an, self._ap, self._bm, self._bq)
1247class MeijerUnShiftC(Operator):
1248 """ Decrement a lower b index. """
1249 # XXX this is "essentially" the same as MeijerUnShiftA. This "essentially"
1250 # can be made rigorous using the functional equation G(1/z) = G'(z),
1251 # where G' denotes a G function of slightly altered parameters.
1252 # However, sorting out the details seems harder than just coding it
1253 # again.
1255 def __init__(self, an, ap, bm, bq, i, z):
1256 """ Note: i counts from zero! """
1257 an, ap, bm, bq, i = list(map(sympify, [an, ap, bm, bq, i]))
1259 self._an = an
1260 self._ap = ap
1261 self._bm = bm
1262 self._bq = bq
1263 self._i = i
1265 an = list(an)
1266 ap = list(ap)
1267 bm = list(bm)
1268 bq = list(bq)
1269 bi = bq.pop(i) - 1
1271 m = Poly(1, _x)
1272 for b in bm:
1273 m *= Poly(b - _x, _x)
1274 for b in bq:
1275 m *= Poly(_x - b, _x)
1277 C = Dummy('C')
1278 D = Poly(bi + C, C)
1279 n = Poly(z, C)
1280 for a in an:
1281 n *= (D + 1 - a)
1282 for a in ap:
1283 n *= (-D + a - 1)
1285 b0 = n.nth(0)
1286 if b0 == 0:
1287 raise ValueError('Cannot decrement lower b index (cancels)')
1289 n = Poly(Poly(n.all_coeffs()[:-1], C).as_expr().subs(C, _x - bi), _x)
1291 self._poly = Poly((m - n)/b0, _x)
1293 def __str__(self):
1294 return '<Decrement lower b index #%s of %s, %s, %s, %s.>' % (self._i,
1295 self._an, self._ap, self._bm, self._bq)
1298class MeijerUnShiftD(Operator):
1299 """ Increment a lower a index. """
1300 # XXX This is essentially the same as MeijerUnShiftA.
1301 # See comment at MeijerUnShiftC.
1303 def __init__(self, an, ap, bm, bq, i, z):
1304 """ Note: i counts from zero! """
1305 an, ap, bm, bq, i = list(map(sympify, [an, ap, bm, bq, i]))
1307 self._an = an
1308 self._ap = ap
1309 self._bm = bm
1310 self._bq = bq
1311 self._i = i
1313 an = list(an)
1314 ap = list(ap)
1315 bm = list(bm)
1316 bq = list(bq)
1317 ai = ap.pop(i) + 1
1319 m = Poly(z, _x)
1320 for a in an:
1321 m *= Poly(1 - a + _x, _x)
1322 for a in ap:
1323 m *= Poly(a - 1 - _x, _x)
1325 B = Dummy('B') # - this is the shift operator `D_I`
1326 D = Poly(ai - 1 - B, B)
1327 n = Poly(1, B)
1328 for b in bm:
1329 n *= (-D + b)
1330 for b in bq:
1331 n *= (D - b)
1333 b0 = n.nth(0)
1334 if b0 == 0:
1335 raise ValueError('Cannot increment lower a index (cancels)')
1337 n = Poly(Poly(n.all_coeffs()[:-1], B).as_expr().subs(
1338 B, ai - 1 - _x), _x)
1340 self._poly = Poly((m - n)/b0, _x)
1342 def __str__(self):
1343 return '<Increment lower a index #%s of %s, %s, %s, %s.>' % (self._i,
1344 self._an, self._ap, self._bm, self._bq)
1347class ReduceOrder(Operator):
1348 """ Reduce Order by cancelling an upper and a lower index. """
1350 def __new__(cls, ai, bj):
1351 """ For convenience if reduction is not possible, return None. """
1352 ai = sympify(ai)
1353 bj = sympify(bj)
1354 n = ai - bj
1355 if not n.is_Integer or n < 0:
1356 return None
1357 if bj.is_integer and bj.is_nonpositive:
1358 return None
1360 expr = Operator.__new__(cls)
1362 p = S.One
1363 for k in range(n):
1364 p *= (_x + bj + k)/(bj + k)
1366 expr._poly = Poly(p, _x)
1367 expr._a = ai
1368 expr._b = bj
1370 return expr
1372 @classmethod
1373 def _meijer(cls, b, a, sign):
1374 """ Cancel b + sign*s and a + sign*s
1375 This is for meijer G functions. """
1376 b = sympify(b)
1377 a = sympify(a)
1378 n = b - a
1379 if n.is_negative or not n.is_Integer:
1380 return None
1382 expr = Operator.__new__(cls)
1384 p = S.One
1385 for k in range(n):
1386 p *= (sign*_x + a + k)
1388 expr._poly = Poly(p, _x)
1389 if sign == -1:
1390 expr._a = b
1391 expr._b = a
1392 else:
1393 expr._b = Add(1, a - 1, evaluate=False)
1394 expr._a = Add(1, b - 1, evaluate=False)
1396 return expr
1398 @classmethod
1399 def meijer_minus(cls, b, a):
1400 return cls._meijer(b, a, -1)
1402 @classmethod
1403 def meijer_plus(cls, a, b):
1404 return cls._meijer(1 - a, 1 - b, 1)
1406 def __str__(self):
1407 return '<Reduce order by cancelling upper %s with lower %s.>' % \
1408 (self._a, self._b)
1411def _reduce_order(ap, bq, gen, key):
1412 """ Order reduction algorithm used in Hypergeometric and Meijer G """
1413 ap = list(ap)
1414 bq = list(bq)
1416 ap.sort(key=key)
1417 bq.sort(key=key)
1419 nap = []
1420 # we will edit bq in place
1421 operators = []
1422 for a in ap:
1423 op = None
1424 for i in range(len(bq)):
1425 op = gen(a, bq[i])
1426 if op is not None:
1427 bq.pop(i)
1428 break
1429 if op is None:
1430 nap.append(a)
1431 else:
1432 operators.append(op)
1434 return nap, bq, operators
1437def reduce_order(func):
1438 """
1439 Given the hypergeometric function ``func``, find a sequence of operators to
1440 reduces order as much as possible.
1442 Explanation
1443 ===========
1445 Return (newfunc, [operators]), where applying the operators to the
1446 hypergeometric function newfunc yields func.
1448 Examples
1449 ========
1451 >>> from sympy.simplify.hyperexpand import reduce_order, Hyper_Function
1452 >>> reduce_order(Hyper_Function((1, 2), (3, 4)))
1453 (Hyper_Function((1, 2), (3, 4)), [])
1454 >>> reduce_order(Hyper_Function((1,), (1,)))
1455 (Hyper_Function((), ()), [<Reduce order by cancelling upper 1 with lower 1.>])
1456 >>> reduce_order(Hyper_Function((2, 4), (3, 3)))
1457 (Hyper_Function((2,), (3,)), [<Reduce order by cancelling
1458 upper 4 with lower 3.>])
1459 """
1460 nap, nbq, operators = _reduce_order(func.ap, func.bq, ReduceOrder, default_sort_key)
1462 return Hyper_Function(Tuple(*nap), Tuple(*nbq)), operators
1465def reduce_order_meijer(func):
1466 """
1467 Given the Meijer G function parameters, ``func``, find a sequence of
1468 operators that reduces order as much as possible.
1470 Return newfunc, [operators].
1472 Examples
1473 ========
1475 >>> from sympy.simplify.hyperexpand import (reduce_order_meijer,
1476 ... G_Function)
1477 >>> reduce_order_meijer(G_Function([3, 4], [5, 6], [3, 4], [1, 2]))[0]
1478 G_Function((4, 3), (5, 6), (3, 4), (2, 1))
1479 >>> reduce_order_meijer(G_Function([3, 4], [5, 6], [3, 4], [1, 8]))[0]
1480 G_Function((3,), (5, 6), (3, 4), (1,))
1481 >>> reduce_order_meijer(G_Function([3, 4], [5, 6], [7, 5], [1, 5]))[0]
1482 G_Function((3,), (), (), (1,))
1483 >>> reduce_order_meijer(G_Function([3, 4], [5, 6], [7, 5], [5, 3]))[0]
1484 G_Function((), (), (), ())
1485 """
1487 nan, nbq, ops1 = _reduce_order(func.an, func.bq, ReduceOrder.meijer_plus,
1488 lambda x: default_sort_key(-x))
1489 nbm, nap, ops2 = _reduce_order(func.bm, func.ap, ReduceOrder.meijer_minus,
1490 default_sort_key)
1492 return G_Function(nan, nap, nbm, nbq), ops1 + ops2
1495def make_derivative_operator(M, z):
1496 """ Create a derivative operator, to be passed to Operator.apply. """
1497 def doit(C):
1498 r = z*C.diff(z) + C*M
1499 r = r.applyfunc(make_simp(z))
1500 return r
1501 return doit
1504def apply_operators(obj, ops, op):
1505 """
1506 Apply the list of operators ``ops`` to object ``obj``, substituting
1507 ``op`` for the generator.
1508 """
1509 res = obj
1510 for o in reversed(ops):
1511 res = o.apply(res, op)
1512 return res
1515def devise_plan(target, origin, z):
1516 """
1517 Devise a plan (consisting of shift and un-shift operators) to be applied
1518 to the hypergeometric function ``target`` to yield ``origin``.
1519 Returns a list of operators.
1521 Examples
1522 ========
1524 >>> from sympy.simplify.hyperexpand import devise_plan, Hyper_Function
1525 >>> from sympy.abc import z
1527 Nothing to do:
1529 >>> devise_plan(Hyper_Function((1, 2), ()), Hyper_Function((1, 2), ()), z)
1530 []
1531 >>> devise_plan(Hyper_Function((), (1, 2)), Hyper_Function((), (1, 2)), z)
1532 []
1534 Very simple plans:
1536 >>> devise_plan(Hyper_Function((2,), ()), Hyper_Function((1,), ()), z)
1537 [<Increment upper 1.>]
1538 >>> devise_plan(Hyper_Function((), (2,)), Hyper_Function((), (1,)), z)
1539 [<Increment lower index #0 of [], [1].>]
1541 Several buckets:
1543 >>> from sympy import S
1544 >>> devise_plan(Hyper_Function((1, S.Half), ()),
1545 ... Hyper_Function((2, S('3/2')), ()), z) #doctest: +NORMALIZE_WHITESPACE
1546 [<Decrement upper index #0 of [3/2, 1], [].>,
1547 <Decrement upper index #0 of [2, 3/2], [].>]
1549 A slightly more complicated plan:
1551 >>> devise_plan(Hyper_Function((1, 3), ()), Hyper_Function((2, 2), ()), z)
1552 [<Increment upper 2.>, <Decrement upper index #0 of [2, 2], [].>]
1554 Another more complicated plan: (note that the ap have to be shifted first!)
1556 >>> devise_plan(Hyper_Function((1, -1), (2,)), Hyper_Function((3, -2), (4,)), z)
1557 [<Decrement lower 3.>, <Decrement lower 4.>,
1558 <Decrement upper index #1 of [-1, 2], [4].>,
1559 <Decrement upper index #1 of [-1, 3], [4].>, <Increment upper -2.>]
1560 """
1561 abuckets, bbuckets, nabuckets, nbbuckets = [sift(params, _mod1) for
1562 params in (target.ap, target.bq, origin.ap, origin.bq)]
1564 if len(list(abuckets.keys())) != len(list(nabuckets.keys())) or \
1565 len(list(bbuckets.keys())) != len(list(nbbuckets.keys())):
1566 raise ValueError('%s not reachable from %s' % (target, origin))
1568 ops = []
1570 def do_shifts(fro, to, inc, dec):
1571 ops = []
1572 for i in range(len(fro)):
1573 if to[i] - fro[i] > 0:
1574 sh = inc
1575 ch = 1
1576 else:
1577 sh = dec
1578 ch = -1
1580 while to[i] != fro[i]:
1581 ops += [sh(fro, i)]
1582 fro[i] += ch
1584 return ops
1586 def do_shifts_a(nal, nbk, al, aother, bother):
1587 """ Shift us from (nal, nbk) to (al, nbk). """
1588 return do_shifts(nal, al, lambda p, i: ShiftA(p[i]),
1589 lambda p, i: UnShiftA(p + aother, nbk + bother, i, z))
1591 def do_shifts_b(nal, nbk, bk, aother, bother):
1592 """ Shift us from (nal, nbk) to (nal, bk). """
1593 return do_shifts(nbk, bk,
1594 lambda p, i: UnShiftB(nal + aother, p + bother, i, z),
1595 lambda p, i: ShiftB(p[i]))
1597 for r in sorted(list(abuckets.keys()) + list(bbuckets.keys()), key=default_sort_key):
1598 al = ()
1599 nal = ()
1600 bk = ()
1601 nbk = ()
1602 if r in abuckets:
1603 al = abuckets[r]
1604 nal = nabuckets[r]
1605 if r in bbuckets:
1606 bk = bbuckets[r]
1607 nbk = nbbuckets[r]
1608 if len(al) != len(nal) or len(bk) != len(nbk):
1609 raise ValueError('%s not reachable from %s' % (target, origin))
1611 al, nal, bk, nbk = [sorted(w, key=default_sort_key)
1612 for w in [al, nal, bk, nbk]]
1614 def others(dic, key):
1615 l = []
1616 for k, value in dic.items():
1617 if k != key:
1618 l += list(dic[k])
1619 return l
1620 aother = others(nabuckets, r)
1621 bother = others(nbbuckets, r)
1623 if len(al) == 0:
1624 # there can be no complications, just shift the bs as we please
1625 ops += do_shifts_b([], nbk, bk, aother, bother)
1626 elif len(bk) == 0:
1627 # there can be no complications, just shift the as as we please
1628 ops += do_shifts_a(nal, [], al, aother, bother)
1629 else:
1630 namax = nal[-1]
1631 amax = al[-1]
1633 if nbk[0] - namax <= 0 or bk[0] - amax <= 0:
1634 raise ValueError('Non-suitable parameters.')
1636 if namax - amax > 0:
1637 # we are going to shift down - first do the as, then the bs
1638 ops += do_shifts_a(nal, nbk, al, aother, bother)
1639 ops += do_shifts_b(al, nbk, bk, aother, bother)
1640 else:
1641 # we are going to shift up - first do the bs, then the as
1642 ops += do_shifts_b(nal, nbk, bk, aother, bother)
1643 ops += do_shifts_a(nal, bk, al, aother, bother)
1645 nabuckets[r] = al
1646 nbbuckets[r] = bk
1648 ops.reverse()
1649 return ops
1652def try_shifted_sum(func, z):
1653 """ Try to recognise a hypergeometric sum that starts from k > 0. """
1654 abuckets, bbuckets = sift(func.ap, _mod1), sift(func.bq, _mod1)
1655 if len(abuckets[S.Zero]) != 1:
1656 return None
1657 r = abuckets[S.Zero][0]
1658 if r <= 0:
1659 return None
1660 if S.Zero not in bbuckets:
1661 return None
1662 l = list(bbuckets[S.Zero])
1663 l.sort()
1664 k = l[0]
1665 if k <= 0:
1666 return None
1668 nap = list(func.ap)
1669 nap.remove(r)
1670 nbq = list(func.bq)
1671 nbq.remove(k)
1672 k -= 1
1673 nap = [x - k for x in nap]
1674 nbq = [x - k for x in nbq]
1676 ops = []
1677 for n in range(r - 1):
1678 ops.append(ShiftA(n + 1))
1679 ops.reverse()
1681 fac = factorial(k)/z**k
1682 fac *= Mul(*[rf(b, k) for b in nbq])
1683 fac /= Mul(*[rf(a, k) for a in nap])
1685 ops += [MultOperator(fac)]
1687 p = 0
1688 for n in range(k):
1689 m = z**n/factorial(n)
1690 m *= Mul(*[rf(a, n) for a in nap])
1691 m /= Mul(*[rf(b, n) for b in nbq])
1692 p += m
1694 return Hyper_Function(nap, nbq), ops, -p
1697def try_polynomial(func, z):
1698 """ Recognise polynomial cases. Returns None if not such a case.
1699 Requires order to be fully reduced. """
1700 abuckets, bbuckets = sift(func.ap, _mod1), sift(func.bq, _mod1)
1701 a0 = abuckets[S.Zero]
1702 b0 = bbuckets[S.Zero]
1703 a0.sort()
1704 b0.sort()
1705 al0 = [x for x in a0 if x <= 0]
1706 bl0 = [x for x in b0 if x <= 0]
1708 if bl0 and all(a < bl0[-1] for a in al0):
1709 return oo
1710 if not al0:
1711 return None
1713 a = al0[-1]
1714 fac = 1
1715 res = S.One
1716 for n in Tuple(*list(range(-a))):
1717 fac *= z
1718 fac /= n + 1
1719 fac *= Mul(*[a + n for a in func.ap])
1720 fac /= Mul(*[b + n for b in func.bq])
1721 res += fac
1722 return res
1725def try_lerchphi(func):
1726 """
1727 Try to find an expression for Hyper_Function ``func`` in terms of Lerch
1728 Transcendents.
1730 Return None if no such expression can be found.
1731 """
1732 # This is actually quite simple, and is described in Roach's paper,
1733 # section 18.
1734 # We don't need to implement the reduction to polylog here, this
1735 # is handled by expand_func.
1737 # First we need to figure out if the summation coefficient is a rational
1738 # function of the summation index, and construct that rational function.
1739 abuckets, bbuckets = sift(func.ap, _mod1), sift(func.bq, _mod1)
1741 paired = {}
1742 for key, value in abuckets.items():
1743 if key != 0 and key not in bbuckets:
1744 return None
1745 bvalue = bbuckets[key]
1746 paired[key] = (list(value), list(bvalue))
1747 bbuckets.pop(key, None)
1748 if bbuckets != {}:
1749 return None
1750 if S.Zero not in abuckets:
1751 return None
1752 aints, bints = paired[S.Zero]
1753 # Account for the additional n! in denominator
1754 paired[S.Zero] = (aints, bints + [1])
1756 t = Dummy('t')
1757 numer = S.One
1758 denom = S.One
1759 for key, (avalue, bvalue) in paired.items():
1760 if len(avalue) != len(bvalue):
1761 return None
1762 # Note that since order has been reduced fully, all the b are
1763 # bigger than all the a they differ from by an integer. In particular
1764 # if there are any negative b left, this function is not well-defined.
1765 for a, b in zip(avalue, bvalue):
1766 if (a - b).is_positive:
1767 k = a - b
1768 numer *= rf(b + t, k)
1769 denom *= rf(b, k)
1770 else:
1771 k = b - a
1772 numer *= rf(a, k)
1773 denom *= rf(a + t, k)
1775 # Now do a partial fraction decomposition.
1776 # We assemble two structures: a list monomials of pairs (a, b) representing
1777 # a*t**b (b a non-negative integer), and a dict terms, where
1778 # terms[a] = [(b, c)] means that there is a term b/(t-a)**c.
1779 part = apart(numer/denom, t)
1780 args = Add.make_args(part)
1781 monomials = []
1782 terms = {}
1783 for arg in args:
1784 numer, denom = arg.as_numer_denom()
1785 if not denom.has(t):
1786 p = Poly(numer, t)
1787 if not p.is_monomial:
1788 raise TypeError("p should be monomial")
1789 ((b, ), a) = p.LT()
1790 monomials += [(a/denom, b)]
1791 continue
1792 if numer.has(t):
1793 raise NotImplementedError('Need partial fraction decomposition'
1794 ' with linear denominators')
1795 indep, [dep] = denom.as_coeff_mul(t)
1796 n = 1
1797 if dep.is_Pow:
1798 n = dep.exp
1799 dep = dep.base
1800 if dep == t:
1801 a == 0
1802 elif dep.is_Add:
1803 a, tmp = dep.as_independent(t)
1804 b = 1
1805 if tmp != t:
1806 b, _ = tmp.as_independent(t)
1807 if dep != b*t + a:
1808 raise NotImplementedError('unrecognised form %s' % dep)
1809 a /= b
1810 indep *= b**n
1811 else:
1812 raise NotImplementedError('unrecognised form of partial fraction')
1813 terms.setdefault(a, []).append((numer/indep, n))
1815 # Now that we have this information, assemble our formula. All the
1816 # monomials yield rational functions and go into one basis element.
1817 # The terms[a] are related by differentiation. If the largest exponent is
1818 # n, we need lerchphi(z, k, a) for k = 1, 2, ..., n.
1819 # deriv maps a basis to its derivative, expressed as a C(z)-linear
1820 # combination of other basis elements.
1821 deriv = {}
1822 coeffs = {}
1823 z = Dummy('z')
1824 monomials.sort(key=lambda x: x[1])
1825 mon = {0: 1/(1 - z)}
1826 if monomials:
1827 for k in range(monomials[-1][1]):
1828 mon[k + 1] = z*mon[k].diff(z)
1829 for a, n in monomials:
1830 coeffs.setdefault(S.One, []).append(a*mon[n])
1831 for a, l in terms.items():
1832 for c, k in l:
1833 coeffs.setdefault(lerchphi(z, k, a), []).append(c)
1834 l.sort(key=lambda x: x[1])
1835 for k in range(2, l[-1][1] + 1):
1836 deriv[lerchphi(z, k, a)] = [(-a, lerchphi(z, k, a)),
1837 (1, lerchphi(z, k - 1, a))]
1838 deriv[lerchphi(z, 1, a)] = [(-a, lerchphi(z, 1, a)),
1839 (1/(1 - z), S.One)]
1840 trans = {}
1841 for n, b in enumerate([S.One] + list(deriv.keys())):
1842 trans[b] = n
1843 basis = [expand_func(b) for (b, _) in sorted(trans.items(),
1844 key=lambda x:x[1])]
1845 B = Matrix(basis)
1846 C = Matrix([[0]*len(B)])
1847 for b, c in coeffs.items():
1848 C[trans[b]] = Add(*c)
1849 M = zeros(len(B))
1850 for b, l in deriv.items():
1851 for c, b2 in l:
1852 M[trans[b], trans[b2]] = c
1853 return Formula(func, z, None, [], B, C, M)
1856def build_hypergeometric_formula(func):
1857 """
1858 Create a formula object representing the hypergeometric function ``func``.
1860 """
1861 # We know that no `ap` are negative integers, otherwise "detect poly"
1862 # would have kicked in. However, `ap` could be empty. In this case we can
1863 # use a different basis.
1864 # I'm not aware of a basis that works in all cases.
1865 z = Dummy('z')
1866 if func.ap:
1867 afactors = [_x + a for a in func.ap]
1868 bfactors = [_x + b - 1 for b in func.bq]
1869 expr = _x*Mul(*bfactors) - z*Mul(*afactors)
1870 poly = Poly(expr, _x)
1871 n = poly.degree()
1872 basis = []
1873 M = zeros(n)
1874 for k in range(n):
1875 a = func.ap[0] + k
1876 basis += [hyper([a] + list(func.ap[1:]), func.bq, z)]
1877 if k < n - 1:
1878 M[k, k] = -a
1879 M[k, k + 1] = a
1880 B = Matrix(basis)
1881 C = Matrix([[1] + [0]*(n - 1)])
1882 derivs = [eye(n)]
1883 for k in range(n):
1884 derivs.append(M*derivs[k])
1885 l = poly.all_coeffs()
1886 l.reverse()
1887 res = [0]*n
1888 for k, c in enumerate(l):
1889 for r, d in enumerate(C*derivs[k]):
1890 res[r] += c*d
1891 for k, c in enumerate(res):
1892 M[n - 1, k] = -c/derivs[n - 1][0, n - 1]/poly.all_coeffs()[0]
1893 return Formula(func, z, None, [], B, C, M)
1894 else:
1895 # Since there are no `ap`, none of the `bq` can be non-positive
1896 # integers.
1897 basis = []
1898 bq = list(func.bq[:])
1899 for i in range(len(bq)):
1900 basis += [hyper([], bq, z)]
1901 bq[i] += 1
1902 basis += [hyper([], bq, z)]
1903 B = Matrix(basis)
1904 n = len(B)
1905 C = Matrix([[1] + [0]*(n - 1)])
1906 M = zeros(n)
1907 M[0, n - 1] = z/Mul(*func.bq)
1908 for k in range(1, n):
1909 M[k, k - 1] = func.bq[k - 1]
1910 M[k, k] = -func.bq[k - 1]
1911 return Formula(func, z, None, [], B, C, M)
1914def hyperexpand_special(ap, bq, z):
1915 """
1916 Try to find a closed-form expression for hyper(ap, bq, z), where ``z``
1917 is supposed to be a "special" value, e.g. 1.
1919 This function tries various of the classical summation formulae
1920 (Gauss, Saalschuetz, etc).
1921 """
1922 # This code is very ad-hoc. There are many clever algorithms
1923 # (notably Zeilberger's) related to this problem.
1924 # For now we just want a few simple cases to work.
1925 p, q = len(ap), len(bq)
1926 z_ = z
1927 z = unpolarify(z)
1928 if z == 0:
1929 return S.One
1930 from sympy.simplify.simplify import simplify
1931 if p == 2 and q == 1:
1932 # 2F1
1933 a, b, c = ap + bq
1934 if z == 1:
1935 # Gauss
1936 return gamma(c - a - b)*gamma(c)/gamma(c - a)/gamma(c - b)
1937 if z == -1 and simplify(b - a + c) == 1:
1938 b, a = a, b
1939 if z == -1 and simplify(a - b + c) == 1:
1940 # Kummer
1941 if b.is_integer and b.is_negative:
1942 return 2*cos(pi*b/2)*gamma(-b)*gamma(b - a + 1) \
1943 /gamma(-b/2)/gamma(b/2 - a + 1)
1944 else:
1945 return gamma(b/2 + 1)*gamma(b - a + 1) \
1946 /gamma(b + 1)/gamma(b/2 - a + 1)
1947 # TODO tons of more formulae
1948 # investigate what algorithms exist
1949 return hyper(ap, bq, z_)
1951_collection = None
1954def _hyperexpand(func, z, ops0=[], z0=Dummy('z0'), premult=1, prem=0,
1955 rewrite='default'):
1956 """
1957 Try to find an expression for the hypergeometric function ``func``.
1959 Explanation
1960 ===========
1962 The result is expressed in terms of a dummy variable ``z0``. Then it
1963 is multiplied by ``premult``. Then ``ops0`` is applied.
1964 ``premult`` must be a*z**prem for some a independent of ``z``.
1965 """
1967 if z.is_zero:
1968 return S.One
1970 from sympy.simplify.simplify import simplify
1972 z = polarify(z, subs=False)
1973 if rewrite == 'default':
1974 rewrite = 'nonrepsmall'
1976 def carryout_plan(f, ops):
1977 C = apply_operators(f.C.subs(f.z, z0), ops,
1978 make_derivative_operator(f.M.subs(f.z, z0), z0))
1979 C = apply_operators(C, ops0,
1980 make_derivative_operator(f.M.subs(f.z, z0)
1981 + prem*eye(f.M.shape[0]), z0))
1983 if premult == 1:
1984 C = C.applyfunc(make_simp(z0))
1985 r = reduce(lambda s,m: s+m[0]*m[1], zip(C, f.B.subs(f.z, z0)), S.Zero)*premult
1986 res = r.subs(z0, z)
1987 if rewrite:
1988 res = res.rewrite(rewrite)
1989 return res
1991 # TODO
1992 # The following would be possible:
1993 # *) PFD Duplication (see Kelly Roach's paper)
1994 # *) In a similar spirit, try_lerchphi() can be generalised considerably.
1996 global _collection
1997 if _collection is None:
1998 _collection = FormulaCollection()
2000 debug('Trying to expand hypergeometric function ', func)
2002 # First reduce order as much as possible.
2003 func, ops = reduce_order(func)
2004 if ops:
2005 debug(' Reduced order to ', func)
2006 else:
2007 debug(' Could not reduce order.')
2009 # Now try polynomial cases
2010 res = try_polynomial(func, z0)
2011 if res is not None:
2012 debug(' Recognised polynomial.')
2013 p = apply_operators(res, ops, lambda f: z0*f.diff(z0))
2014 p = apply_operators(p*premult, ops0, lambda f: z0*f.diff(z0))
2015 return unpolarify(simplify(p).subs(z0, z))
2017 # Try to recognise a shifted sum.
2018 p = S.Zero
2019 res = try_shifted_sum(func, z0)
2020 if res is not None:
2021 func, nops, p = res
2022 debug(' Recognised shifted sum, reduced order to ', func)
2023 ops += nops
2025 # apply the plan for poly
2026 p = apply_operators(p, ops, lambda f: z0*f.diff(z0))
2027 p = apply_operators(p*premult, ops0, lambda f: z0*f.diff(z0))
2028 p = simplify(p).subs(z0, z)
2030 # Try special expansions early.
2031 if unpolarify(z) in [1, -1] and (len(func.ap), len(func.bq)) == (2, 1):
2032 f = build_hypergeometric_formula(func)
2033 r = carryout_plan(f, ops).replace(hyper, hyperexpand_special)
2034 if not r.has(hyper):
2035 return r + p
2037 # Try to find a formula in our collection
2038 formula = _collection.lookup_origin(func)
2040 # Now try a lerch phi formula
2041 if formula is None:
2042 formula = try_lerchphi(func)
2044 if formula is None:
2045 debug(' Could not find an origin. ',
2046 'Will return answer in terms of '
2047 'simpler hypergeometric functions.')
2048 formula = build_hypergeometric_formula(func)
2050 debug(' Found an origin: ', formula.closed_form, ' ', formula.func)
2052 # We need to find the operators that convert formula into func.
2053 ops += devise_plan(func, formula.func, z0)
2055 # Now carry out the plan.
2056 r = carryout_plan(formula, ops) + p
2058 return powdenest(r, polar=True).replace(hyper, hyperexpand_special)
2061def devise_plan_meijer(fro, to, z):
2062 """
2063 Find operators to convert G-function ``fro`` into G-function ``to``.
2065 Explanation
2066 ===========
2068 It is assumed that ``fro`` and ``to`` have the same signatures, and that in fact
2069 any corresponding pair of parameters differs by integers, and a direct path
2070 is possible. I.e. if there are parameters a1 b1 c1 and a2 b2 c2 it is
2071 assumed that a1 can be shifted to a2, etc. The only thing this routine
2072 determines is the order of shifts to apply, nothing clever will be tried.
2073 It is also assumed that ``fro`` is suitable.
2075 Examples
2076 ========
2078 >>> from sympy.simplify.hyperexpand import (devise_plan_meijer,
2079 ... G_Function)
2080 >>> from sympy.abc import z
2082 Empty plan:
2084 >>> devise_plan_meijer(G_Function([1], [2], [3], [4]),
2085 ... G_Function([1], [2], [3], [4]), z)
2086 []
2088 Very simple plans:
2090 >>> devise_plan_meijer(G_Function([0], [], [], []),
2091 ... G_Function([1], [], [], []), z)
2092 [<Increment upper a index #0 of [0], [], [], [].>]
2093 >>> devise_plan_meijer(G_Function([0], [], [], []),
2094 ... G_Function([-1], [], [], []), z)
2095 [<Decrement upper a=0.>]
2096 >>> devise_plan_meijer(G_Function([], [1], [], []),
2097 ... G_Function([], [2], [], []), z)
2098 [<Increment lower a index #0 of [], [1], [], [].>]
2100 Slightly more complicated plans:
2102 >>> devise_plan_meijer(G_Function([0], [], [], []),
2103 ... G_Function([2], [], [], []), z)
2104 [<Increment upper a index #0 of [1], [], [], [].>,
2105 <Increment upper a index #0 of [0], [], [], [].>]
2106 >>> devise_plan_meijer(G_Function([0], [], [0], []),
2107 ... G_Function([-1], [], [1], []), z)
2108 [<Increment upper b=0.>, <Decrement upper a=0.>]
2110 Order matters:
2112 >>> devise_plan_meijer(G_Function([0], [], [0], []),
2113 ... G_Function([1], [], [1], []), z)
2114 [<Increment upper a index #0 of [0], [], [1], [].>, <Increment upper b=0.>]
2115 """
2116 # TODO for now, we use the following simple heuristic: inverse-shift
2117 # when possible, shift otherwise. Give up if we cannot make progress.
2119 def try_shift(f, t, shifter, diff, counter):
2120 """ Try to apply ``shifter`` in order to bring some element in ``f``
2121 nearer to its counterpart in ``to``. ``diff`` is +/- 1 and
2122 determines the effect of ``shifter``. Counter is a list of elements
2123 blocking the shift.
2125 Return an operator if change was possible, else None.
2126 """
2127 for idx, (a, b) in enumerate(zip(f, t)):
2128 if (
2129 (a - b).is_integer and (b - a)/diff > 0 and
2130 all(a != x for x in counter)):
2131 sh = shifter(idx)
2132 f[idx] += diff
2133 return sh
2134 fan = list(fro.an)
2135 fap = list(fro.ap)
2136 fbm = list(fro.bm)
2137 fbq = list(fro.bq)
2138 ops = []
2139 change = True
2140 while change:
2141 change = False
2142 op = try_shift(fan, to.an,
2143 lambda i: MeijerUnShiftB(fan, fap, fbm, fbq, i, z),
2144 1, fbm + fbq)
2145 if op is not None:
2146 ops += [op]
2147 change = True
2148 continue
2149 op = try_shift(fap, to.ap,
2150 lambda i: MeijerUnShiftD(fan, fap, fbm, fbq, i, z),
2151 1, fbm + fbq)
2152 if op is not None:
2153 ops += [op]
2154 change = True
2155 continue
2156 op = try_shift(fbm, to.bm,
2157 lambda i: MeijerUnShiftA(fan, fap, fbm, fbq, i, z),
2158 -1, fan + fap)
2159 if op is not None:
2160 ops += [op]
2161 change = True
2162 continue
2163 op = try_shift(fbq, to.bq,
2164 lambda i: MeijerUnShiftC(fan, fap, fbm, fbq, i, z),
2165 -1, fan + fap)
2166 if op is not None:
2167 ops += [op]
2168 change = True
2169 continue
2170 op = try_shift(fan, to.an, lambda i: MeijerShiftB(fan[i]), -1, [])
2171 if op is not None:
2172 ops += [op]
2173 change = True
2174 continue
2175 op = try_shift(fap, to.ap, lambda i: MeijerShiftD(fap[i]), -1, [])
2176 if op is not None:
2177 ops += [op]
2178 change = True
2179 continue
2180 op = try_shift(fbm, to.bm, lambda i: MeijerShiftA(fbm[i]), 1, [])
2181 if op is not None:
2182 ops += [op]
2183 change = True
2184 continue
2185 op = try_shift(fbq, to.bq, lambda i: MeijerShiftC(fbq[i]), 1, [])
2186 if op is not None:
2187 ops += [op]
2188 change = True
2189 continue
2190 if fan != list(to.an) or fap != list(to.ap) or fbm != list(to.bm) or \
2191 fbq != list(to.bq):
2192 raise NotImplementedError('Could not devise plan.')
2193 ops.reverse()
2194 return ops
2196_meijercollection = None
2199def _meijergexpand(func, z0, allow_hyper=False, rewrite='default',
2200 place=None):
2201 """
2202 Try to find an expression for the Meijer G function specified
2203 by the G_Function ``func``. If ``allow_hyper`` is True, then returning
2204 an expression in terms of hypergeometric functions is allowed.
2206 Currently this just does Slater's theorem.
2207 If expansions exist both at zero and at infinity, ``place``
2208 can be set to ``0`` or ``zoo`` for the preferred choice.
2209 """
2210 global _meijercollection
2211 if _meijercollection is None:
2212 _meijercollection = MeijerFormulaCollection()
2213 if rewrite == 'default':
2214 rewrite = None
2216 func0 = func
2217 debug('Try to expand Meijer G function corresponding to ', func)
2219 # We will play games with analytic continuation - rather use a fresh symbol
2220 z = Dummy('z')
2222 func, ops = reduce_order_meijer(func)
2223 if ops:
2224 debug(' Reduced order to ', func)
2225 else:
2226 debug(' Could not reduce order.')
2228 # Try to find a direct formula
2229 f = _meijercollection.lookup_origin(func)
2230 if f is not None:
2231 debug(' Found a Meijer G formula: ', f.func)
2232 ops += devise_plan_meijer(f.func, func, z)
2234 # Now carry out the plan.
2235 C = apply_operators(f.C.subs(f.z, z), ops,
2236 make_derivative_operator(f.M.subs(f.z, z), z))
2238 C = C.applyfunc(make_simp(z))
2239 r = C*f.B.subs(f.z, z)
2240 r = r[0].subs(z, z0)
2241 return powdenest(r, polar=True)
2243 debug(" Could not find a direct formula. Trying Slater's theorem.")
2245 # TODO the following would be possible:
2246 # *) Paired Index Theorems
2247 # *) PFD Duplication
2248 # (See Kelly Roach's paper for details on either.)
2249 #
2250 # TODO Also, we tend to create combinations of gamma functions that can be
2251 # simplified.
2253 def can_do(pbm, pap):
2254 """ Test if slater applies. """
2255 for i in pbm:
2256 if len(pbm[i]) > 1:
2257 l = 0
2258 if i in pap:
2259 l = len(pap[i])
2260 if l + 1 < len(pbm[i]):
2261 return False
2262 return True
2264 def do_slater(an, bm, ap, bq, z, zfinal):
2265 # zfinal is the value that will eventually be substituted for z.
2266 # We pass it to _hyperexpand to improve performance.
2267 func = G_Function(an, bm, ap, bq)
2268 _, pbm, pap, _ = func.compute_buckets()
2269 if not can_do(pbm, pap):
2270 return S.Zero, False
2272 cond = len(an) + len(ap) < len(bm) + len(bq)
2273 if len(an) + len(ap) == len(bm) + len(bq):
2274 cond = abs(z) < 1
2275 if cond is False:
2276 return S.Zero, False
2278 res = S.Zero
2279 for m in pbm:
2280 if len(pbm[m]) == 1:
2281 bh = pbm[m][0]
2282 fac = 1
2283 bo = list(bm)
2284 bo.remove(bh)
2285 for bj in bo:
2286 fac *= gamma(bj - bh)
2287 for aj in an:
2288 fac *= gamma(1 + bh - aj)
2289 for bj in bq:
2290 fac /= gamma(1 + bh - bj)
2291 for aj in ap:
2292 fac /= gamma(aj - bh)
2293 nap = [1 + bh - a for a in list(an) + list(ap)]
2294 nbq = [1 + bh - b for b in list(bo) + list(bq)]
2296 k = polar_lift(S.NegativeOne**(len(ap) - len(bm)))
2297 harg = k*zfinal
2298 # NOTE even though k "is" +-1, this has to be t/k instead of
2299 # t*k ... we are using polar numbers for consistency!
2300 premult = (t/k)**bh
2301 hyp = _hyperexpand(Hyper_Function(nap, nbq), harg, ops,
2302 t, premult, bh, rewrite=None)
2303 res += fac * hyp
2304 else:
2305 b_ = pbm[m][0]
2306 ki = [bi - b_ for bi in pbm[m][1:]]
2307 u = len(ki)
2308 li = [ai - b_ for ai in pap[m][:u + 1]]
2309 bo = list(bm)
2310 for b in pbm[m]:
2311 bo.remove(b)
2312 ao = list(ap)
2313 for a in pap[m][:u]:
2314 ao.remove(a)
2315 lu = li[-1]
2316 di = [l - k for (l, k) in zip(li, ki)]
2318 # We first work out the integrand:
2319 s = Dummy('s')
2320 integrand = z**s
2321 for b in bm:
2322 if not Mod(b, 1) and b.is_Number:
2323 b = int(round(b))
2324 integrand *= gamma(b - s)
2325 for a in an:
2326 integrand *= gamma(1 - a + s)
2327 for b in bq:
2328 integrand /= gamma(1 - b + s)
2329 for a in ap:
2330 integrand /= gamma(a - s)
2332 # Now sum the finitely many residues:
2333 # XXX This speeds up some cases - is it a good idea?
2334 integrand = expand_func(integrand)
2335 for r in range(int(round(lu))):
2336 resid = residue(integrand, s, b_ + r)
2337 resid = apply_operators(resid, ops, lambda f: z*f.diff(z))
2338 res -= resid
2340 # Now the hypergeometric term.
2341 au = b_ + lu
2342 k = polar_lift(S.NegativeOne**(len(ao) + len(bo) + 1))
2343 harg = k*zfinal
2344 premult = (t/k)**au
2345 nap = [1 + au - a for a in list(an) + list(ap)] + [1]
2346 nbq = [1 + au - b for b in list(bm) + list(bq)]
2348 hyp = _hyperexpand(Hyper_Function(nap, nbq), harg, ops,
2349 t, premult, au, rewrite=None)
2351 C = S.NegativeOne**(lu)/factorial(lu)
2352 for i in range(u):
2353 C *= S.NegativeOne**di[i]/rf(lu - li[i] + 1, di[i])
2354 for a in an:
2355 C *= gamma(1 - a + au)
2356 for b in bo:
2357 C *= gamma(b - au)
2358 for a in ao:
2359 C /= gamma(a - au)
2360 for b in bq:
2361 C /= gamma(1 - b + au)
2363 res += C*hyp
2365 return res, cond
2367 t = Dummy('t')
2368 slater1, cond1 = do_slater(func.an, func.bm, func.ap, func.bq, z, z0)
2370 def tr(l):
2371 return [1 - x for x in l]
2373 for op in ops:
2374 op._poly = Poly(op._poly.subs({z: 1/t, _x: -_x}), _x)
2375 slater2, cond2 = do_slater(tr(func.bm), tr(func.an), tr(func.bq), tr(func.ap),
2376 t, 1/z0)
2378 slater1 = powdenest(slater1.subs(z, z0), polar=True)
2379 slater2 = powdenest(slater2.subs(t, 1/z0), polar=True)
2380 if not isinstance(cond2, bool):
2381 cond2 = cond2.subs(t, 1/z)
2383 m = func(z)
2384 if m.delta > 0 or \
2385 (m.delta == 0 and len(m.ap) == len(m.bq) and
2386 (re(m.nu) < -1) is not False and polar_lift(z0) == polar_lift(1)):
2387 # The condition delta > 0 means that the convergence region is
2388 # connected. Any expression we find can be continued analytically
2389 # to the entire convergence region.
2390 # The conditions delta==0, p==q, re(nu) < -1 imply that G is continuous
2391 # on the positive reals, so the values at z=1 agree.
2392 if cond1 is not False:
2393 cond1 = True
2394 if cond2 is not False:
2395 cond2 = True
2397 if cond1 is True:
2398 slater1 = slater1.rewrite(rewrite or 'nonrep')
2399 else:
2400 slater1 = slater1.rewrite(rewrite or 'nonrepsmall')
2401 if cond2 is True:
2402 slater2 = slater2.rewrite(rewrite or 'nonrep')
2403 else:
2404 slater2 = slater2.rewrite(rewrite or 'nonrepsmall')
2406 if cond1 is not False and cond2 is not False:
2407 # If one condition is False, there is no choice.
2408 if place == 0:
2409 cond2 = False
2410 if place == zoo:
2411 cond1 = False
2413 if not isinstance(cond1, bool):
2414 cond1 = cond1.subs(z, z0)
2415 if not isinstance(cond2, bool):
2416 cond2 = cond2.subs(z, z0)
2418 def weight(expr, cond):
2419 if cond is True:
2420 c0 = 0
2421 elif cond is False:
2422 c0 = 1
2423 else:
2424 c0 = 2
2425 if expr.has(oo, zoo, -oo, nan):
2426 # XXX this actually should not happen, but consider
2427 # S('meijerg(((0, -1/2, 0, -1/2, 1/2), ()), ((0,),
2428 # (-1/2, -1/2, -1/2, -1)), exp_polar(I*pi))/4')
2429 c0 = 3
2430 return (c0, expr.count(hyper), expr.count_ops())
2432 w1 = weight(slater1, cond1)
2433 w2 = weight(slater2, cond2)
2434 if min(w1, w2) <= (0, 1, oo):
2435 if w1 < w2:
2436 return slater1
2437 else:
2438 return slater2
2439 if max(w1[0], w2[0]) <= 1 and max(w1[1], w2[1]) <= 1:
2440 return Piecewise((slater1, cond1), (slater2, cond2), (func0(z0), True))
2442 # We couldn't find an expression without hypergeometric functions.
2443 # TODO it would be helpful to give conditions under which the integral
2444 # is known to diverge.
2445 r = Piecewise((slater1, cond1), (slater2, cond2), (func0(z0), True))
2446 if r.has(hyper) and not allow_hyper:
2447 debug(' Could express using hypergeometric functions, '
2448 'but not allowed.')
2449 if not r.has(hyper) or allow_hyper:
2450 return r
2452 return func0(z0)
2455def hyperexpand(f, allow_hyper=False, rewrite='default', place=None):
2456 """
2457 Expand hypergeometric functions. If allow_hyper is True, allow partial
2458 simplification (that is a result different from input,
2459 but still containing hypergeometric functions).
2461 If a G-function has expansions both at zero and at infinity,
2462 ``place`` can be set to ``0`` or ``zoo`` to indicate the
2463 preferred choice.
2465 Examples
2466 ========
2468 >>> from sympy.simplify.hyperexpand import hyperexpand
2469 >>> from sympy.functions import hyper
2470 >>> from sympy.abc import z
2471 >>> hyperexpand(hyper([], [], z))
2472 exp(z)
2474 Non-hyperegeometric parts of the expression and hypergeometric expressions
2475 that are not recognised are left unchanged:
2477 >>> hyperexpand(1 + hyper([1, 1, 1], [], z))
2478 hyper((1, 1, 1), (), z) + 1
2479 """
2480 f = sympify(f)
2482 def do_replace(ap, bq, z):
2483 r = _hyperexpand(Hyper_Function(ap, bq), z, rewrite=rewrite)
2484 if r is None:
2485 return hyper(ap, bq, z)
2486 else:
2487 return r
2489 def do_meijer(ap, bq, z):
2490 r = _meijergexpand(G_Function(ap[0], ap[1], bq[0], bq[1]), z,
2491 allow_hyper, rewrite=rewrite, place=place)
2492 if not r.has(nan, zoo, oo, -oo):
2493 return r
2494 return f.replace(hyper, do_replace).replace(meijerg, do_meijer)