Coverage for /usr/lib/python3/dist-packages/sympy/integrals/meijerint.py: 8%
1134 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2Integrate functions by rewriting them as Meijer G-functions.
4There are three user-visible functions that can be used by other parts of the
5sympy library to solve various integration problems:
7- meijerint_indefinite
8- meijerint_definite
9- meijerint_inversion
11They can be used to compute, respectively, indefinite integrals, definite
12integrals over intervals of the real line, and inverse laplace-type integrals
13(from c-I*oo to c+I*oo). See the respective docstrings for details.
15The main references for this are:
17[L] Luke, Y. L. (1969), The Special Functions and Their Approximations,
18 Volume 1
20[R] Kelly B. Roach. Meijer G Function Representations.
21 In: Proceedings of the 1997 International Symposium on Symbolic and
22 Algebraic Computation, pages 205-211, New York, 1997. ACM.
24[P] A. P. Prudnikov, Yu. A. Brychkov and O. I. Marichev (1990).
25 Integrals and Series: More Special Functions, Vol. 3,.
26 Gordon and Breach Science Publisher
27"""
29from __future__ import annotations
30import itertools
32from sympy import SYMPY_DEBUG
33from sympy.core import S, Expr
34from sympy.core.add import Add
35from sympy.core.basic import Basic
36from sympy.core.cache import cacheit
37from sympy.core.containers import Tuple
38from sympy.core.exprtools import factor_terms
39from sympy.core.function import (expand, expand_mul, expand_power_base,
40 expand_trig, Function)
41from sympy.core.mul import Mul
42from sympy.core.numbers import ilcm, Rational, pi
43from sympy.core.relational import Eq, Ne, _canonical_coeff
44from sympy.core.sorting import default_sort_key, ordered
45from sympy.core.symbol import Dummy, symbols, Wild, Symbol
46from sympy.core.sympify import sympify
47from sympy.functions.combinatorial.factorials import factorial
48from sympy.functions.elementary.complexes import (re, im, arg, Abs, sign,
49 unpolarify, polarify, polar_lift, principal_branch, unbranched_argument,
50 periodic_argument)
51from sympy.functions.elementary.exponential import exp, exp_polar, log
52from sympy.functions.elementary.integers import ceiling
53from sympy.functions.elementary.hyperbolic import (cosh, sinh,
54 _rewrite_hyperbolics_as_exp, HyperbolicFunction)
55from sympy.functions.elementary.miscellaneous import sqrt
56from sympy.functions.elementary.piecewise import Piecewise, piecewise_fold
57from sympy.functions.elementary.trigonometric import (cos, sin, sinc,
58 TrigonometricFunction)
59from sympy.functions.special.bessel import besselj, bessely, besseli, besselk
60from sympy.functions.special.delta_functions import DiracDelta, Heaviside
61from sympy.functions.special.elliptic_integrals import elliptic_k, elliptic_e
62from sympy.functions.special.error_functions import (erf, erfc, erfi, Ei,
63 expint, Si, Ci, Shi, Chi, fresnels, fresnelc)
64from sympy.functions.special.gamma_functions import gamma
65from sympy.functions.special.hyper import hyper, meijerg
66from sympy.functions.special.singularity_functions import SingularityFunction
67from .integrals import Integral
68from sympy.logic.boolalg import And, Or, BooleanAtom, Not, BooleanFunction
69from sympy.polys import cancel, factor
70from sympy.utilities.iterables import multiset_partitions
71from sympy.utilities.misc import debug as _debug
72from sympy.utilities.misc import debugf as _debugf
74# keep this at top for easy reference
75z = Dummy('z')
78def _has(res, *f):
79 # return True if res has f; in the case of Piecewise
80 # only return True if *all* pieces have f
81 res = piecewise_fold(res)
82 if getattr(res, 'is_Piecewise', False):
83 return all(_has(i, *f) for i in res.args)
84 return res.has(*f)
87def _create_lookup_table(table):
88 """ Add formulae for the function -> meijerg lookup table. """
89 def wild(n):
90 return Wild(n, exclude=[z])
91 p, q, a, b, c = list(map(wild, 'pqabc'))
92 n = Wild('n', properties=[lambda x: x.is_Integer and x > 0])
93 t = p*z**q
95 def add(formula, an, ap, bm, bq, arg=t, fac=S.One, cond=True, hint=True):
96 table.setdefault(_mytype(formula, z), []).append((formula,
97 [(fac, meijerg(an, ap, bm, bq, arg))], cond, hint))
99 def addi(formula, inst, cond, hint=True):
100 table.setdefault(
101 _mytype(formula, z), []).append((formula, inst, cond, hint))
103 def constant(a):
104 return [(a, meijerg([1], [], [], [0], z)),
105 (a, meijerg([], [1], [0], [], z))]
106 table[()] = [(a, constant(a), True, True)]
108 # [P], Section 8.
109 class IsNonPositiveInteger(Function):
111 @classmethod
112 def eval(cls, arg):
113 arg = unpolarify(arg)
114 if arg.is_Integer is True:
115 return arg <= 0
117 # Section 8.4.2
118 # TODO this needs more polar_lift (c/f entry for exp)
119 add(Heaviside(t - b)*(t - b)**(a - 1), [a], [], [], [0], t/b,
120 gamma(a)*b**(a - 1), And(b > 0))
121 add(Heaviside(b - t)*(b - t)**(a - 1), [], [a], [0], [], t/b,
122 gamma(a)*b**(a - 1), And(b > 0))
123 add(Heaviside(z - (b/p)**(1/q))*(t - b)**(a - 1), [a], [], [], [0], t/b,
124 gamma(a)*b**(a - 1), And(b > 0))
125 add(Heaviside((b/p)**(1/q) - z)*(b - t)**(a - 1), [], [a], [0], [], t/b,
126 gamma(a)*b**(a - 1), And(b > 0))
127 add((b + t)**(-a), [1 - a], [], [0], [], t/b, b**(-a)/gamma(a),
128 hint=Not(IsNonPositiveInteger(a)))
129 add(Abs(b - t)**(-a), [1 - a], [(1 - a)/2], [0], [(1 - a)/2], t/b,
130 2*sin(pi*a/2)*gamma(1 - a)*Abs(b)**(-a), re(a) < 1)
131 add((t**a - b**a)/(t - b), [0, a], [], [0, a], [], t/b,
132 b**(a - 1)*sin(a*pi)/pi)
134 # 12
135 def A1(r, sign, nu):
136 return pi**Rational(-1, 2)*(-sign*nu/2)**(1 - 2*r)
138 def tmpadd(r, sgn):
139 # XXX the a**2 is bad for matching
140 add((sqrt(a**2 + t) + sgn*a)**b/(a**2 + t)**r,
141 [(1 + b)/2, 1 - 2*r + b/2], [],
142 [(b - sgn*b)/2], [(b + sgn*b)/2], t/a**2,
143 a**(b - 2*r)*A1(r, sgn, b))
144 tmpadd(0, 1)
145 tmpadd(0, -1)
146 tmpadd(S.Half, 1)
147 tmpadd(S.Half, -1)
149 # 13
150 def tmpadd(r, sgn):
151 add((sqrt(a + p*z**q) + sgn*sqrt(p)*z**(q/2))**b/(a + p*z**q)**r,
152 [1 - r + sgn*b/2], [1 - r - sgn*b/2], [0, S.Half], [],
153 p*z**q/a, a**(b/2 - r)*A1(r, sgn, b))
154 tmpadd(0, 1)
155 tmpadd(0, -1)
156 tmpadd(S.Half, 1)
157 tmpadd(S.Half, -1)
158 # (those after look obscure)
160 # Section 8.4.3
161 add(exp(polar_lift(-1)*t), [], [], [0], [])
163 # TODO can do sin^n, sinh^n by expansion ... where?
164 # 8.4.4 (hyperbolic functions)
165 add(sinh(t), [], [1], [S.Half], [1, 0], t**2/4, pi**Rational(3, 2))
166 add(cosh(t), [], [S.Half], [0], [S.Half, S.Half], t**2/4, pi**Rational(3, 2))
168 # Section 8.4.5
169 # TODO can do t + a. but can also do by expansion... (XXX not really)
170 add(sin(t), [], [], [S.Half], [0], t**2/4, sqrt(pi))
171 add(cos(t), [], [], [0], [S.Half], t**2/4, sqrt(pi))
173 # Section 8.4.6 (sinc function)
174 add(sinc(t), [], [], [0], [Rational(-1, 2)], t**2/4, sqrt(pi)/2)
176 # Section 8.5.5
177 def make_log1(subs):
178 N = subs[n]
179 return [(S.NegativeOne**N*factorial(N),
180 meijerg([], [1]*(N + 1), [0]*(N + 1), [], t))]
182 def make_log2(subs):
183 N = subs[n]
184 return [(factorial(N),
185 meijerg([1]*(N + 1), [], [], [0]*(N + 1), t))]
186 # TODO these only hold for positive p, and can be made more general
187 # but who uses log(x)*Heaviside(a-x) anyway ...
188 # TODO also it would be nice to derive them recursively ...
189 addi(log(t)**n*Heaviside(1 - t), make_log1, True)
190 addi(log(t)**n*Heaviside(t - 1), make_log2, True)
192 def make_log3(subs):
193 return make_log1(subs) + make_log2(subs)
194 addi(log(t)**n, make_log3, True)
195 addi(log(t + a),
196 constant(log(a)) + [(S.One, meijerg([1, 1], [], [1], [0], t/a))],
197 True)
198 addi(log(Abs(t - a)), constant(log(Abs(a))) +
199 [(pi, meijerg([1, 1], [S.Half], [1], [0, S.Half], t/a))],
200 True)
201 # TODO log(x)/(x+a) and log(x)/(x-1) can also be done. should they
202 # be derivable?
203 # TODO further formulae in this section seem obscure
205 # Sections 8.4.9-10
206 # TODO
208 # Section 8.4.11
209 addi(Ei(t),
210 constant(-S.ImaginaryUnit*pi) + [(S.NegativeOne, meijerg([], [1], [0, 0], [],
211 t*polar_lift(-1)))],
212 True)
214 # Section 8.4.12
215 add(Si(t), [1], [], [S.Half], [0, 0], t**2/4, sqrt(pi)/2)
216 add(Ci(t), [], [1], [0, 0], [S.Half], t**2/4, -sqrt(pi)/2)
218 # Section 8.4.13
219 add(Shi(t), [S.Half], [], [0], [Rational(-1, 2), Rational(-1, 2)], polar_lift(-1)*t**2/4,
220 t*sqrt(pi)/4)
221 add(Chi(t), [], [S.Half, 1], [0, 0], [S.Half, S.Half], t**2/4, -
222 pi**S('3/2')/2)
224 # generalized exponential integral
225 add(expint(a, t), [], [a], [a - 1, 0], [], t)
227 # Section 8.4.14
228 add(erf(t), [1], [], [S.Half], [0], t**2, 1/sqrt(pi))
229 # TODO exp(-x)*erf(I*x) does not work
230 add(erfc(t), [], [1], [0, S.Half], [], t**2, 1/sqrt(pi))
231 # This formula for erfi(z) yields a wrong(?) minus sign
232 #add(erfi(t), [1], [], [S.Half], [0], -t**2, I/sqrt(pi))
233 add(erfi(t), [S.Half], [], [0], [Rational(-1, 2)], -t**2, t/sqrt(pi))
235 # Fresnel Integrals
236 add(fresnels(t), [1], [], [Rational(3, 4)], [0, Rational(1, 4)], pi**2*t**4/16, S.Half)
237 add(fresnelc(t), [1], [], [Rational(1, 4)], [0, Rational(3, 4)], pi**2*t**4/16, S.Half)
239 ##### bessel-type functions #####
240 # Section 8.4.19
241 add(besselj(a, t), [], [], [a/2], [-a/2], t**2/4)
243 # all of the following are derivable
244 #add(sin(t)*besselj(a, t), [Rational(1, 4), Rational(3, 4)], [], [(1+a)/2],
245 # [-a/2, a/2, (1-a)/2], t**2, 1/sqrt(2))
246 #add(cos(t)*besselj(a, t), [Rational(1, 4), Rational(3, 4)], [], [a/2],
247 # [-a/2, (1+a)/2, (1-a)/2], t**2, 1/sqrt(2))
248 #add(besselj(a, t)**2, [S.Half], [], [a], [-a, 0], t**2, 1/sqrt(pi))
249 #add(besselj(a, t)*besselj(b, t), [0, S.Half], [], [(a + b)/2],
250 # [-(a+b)/2, (a - b)/2, (b - a)/2], t**2, 1/sqrt(pi))
252 # Section 8.4.20
253 add(bessely(a, t), [], [-(a + 1)/2], [a/2, -a/2], [-(a + 1)/2], t**2/4)
255 # TODO all of the following should be derivable
256 #add(sin(t)*bessely(a, t), [Rational(1, 4), Rational(3, 4)], [(1 - a - 1)/2],
257 # [(1 + a)/2, (1 - a)/2], [(1 - a - 1)/2, (1 - 1 - a)/2, (1 - 1 + a)/2],
258 # t**2, 1/sqrt(2))
259 #add(cos(t)*bessely(a, t), [Rational(1, 4), Rational(3, 4)], [(0 - a - 1)/2],
260 # [(0 + a)/2, (0 - a)/2], [(0 - a - 1)/2, (1 - 0 - a)/2, (1 - 0 + a)/2],
261 # t**2, 1/sqrt(2))
262 #add(besselj(a, t)*bessely(b, t), [0, S.Half], [(a - b - 1)/2],
263 # [(a + b)/2, (a - b)/2], [(a - b - 1)/2, -(a + b)/2, (b - a)/2],
264 # t**2, 1/sqrt(pi))
265 #addi(bessely(a, t)**2,
266 # [(2/sqrt(pi), meijerg([], [S.Half, S.Half - a], [0, a, -a],
267 # [S.Half - a], t**2)),
268 # (1/sqrt(pi), meijerg([S.Half], [], [a], [-a, 0], t**2))],
269 # True)
270 #addi(bessely(a, t)*bessely(b, t),
271 # [(2/sqrt(pi), meijerg([], [0, S.Half, (1 - a - b)/2],
272 # [(a + b)/2, (a - b)/2, (b - a)/2, -(a + b)/2],
273 # [(1 - a - b)/2], t**2)),
274 # (1/sqrt(pi), meijerg([0, S.Half], [], [(a + b)/2],
275 # [-(a + b)/2, (a - b)/2, (b - a)/2], t**2))],
276 # True)
278 # Section 8.4.21 ?
279 # Section 8.4.22
280 add(besseli(a, t), [], [(1 + a)/2], [a/2], [-a/2, (1 + a)/2], t**2/4, pi)
281 # TODO many more formulas. should all be derivable
283 # Section 8.4.23
284 add(besselk(a, t), [], [], [a/2, -a/2], [], t**2/4, S.Half)
285 # TODO many more formulas. should all be derivable
287 # Complete elliptic integrals K(z) and E(z)
288 add(elliptic_k(t), [S.Half, S.Half], [], [0], [0], -t, S.Half)
289 add(elliptic_e(t), [S.Half, 3*S.Half], [], [0], [0], -t, Rational(-1, 2)/2)
292####################################################################
293# First some helper functions.
294####################################################################
296from sympy.utilities.timeutils import timethis
297timeit = timethis('meijerg')
300def _mytype(f: Basic, x: Symbol) -> tuple[type[Basic], ...]:
301 """ Create a hashable entity describing the type of f. """
302 def key(x: type[Basic]) -> tuple[int, int, str]:
303 return x.class_key()
305 if x not in f.free_symbols:
306 return ()
307 elif f.is_Function:
308 return type(f),
309 return tuple(sorted((t for a in f.args for t in _mytype(a, x)), key=key))
312class _CoeffExpValueError(ValueError):
313 """
314 Exception raised by _get_coeff_exp, for internal use only.
315 """
316 pass
319def _get_coeff_exp(expr, x):
320 """
321 When expr is known to be of the form c*x**b, with c and/or b possibly 1,
322 return c, b.
324 Examples
325 ========
327 >>> from sympy.abc import x, a, b
328 >>> from sympy.integrals.meijerint import _get_coeff_exp
329 >>> _get_coeff_exp(a*x**b, x)
330 (a, b)
331 >>> _get_coeff_exp(x, x)
332 (1, 1)
333 >>> _get_coeff_exp(2*x, x)
334 (2, 1)
335 >>> _get_coeff_exp(x**3, x)
336 (1, 3)
337 """
338 from sympy.simplify import powsimp
339 (c, m) = expand_power_base(powsimp(expr)).as_coeff_mul(x)
340 if not m:
341 return c, S.Zero
342 [m] = m
343 if m.is_Pow:
344 if m.base != x:
345 raise _CoeffExpValueError('expr not of form a*x**b')
346 return c, m.exp
347 elif m == x:
348 return c, S.One
349 else:
350 raise _CoeffExpValueError('expr not of form a*x**b: %s' % expr)
353def _exponents(expr, x):
354 """
355 Find the exponents of ``x`` (not including zero) in ``expr``.
357 Examples
358 ========
360 >>> from sympy.integrals.meijerint import _exponents
361 >>> from sympy.abc import x, y
362 >>> from sympy import sin
363 >>> _exponents(x, x)
364 {1}
365 >>> _exponents(x**2, x)
366 {2}
367 >>> _exponents(x**2 + x, x)
368 {1, 2}
369 >>> _exponents(x**3*sin(x + x**y) + 1/x, x)
370 {-1, 1, 3, y}
371 """
372 def _exponents_(expr, x, res):
373 if expr == x:
374 res.update([1])
375 return
376 if expr.is_Pow and expr.base == x:
377 res.update([expr.exp])
378 return
379 for argument in expr.args:
380 _exponents_(argument, x, res)
381 res = set()
382 _exponents_(expr, x, res)
383 return res
386def _functions(expr, x):
387 """ Find the types of functions in expr, to estimate the complexity. """
388 return {e.func for e in expr.atoms(Function) if x in e.free_symbols}
391def _find_splitting_points(expr, x):
392 """
393 Find numbers a such that a linear substitution x -> x + a would
394 (hopefully) simplify expr.
396 Examples
397 ========
399 >>> from sympy.integrals.meijerint import _find_splitting_points as fsp
400 >>> from sympy import sin
401 >>> from sympy.abc import x
402 >>> fsp(x, x)
403 {0}
404 >>> fsp((x-1)**3, x)
405 {1}
406 >>> fsp(sin(x+3)*x, x)
407 {-3, 0}
408 """
409 p, q = [Wild(n, exclude=[x]) for n in 'pq']
411 def compute_innermost(expr, res):
412 if not isinstance(expr, Expr):
413 return
414 m = expr.match(p*x + q)
415 if m and m[p] != 0:
416 res.add(-m[q]/m[p])
417 return
418 if expr.is_Atom:
419 return
420 for argument in expr.args:
421 compute_innermost(argument, res)
422 innermost = set()
423 compute_innermost(expr, innermost)
424 return innermost
427def _split_mul(f, x):
428 """
429 Split expression ``f`` into fac, po, g, where fac is a constant factor,
430 po = x**s for some s independent of s, and g is "the rest".
432 Examples
433 ========
435 >>> from sympy.integrals.meijerint import _split_mul
436 >>> from sympy import sin
437 >>> from sympy.abc import s, x
438 >>> _split_mul((3*x)**s*sin(x**2)*x, x)
439 (3**s, x*x**s, sin(x**2))
440 """
441 fac = S.One
442 po = S.One
443 g = S.One
444 f = expand_power_base(f)
446 args = Mul.make_args(f)
447 for a in args:
448 if a == x:
449 po *= x
450 elif x not in a.free_symbols:
451 fac *= a
452 else:
453 if a.is_Pow and x not in a.exp.free_symbols:
454 c, t = a.base.as_coeff_mul(x)
455 if t != (x,):
456 c, t = expand_mul(a.base).as_coeff_mul(x)
457 if t == (x,):
458 po *= x**a.exp
459 fac *= unpolarify(polarify(c**a.exp, subs=False))
460 continue
461 g *= a
463 return fac, po, g
466def _mul_args(f):
467 """
468 Return a list ``L`` such that ``Mul(*L) == f``.
470 If ``f`` is not a ``Mul`` or ``Pow``, ``L=[f]``.
471 If ``f=g**n`` for an integer ``n``, ``L=[g]*n``.
472 If ``f`` is a ``Mul``, ``L`` comes from applying ``_mul_args`` to all factors of ``f``.
473 """
474 args = Mul.make_args(f)
475 gs = []
476 for g in args:
477 if g.is_Pow and g.exp.is_Integer:
478 n = g.exp
479 base = g.base
480 if n < 0:
481 n = -n
482 base = 1/base
483 gs += [base]*n
484 else:
485 gs.append(g)
486 return gs
489def _mul_as_two_parts(f):
490 """
491 Find all the ways to split ``f`` into a product of two terms.
492 Return None on failure.
494 Explanation
495 ===========
497 Although the order is canonical from multiset_partitions, this is
498 not necessarily the best order to process the terms. For example,
499 if the case of len(gs) == 2 is removed and multiset is allowed to
500 sort the terms, some tests fail.
502 Examples
503 ========
505 >>> from sympy.integrals.meijerint import _mul_as_two_parts
506 >>> from sympy import sin, exp, ordered
507 >>> from sympy.abc import x
508 >>> list(ordered(_mul_as_two_parts(x*sin(x)*exp(x))))
509 [(x, exp(x)*sin(x)), (x*exp(x), sin(x)), (x*sin(x), exp(x))]
510 """
512 gs = _mul_args(f)
513 if len(gs) < 2:
514 return None
515 if len(gs) == 2:
516 return [tuple(gs)]
517 return [(Mul(*x), Mul(*y)) for (x, y) in multiset_partitions(gs, 2)]
520def _inflate_g(g, n):
521 """ Return C, h such that h is a G function of argument z**n and
522 g = C*h. """
523 # TODO should this be a method of meijerg?
524 # See: [L, page 150, equation (5)]
525 def inflate(params, n):
526 """ (a1, .., ak) -> (a1/n, (a1+1)/n, ..., (ak + n-1)/n) """
527 return [(a + i)/n for a, i in itertools.product(params, range(n))]
528 v = S(len(g.ap) - len(g.bq))
529 C = n**(1 + g.nu + v/2)
530 C /= (2*pi)**((n - 1)*g.delta)
531 return C, meijerg(inflate(g.an, n), inflate(g.aother, n),
532 inflate(g.bm, n), inflate(g.bother, n),
533 g.argument**n * n**(n*v))
536def _flip_g(g):
537 """ Turn the G function into one of inverse argument
538 (i.e. G(1/x) -> G'(x)) """
539 # See [L], section 5.2
540 def tr(l):
541 return [1 - a for a in l]
542 return meijerg(tr(g.bm), tr(g.bother), tr(g.an), tr(g.aother), 1/g.argument)
545def _inflate_fox_h(g, a):
546 r"""
547 Let d denote the integrand in the definition of the G function ``g``.
548 Consider the function H which is defined in the same way, but with
549 integrand d/Gamma(a*s) (contour conventions as usual).
551 If ``a`` is rational, the function H can be written as C*G, for a constant C
552 and a G-function G.
554 This function returns C, G.
555 """
556 if a < 0:
557 return _inflate_fox_h(_flip_g(g), -a)
558 p = S(a.p)
559 q = S(a.q)
560 # We use the substitution s->qs, i.e. inflate g by q. We are left with an
561 # extra factor of Gamma(p*s), for which we use Gauss' multiplication
562 # theorem.
563 D, g = _inflate_g(g, q)
564 z = g.argument
565 D /= (2*pi)**((1 - p)/2)*p**Rational(-1, 2)
566 z /= p**p
567 bs = [(n + 1)/p for n in range(p)]
568 return D, meijerg(g.an, g.aother, g.bm, list(g.bother) + bs, z)
570_dummies: dict[tuple[str, str], Dummy] = {}
573def _dummy(name, token, expr, **kwargs):
574 """
575 Return a dummy. This will return the same dummy if the same token+name is
576 requested more than once, and it is not already in expr.
577 This is for being cache-friendly.
578 """
579 d = _dummy_(name, token, **kwargs)
580 if d in expr.free_symbols:
581 return Dummy(name, **kwargs)
582 return d
585def _dummy_(name, token, **kwargs):
586 """
587 Return a dummy associated to name and token. Same effect as declaring
588 it globally.
589 """
590 global _dummies
591 if not (name, token) in _dummies:
592 _dummies[(name, token)] = Dummy(name, **kwargs)
593 return _dummies[(name, token)]
596def _is_analytic(f, x):
597 """ Check if f(x), when expressed using G functions on the positive reals,
598 will in fact agree with the G functions almost everywhere """
599 return not any(x in expr.free_symbols for expr in f.atoms(Heaviside, Abs))
602def _condsimp(cond, first=True):
603 """
604 Do naive simplifications on ``cond``.
606 Explanation
607 ===========
609 Note that this routine is completely ad-hoc, simplification rules being
610 added as need arises rather than following any logical pattern.
612 Examples
613 ========
615 >>> from sympy.integrals.meijerint import _condsimp as simp
616 >>> from sympy import Or, Eq
617 >>> from sympy.abc import x, y
618 >>> simp(Or(x < y, Eq(x, y)))
619 x <= y
620 """
621 if first:
622 cond = cond.replace(lambda _: _.is_Relational, _canonical_coeff)
623 first = False
624 if not isinstance(cond, BooleanFunction):
625 return cond
626 p, q, r = symbols('p q r', cls=Wild)
627 # transforms tests use 0, 4, 5 and 11-14
628 # meijer tests use 0, 2, 11, 14
629 # joint_rv uses 6, 7
630 rules = [
631 (Or(p < q, Eq(p, q)), p <= q), # 0
632 # The next two obviously are instances of a general pattern, but it is
633 # easier to spell out the few cases we care about.
634 (And(Abs(arg(p)) <= pi, Abs(arg(p) - 2*pi) <= pi),
635 Eq(arg(p) - pi, 0)), # 1
636 (And(Abs(2*arg(p) + pi) <= pi, Abs(2*arg(p) - pi) <= pi),
637 Eq(arg(p), 0)), # 2
638 (And(Abs(2*arg(p) + pi) < pi, Abs(2*arg(p) - pi) <= pi),
639 S.false), # 3
640 (And(Abs(arg(p) - pi/2) <= pi/2, Abs(arg(p) + pi/2) <= pi/2),
641 Eq(arg(p), 0)), # 4
642 (And(Abs(arg(p) - pi/2) <= pi/2, Abs(arg(p) + pi/2) < pi/2),
643 S.false), # 5
644 (And(Abs(arg(p**2/2 + 1)) < pi, Ne(Abs(arg(p**2/2 + 1)), pi)),
645 S.true), # 6
646 (Or(Abs(arg(p**2/2 + 1)) < pi, Ne(1/(p**2/2 + 1), 0)),
647 S.true), # 7
648 (And(Abs(unbranched_argument(p)) <= pi,
649 Abs(unbranched_argument(exp_polar(-2*pi*S.ImaginaryUnit)*p)) <= pi),
650 Eq(unbranched_argument(exp_polar(-S.ImaginaryUnit*pi)*p), 0)), # 8
651 (And(Abs(unbranched_argument(p)) <= pi/2,
652 Abs(unbranched_argument(exp_polar(-pi*S.ImaginaryUnit)*p)) <= pi/2),
653 Eq(unbranched_argument(exp_polar(-S.ImaginaryUnit*pi/2)*p), 0)), # 9
654 (Or(p <= q, And(p < q, r)), p <= q), # 10
655 (Ne(p**2, 1) & (p**2 > 1), p**2 > 1), # 11
656 (Ne(1/p, 1) & (cos(Abs(arg(p)))*Abs(p) > 1), Abs(p) > 1), # 12
657 (Ne(p, 2) & (cos(Abs(arg(p)))*Abs(p) > 2), Abs(p) > 2), # 13
658 ((Abs(arg(p)) < pi/2) & (cos(Abs(arg(p)))*sqrt(Abs(p**2)) > 1), p**2 > 1), # 14
659 ]
660 cond = cond.func(*[_condsimp(_, first) for _ in cond.args])
661 change = True
662 while change:
663 change = False
664 for irule, (fro, to) in enumerate(rules):
665 if fro.func != cond.func:
666 continue
667 for n, arg1 in enumerate(cond.args):
668 if r in fro.args[0].free_symbols:
669 m = arg1.match(fro.args[1])
670 num = 1
671 else:
672 num = 0
673 m = arg1.match(fro.args[0])
674 if not m:
675 continue
676 otherargs = [x.subs(m) for x in fro.args[:num] + fro.args[num + 1:]]
677 otherlist = [n]
678 for arg2 in otherargs:
679 for k, arg3 in enumerate(cond.args):
680 if k in otherlist:
681 continue
682 if arg2 == arg3:
683 otherlist += [k]
684 break
685 if isinstance(arg3, And) and arg2.args[1] == r and \
686 isinstance(arg2, And) and arg2.args[0] in arg3.args:
687 otherlist += [k]
688 break
689 if isinstance(arg3, And) and arg2.args[0] == r and \
690 isinstance(arg2, And) and arg2.args[1] in arg3.args:
691 otherlist += [k]
692 break
693 if len(otherlist) != len(otherargs) + 1:
694 continue
695 newargs = [arg_ for (k, arg_) in enumerate(cond.args)
696 if k not in otherlist] + [to.subs(m)]
697 if SYMPY_DEBUG:
698 if irule not in (0, 2, 4, 5, 6, 7, 11, 12, 13, 14):
699 print('used new rule:', irule)
700 cond = cond.func(*newargs)
701 change = True
702 break
704 # final tweak
705 def rel_touchup(rel):
706 if rel.rel_op != '==' or rel.rhs != 0:
707 return rel
709 # handle Eq(*, 0)
710 LHS = rel.lhs
711 m = LHS.match(arg(p)**q)
712 if not m:
713 m = LHS.match(unbranched_argument(polar_lift(p)**q))
714 if not m:
715 if isinstance(LHS, periodic_argument) and not LHS.args[0].is_polar \
716 and LHS.args[1] is S.Infinity:
717 return (LHS.args[0] > 0)
718 return rel
719 return (m[p] > 0)
720 cond = cond.replace(lambda _: _.is_Relational, rel_touchup)
721 if SYMPY_DEBUG:
722 print('_condsimp: ', cond)
723 return cond
725def _eval_cond(cond):
726 """ Re-evaluate the conditions. """
727 if isinstance(cond, bool):
728 return cond
729 return _condsimp(cond.doit())
731####################################################################
732# Now the "backbone" functions to do actual integration.
733####################################################################
736def _my_principal_branch(expr, period, full_pb=False):
737 """ Bring expr nearer to its principal branch by removing superfluous
738 factors.
739 This function does *not* guarantee to yield the principal branch,
740 to avoid introducing opaque principal_branch() objects,
741 unless full_pb=True. """
742 res = principal_branch(expr, period)
743 if not full_pb:
744 res = res.replace(principal_branch, lambda x, y: x)
745 return res
748def _rewrite_saxena_1(fac, po, g, x):
749 """
750 Rewrite the integral fac*po*g dx, from zero to infinity, as
751 integral fac*G, where G has argument a*x. Note po=x**s.
752 Return fac, G.
753 """
754 _, s = _get_coeff_exp(po, x)
755 a, b = _get_coeff_exp(g.argument, x)
756 period = g.get_period()
757 a = _my_principal_branch(a, period)
759 # We substitute t = x**b.
760 C = fac/(Abs(b)*a**((s + 1)/b - 1))
761 # Absorb a factor of (at)**((1 + s)/b - 1).
763 def tr(l):
764 return [a + (1 + s)/b - 1 for a in l]
765 return C, meijerg(tr(g.an), tr(g.aother), tr(g.bm), tr(g.bother),
766 a*x)
769def _check_antecedents_1(g, x, helper=False):
770 r"""
771 Return a condition under which the mellin transform of g exists.
772 Any power of x has already been absorbed into the G function,
773 so this is just $\int_0^\infty g\, dx$.
775 See [L, section 5.6.1]. (Note that s=1.)
777 If ``helper`` is True, only check if the MT exists at infinity, i.e. if
778 $\int_1^\infty g\, dx$ exists.
779 """
780 # NOTE if you update these conditions, please update the documentation as well
781 delta = g.delta
782 eta, _ = _get_coeff_exp(g.argument, x)
783 m, n, p, q = S([len(g.bm), len(g.an), len(g.ap), len(g.bq)])
785 if p > q:
786 def tr(l):
787 return [1 - x for x in l]
788 return _check_antecedents_1(meijerg(tr(g.bm), tr(g.bother),
789 tr(g.an), tr(g.aother), x/eta),
790 x)
792 tmp = [-re(b) < 1 for b in g.bm] + [1 < 1 - re(a) for a in g.an]
793 cond_3 = And(*tmp)
795 tmp += [-re(b) < 1 for b in g.bother]
796 tmp += [1 < 1 - re(a) for a in g.aother]
797 cond_3_star = And(*tmp)
799 cond_4 = (-re(g.nu) + (q + 1 - p)/2 > q - p)
801 def debug(*msg):
802 _debug(*msg)
804 def debugf(string, arg):
805 _debugf(string, arg)
807 debug('Checking antecedents for 1 function:')
808 debugf(' delta=%s, eta=%s, m=%s, n=%s, p=%s, q=%s',
809 (delta, eta, m, n, p, q))
810 debugf(' ap = %s, %s', (list(g.an), list(g.aother)))
811 debugf(' bq = %s, %s', (list(g.bm), list(g.bother)))
812 debugf(' cond_3=%s, cond_3*=%s, cond_4=%s', (cond_3, cond_3_star, cond_4))
814 conds = []
816 # case 1
817 case1 = []
818 tmp1 = [1 <= n, p < q, 1 <= m]
819 tmp2 = [1 <= p, 1 <= m, Eq(q, p + 1), Not(And(Eq(n, 0), Eq(m, p + 1)))]
820 tmp3 = [1 <= p, Eq(q, p)]
821 for k in range(ceiling(delta/2) + 1):
822 tmp3 += [Ne(Abs(unbranched_argument(eta)), (delta - 2*k)*pi)]
823 tmp = [delta > 0, Abs(unbranched_argument(eta)) < delta*pi]
824 extra = [Ne(eta, 0), cond_3]
825 if helper:
826 extra = []
827 for t in [tmp1, tmp2, tmp3]:
828 case1 += [And(*(t + tmp + extra))]
829 conds += case1
830 debug(' case 1:', case1)
832 # case 2
833 extra = [cond_3]
834 if helper:
835 extra = []
836 case2 = [And(Eq(n, 0), p + 1 <= m, m <= q,
837 Abs(unbranched_argument(eta)) < delta*pi, *extra)]
838 conds += case2
839 debug(' case 2:', case2)
841 # case 3
842 extra = [cond_3, cond_4]
843 if helper:
844 extra = []
845 case3 = [And(p < q, 1 <= m, delta > 0, Eq(Abs(unbranched_argument(eta)), delta*pi),
846 *extra)]
847 case3 += [And(p <= q - 2, Eq(delta, 0), Eq(Abs(unbranched_argument(eta)), 0), *extra)]
848 conds += case3
849 debug(' case 3:', case3)
851 # TODO altered cases 4-7
853 # extra case from wofram functions site:
854 # (reproduced verbatim from Prudnikov, section 2.24.2)
855 # https://functions.wolfram.com/HypergeometricFunctions/MeijerG/21/02/01/
856 case_extra = []
857 case_extra += [Eq(p, q), Eq(delta, 0), Eq(unbranched_argument(eta), 0), Ne(eta, 0)]
858 if not helper:
859 case_extra += [cond_3]
860 s = []
861 for a, b in zip(g.ap, g.bq):
862 s += [b - a]
863 case_extra += [re(Add(*s)) < 0]
864 case_extra = And(*case_extra)
865 conds += [case_extra]
866 debug(' extra case:', [case_extra])
868 case_extra_2 = [And(delta > 0, Abs(unbranched_argument(eta)) < delta*pi)]
869 if not helper:
870 case_extra_2 += [cond_3]
871 case_extra_2 = And(*case_extra_2)
872 conds += [case_extra_2]
873 debug(' second extra case:', [case_extra_2])
875 # TODO This leaves only one case from the three listed by Prudnikov.
876 # Investigate if these indeed cover everything; if so, remove the rest.
878 return Or(*conds)
881def _int0oo_1(g, x):
882 r"""
883 Evaluate $\int_0^\infty g\, dx$ using G functions,
884 assuming the necessary conditions are fulfilled.
886 Examples
887 ========
889 >>> from sympy.abc import a, b, c, d, x, y
890 >>> from sympy import meijerg
891 >>> from sympy.integrals.meijerint import _int0oo_1
892 >>> _int0oo_1(meijerg([a], [b], [c], [d], x*y), x)
893 gamma(-a)*gamma(c + 1)/(y*gamma(-d)*gamma(b + 1))
894 """
895 from sympy.simplify import gammasimp
896 # See [L, section 5.6.1]. Note that s=1.
897 eta, _ = _get_coeff_exp(g.argument, x)
898 res = 1/eta
899 # XXX TODO we should reduce order first
900 for b in g.bm:
901 res *= gamma(b + 1)
902 for a in g.an:
903 res *= gamma(1 - a - 1)
904 for b in g.bother:
905 res /= gamma(1 - b - 1)
906 for a in g.aother:
907 res /= gamma(a + 1)
908 return gammasimp(unpolarify(res))
911def _rewrite_saxena(fac, po, g1, g2, x, full_pb=False):
912 """
913 Rewrite the integral ``fac*po*g1*g2`` from 0 to oo in terms of G
914 functions with argument ``c*x``.
916 Explanation
917 ===========
919 Return C, f1, f2 such that integral C f1 f2 from 0 to infinity equals
920 integral fac ``po``, ``g1``, ``g2`` from 0 to infinity.
922 Examples
923 ========
925 >>> from sympy.integrals.meijerint import _rewrite_saxena
926 >>> from sympy.abc import s, t, m
927 >>> from sympy import meijerg
928 >>> g1 = meijerg([], [], [0], [], s*t)
929 >>> g2 = meijerg([], [], [m/2], [-m/2], t**2/4)
930 >>> r = _rewrite_saxena(1, t**0, g1, g2, t)
931 >>> r[0]
932 s/(4*sqrt(pi))
933 >>> r[1]
934 meijerg(((), ()), ((-1/2, 0), ()), s**2*t/4)
935 >>> r[2]
936 meijerg(((), ()), ((m/2,), (-m/2,)), t/4)
937 """
938 def pb(g):
939 a, b = _get_coeff_exp(g.argument, x)
940 per = g.get_period()
941 return meijerg(g.an, g.aother, g.bm, g.bother,
942 _my_principal_branch(a, per, full_pb)*x**b)
944 _, s = _get_coeff_exp(po, x)
945 _, b1 = _get_coeff_exp(g1.argument, x)
946 _, b2 = _get_coeff_exp(g2.argument, x)
947 if (b1 < 0) == True:
948 b1 = -b1
949 g1 = _flip_g(g1)
950 if (b2 < 0) == True:
951 b2 = -b2
952 g2 = _flip_g(g2)
953 if not b1.is_Rational or not b2.is_Rational:
954 return
955 m1, n1 = b1.p, b1.q
956 m2, n2 = b2.p, b2.q
957 tau = ilcm(m1*n2, m2*n1)
958 r1 = tau//(m1*n2)
959 r2 = tau//(m2*n1)
961 C1, g1 = _inflate_g(g1, r1)
962 C2, g2 = _inflate_g(g2, r2)
963 g1 = pb(g1)
964 g2 = pb(g2)
966 fac *= C1*C2
967 a1, b = _get_coeff_exp(g1.argument, x)
968 a2, _ = _get_coeff_exp(g2.argument, x)
970 # arbitrarily tack on the x**s part to g1
971 # TODO should we try both?
972 exp = (s + 1)/b - 1
973 fac = fac/(Abs(b) * a1**exp)
975 def tr(l):
976 return [a + exp for a in l]
977 g1 = meijerg(tr(g1.an), tr(g1.aother), tr(g1.bm), tr(g1.bother), a1*x)
978 g2 = meijerg(g2.an, g2.aother, g2.bm, g2.bother, a2*x)
980 from sympy.simplify import powdenest
981 return powdenest(fac, polar=True), g1, g2
984def _check_antecedents(g1, g2, x):
985 """ Return a condition under which the integral theorem applies. """
986 # Yes, this is madness.
987 # XXX TODO this is a testing *nightmare*
988 # NOTE if you update these conditions, please update the documentation as well
990 # The following conditions are found in
991 # [P], Section 2.24.1
992 #
993 # They are also reproduced (verbatim!) at
994 # https://functions.wolfram.com/HypergeometricFunctions/MeijerG/21/02/03/
995 #
996 # Note: k=l=r=alpha=1
997 sigma, _ = _get_coeff_exp(g1.argument, x)
998 omega, _ = _get_coeff_exp(g2.argument, x)
999 s, t, u, v = S([len(g1.bm), len(g1.an), len(g1.ap), len(g1.bq)])
1000 m, n, p, q = S([len(g2.bm), len(g2.an), len(g2.ap), len(g2.bq)])
1001 bstar = s + t - (u + v)/2
1002 cstar = m + n - (p + q)/2
1003 rho = g1.nu + (u - v)/2 + 1
1004 mu = g2.nu + (p - q)/2 + 1
1005 phi = q - p - (v - u)
1006 eta = 1 - (v - u) - mu - rho
1007 psi = (pi*(q - m - n) + Abs(unbranched_argument(omega)))/(q - p)
1008 theta = (pi*(v - s - t) + Abs(unbranched_argument(sigma)))/(v - u)
1010 _debug('Checking antecedents:')
1011 _debugf(' sigma=%s, s=%s, t=%s, u=%s, v=%s, b*=%s, rho=%s',
1012 (sigma, s, t, u, v, bstar, rho))
1013 _debugf(' omega=%s, m=%s, n=%s, p=%s, q=%s, c*=%s, mu=%s,',
1014 (omega, m, n, p, q, cstar, mu))
1015 _debugf(' phi=%s, eta=%s, psi=%s, theta=%s', (phi, eta, psi, theta))
1017 def _c1():
1018 for g in [g1, g2]:
1019 for i, j in itertools.product(g.an, g.bm):
1020 diff = i - j
1021 if diff.is_integer and diff.is_positive:
1022 return False
1023 return True
1024 c1 = _c1()
1025 c2 = And(*[re(1 + i + j) > 0 for i in g1.bm for j in g2.bm])
1026 c3 = And(*[re(1 + i + j) < 1 + 1 for i in g1.an for j in g2.an])
1027 c4 = And(*[(p - q)*re(1 + i - 1) - re(mu) > Rational(-3, 2) for i in g1.an])
1028 c5 = And(*[(p - q)*re(1 + i) - re(mu) > Rational(-3, 2) for i in g1.bm])
1029 c6 = And(*[(u - v)*re(1 + i - 1) - re(rho) > Rational(-3, 2) for i in g2.an])
1030 c7 = And(*[(u - v)*re(1 + i) - re(rho) > Rational(-3, 2) for i in g2.bm])
1031 c8 = (Abs(phi) + 2*re((rho - 1)*(q - p) + (v - u)*(q - p) + (mu -
1032 1)*(v - u)) > 0)
1033 c9 = (Abs(phi) - 2*re((rho - 1)*(q - p) + (v - u)*(q - p) + (mu -
1034 1)*(v - u)) > 0)
1035 c10 = (Abs(unbranched_argument(sigma)) < bstar*pi)
1036 c11 = Eq(Abs(unbranched_argument(sigma)), bstar*pi)
1037 c12 = (Abs(unbranched_argument(omega)) < cstar*pi)
1038 c13 = Eq(Abs(unbranched_argument(omega)), cstar*pi)
1040 # The following condition is *not* implemented as stated on the wolfram
1041 # function site. In the book of Prudnikov there is an additional part
1042 # (the And involving re()). However, I only have this book in russian, and
1043 # I don't read any russian. The following condition is what other people
1044 # have told me it means.
1045 # Worryingly, it is different from the condition implemented in REDUCE.
1046 # The REDUCE implementation:
1047 # https://reduce-algebra.svn.sourceforge.net/svnroot/reduce-algebra/trunk/packages/defint/definta.red
1048 # (search for tst14)
1049 # The Wolfram alpha version:
1050 # https://functions.wolfram.com/HypergeometricFunctions/MeijerG/21/02/03/03/0014/
1051 z0 = exp(-(bstar + cstar)*pi*S.ImaginaryUnit)
1052 zos = unpolarify(z0*omega/sigma)
1053 zso = unpolarify(z0*sigma/omega)
1054 if zos == 1/zso:
1055 c14 = And(Eq(phi, 0), bstar + cstar <= 1,
1056 Or(Ne(zos, 1), re(mu + rho + v - u) < 1,
1057 re(mu + rho + q - p) < 1))
1058 else:
1059 def _cond(z):
1060 '''Returns True if abs(arg(1-z)) < pi, avoiding arg(0).
1062 Explanation
1063 ===========
1065 If ``z`` is 1 then arg is NaN. This raises a
1066 TypeError on `NaN < pi`. Previously this gave `False` so
1067 this behavior has been hardcoded here but someone should
1068 check if this NaN is more serious! This NaN is triggered by
1069 test_meijerint() in test_meijerint.py:
1070 `meijerint_definite(exp(x), x, 0, I)`
1071 '''
1072 return z != 1 and Abs(arg(1 - z)) < pi
1074 c14 = And(Eq(phi, 0), bstar - 1 + cstar <= 0,
1075 Or(And(Ne(zos, 1), _cond(zos)),
1076 And(re(mu + rho + v - u) < 1, Eq(zos, 1))))
1078 c14_alt = And(Eq(phi, 0), cstar - 1 + bstar <= 0,
1079 Or(And(Ne(zso, 1), _cond(zso)),
1080 And(re(mu + rho + q - p) < 1, Eq(zso, 1))))
1082 # Since r=k=l=1, in our case there is c14_alt which is the same as calling
1083 # us with (g1, g2) = (g2, g1). The conditions below enumerate all cases
1084 # (i.e. we don't have to try arguments reversed by hand), and indeed try
1085 # all symmetric cases. (i.e. whenever there is a condition involving c14,
1086 # there is also a dual condition which is exactly what we would get when g1,
1087 # g2 were interchanged, *but c14 was unaltered*).
1088 # Hence the following seems correct:
1089 c14 = Or(c14, c14_alt)
1091 '''
1092 When `c15` is NaN (e.g. from `psi` being NaN as happens during
1093 'test_issue_4992' and/or `theta` is NaN as in 'test_issue_6253',
1094 both in `test_integrals.py`) the comparison to 0 formerly gave False
1095 whereas now an error is raised. To keep the old behavior, the value
1096 of NaN is replaced with False but perhaps a closer look at this condition
1097 should be made: XXX how should conditions leading to c15=NaN be handled?
1098 '''
1099 try:
1100 lambda_c = (q - p)*Abs(omega)**(1/(q - p))*cos(psi) \
1101 + (v - u)*Abs(sigma)**(1/(v - u))*cos(theta)
1102 # the TypeError might be raised here, e.g. if lambda_c is NaN
1103 if _eval_cond(lambda_c > 0) != False:
1104 c15 = (lambda_c > 0)
1105 else:
1106 def lambda_s0(c1, c2):
1107 return c1*(q - p)*Abs(omega)**(1/(q - p))*sin(psi) \
1108 + c2*(v - u)*Abs(sigma)**(1/(v - u))*sin(theta)
1109 lambda_s = Piecewise(
1110 ((lambda_s0(+1, +1)*lambda_s0(-1, -1)),
1111 And(Eq(unbranched_argument(sigma), 0), Eq(unbranched_argument(omega), 0))),
1112 (lambda_s0(sign(unbranched_argument(omega)), +1)*lambda_s0(sign(unbranched_argument(omega)), -1),
1113 And(Eq(unbranched_argument(sigma), 0), Ne(unbranched_argument(omega), 0))),
1114 (lambda_s0(+1, sign(unbranched_argument(sigma)))*lambda_s0(-1, sign(unbranched_argument(sigma))),
1115 And(Ne(unbranched_argument(sigma), 0), Eq(unbranched_argument(omega), 0))),
1116 (lambda_s0(sign(unbranched_argument(omega)), sign(unbranched_argument(sigma))), True))
1117 tmp = [lambda_c > 0,
1118 And(Eq(lambda_c, 0), Ne(lambda_s, 0), re(eta) > -1),
1119 And(Eq(lambda_c, 0), Eq(lambda_s, 0), re(eta) > 0)]
1120 c15 = Or(*tmp)
1121 except TypeError:
1122 c15 = False
1123 for cond, i in [(c1, 1), (c2, 2), (c3, 3), (c4, 4), (c5, 5), (c6, 6),
1124 (c7, 7), (c8, 8), (c9, 9), (c10, 10), (c11, 11),
1125 (c12, 12), (c13, 13), (c14, 14), (c15, 15)]:
1126 _debugf(' c%s: %s', (i, cond))
1128 # We will return Or(*conds)
1129 conds = []
1131 def pr(count):
1132 _debugf(' case %s: %s', (count, conds[-1]))
1133 conds += [And(m*n*s*t != 0, bstar.is_positive is True, cstar.is_positive is True, c1, c2, c3, c10,
1134 c12)] # 1
1135 pr(1)
1136 conds += [And(Eq(u, v), Eq(bstar, 0), cstar.is_positive is True, sigma.is_positive is True, re(rho) < 1,
1137 c1, c2, c3, c12)] # 2
1138 pr(2)
1139 conds += [And(Eq(p, q), Eq(cstar, 0), bstar.is_positive is True, omega.is_positive is True, re(mu) < 1,
1140 c1, c2, c3, c10)] # 3
1141 pr(3)
1142 conds += [And(Eq(p, q), Eq(u, v), Eq(bstar, 0), Eq(cstar, 0),
1143 sigma.is_positive is True, omega.is_positive is True, re(mu) < 1, re(rho) < 1,
1144 Ne(sigma, omega), c1, c2, c3)] # 4
1145 pr(4)
1146 conds += [And(Eq(p, q), Eq(u, v), Eq(bstar, 0), Eq(cstar, 0),
1147 sigma.is_positive is True, omega.is_positive is True, re(mu + rho) < 1,
1148 Ne(omega, sigma), c1, c2, c3)] # 5
1149 pr(5)
1150 conds += [And(p > q, s.is_positive is True, bstar.is_positive is True, cstar >= 0,
1151 c1, c2, c3, c5, c10, c13)] # 6
1152 pr(6)
1153 conds += [And(p < q, t.is_positive is True, bstar.is_positive is True, cstar >= 0,
1154 c1, c2, c3, c4, c10, c13)] # 7
1155 pr(7)
1156 conds += [And(u > v, m.is_positive is True, cstar.is_positive is True, bstar >= 0,
1157 c1, c2, c3, c7, c11, c12)] # 8
1158 pr(8)
1159 conds += [And(u < v, n.is_positive is True, cstar.is_positive is True, bstar >= 0,
1160 c1, c2, c3, c6, c11, c12)] # 9
1161 pr(9)
1162 conds += [And(p > q, Eq(u, v), Eq(bstar, 0), cstar >= 0, sigma.is_positive is True,
1163 re(rho) < 1, c1, c2, c3, c5, c13)] # 10
1164 pr(10)
1165 conds += [And(p < q, Eq(u, v), Eq(bstar, 0), cstar >= 0, sigma.is_positive is True,
1166 re(rho) < 1, c1, c2, c3, c4, c13)] # 11
1167 pr(11)
1168 conds += [And(Eq(p, q), u > v, bstar >= 0, Eq(cstar, 0), omega.is_positive is True,
1169 re(mu) < 1, c1, c2, c3, c7, c11)] # 12
1170 pr(12)
1171 conds += [And(Eq(p, q), u < v, bstar >= 0, Eq(cstar, 0), omega.is_positive is True,
1172 re(mu) < 1, c1, c2, c3, c6, c11)] # 13
1173 pr(13)
1174 conds += [And(p < q, u > v, bstar >= 0, cstar >= 0,
1175 c1, c2, c3, c4, c7, c11, c13)] # 14
1176 pr(14)
1177 conds += [And(p > q, u < v, bstar >= 0, cstar >= 0,
1178 c1, c2, c3, c5, c6, c11, c13)] # 15
1179 pr(15)
1180 conds += [And(p > q, u > v, bstar >= 0, cstar >= 0,
1181 c1, c2, c3, c5, c7, c8, c11, c13, c14)] # 16
1182 pr(16)
1183 conds += [And(p < q, u < v, bstar >= 0, cstar >= 0,
1184 c1, c2, c3, c4, c6, c9, c11, c13, c14)] # 17
1185 pr(17)
1186 conds += [And(Eq(t, 0), s.is_positive is True, bstar.is_positive is True, phi.is_positive is True, c1, c2, c10)] # 18
1187 pr(18)
1188 conds += [And(Eq(s, 0), t.is_positive is True, bstar.is_positive is True, phi.is_negative is True, c1, c3, c10)] # 19
1189 pr(19)
1190 conds += [And(Eq(n, 0), m.is_positive is True, cstar.is_positive is True, phi.is_negative is True, c1, c2, c12)] # 20
1191 pr(20)
1192 conds += [And(Eq(m, 0), n.is_positive is True, cstar.is_positive is True, phi.is_positive is True, c1, c3, c12)] # 21
1193 pr(21)
1194 conds += [And(Eq(s*t, 0), bstar.is_positive is True, cstar.is_positive is True,
1195 c1, c2, c3, c10, c12)] # 22
1196 pr(22)
1197 conds += [And(Eq(m*n, 0), bstar.is_positive is True, cstar.is_positive is True,
1198 c1, c2, c3, c10, c12)] # 23
1199 pr(23)
1201 # The following case is from [Luke1969]. As far as I can tell, it is *not*
1202 # covered by Prudnikov's.
1203 # Let G1 and G2 be the two G-functions. Suppose the integral exists from
1204 # 0 to a > 0 (this is easy the easy part), that G1 is exponential decay at
1205 # infinity, and that the mellin transform of G2 exists.
1206 # Then the integral exists.
1207 mt1_exists = _check_antecedents_1(g1, x, helper=True)
1208 mt2_exists = _check_antecedents_1(g2, x, helper=True)
1209 conds += [And(mt2_exists, Eq(t, 0), u < s, bstar.is_positive is True, c10, c1, c2, c3)]
1210 pr('E1')
1211 conds += [And(mt2_exists, Eq(s, 0), v < t, bstar.is_positive is True, c10, c1, c2, c3)]
1212 pr('E2')
1213 conds += [And(mt1_exists, Eq(n, 0), p < m, cstar.is_positive is True, c12, c1, c2, c3)]
1214 pr('E3')
1215 conds += [And(mt1_exists, Eq(m, 0), q < n, cstar.is_positive is True, c12, c1, c2, c3)]
1216 pr('E4')
1218 # Let's short-circuit if this worked ...
1219 # the rest is corner-cases and terrible to read.
1220 r = Or(*conds)
1221 if _eval_cond(r) != False:
1222 return r
1224 conds += [And(m + n > p, Eq(t, 0), Eq(phi, 0), s.is_positive is True, bstar.is_positive is True, cstar.is_negative is True,
1225 Abs(unbranched_argument(omega)) < (m + n - p + 1)*pi,
1226 c1, c2, c10, c14, c15)] # 24
1227 pr(24)
1228 conds += [And(m + n > q, Eq(s, 0), Eq(phi, 0), t.is_positive is True, bstar.is_positive is True, cstar.is_negative is True,
1229 Abs(unbranched_argument(omega)) < (m + n - q + 1)*pi,
1230 c1, c3, c10, c14, c15)] # 25
1231 pr(25)
1232 conds += [And(Eq(p, q - 1), Eq(t, 0), Eq(phi, 0), s.is_positive is True, bstar.is_positive is True,
1233 cstar >= 0, cstar*pi < Abs(unbranched_argument(omega)),
1234 c1, c2, c10, c14, c15)] # 26
1235 pr(26)
1236 conds += [And(Eq(p, q + 1), Eq(s, 0), Eq(phi, 0), t.is_positive is True, bstar.is_positive is True,
1237 cstar >= 0, cstar*pi < Abs(unbranched_argument(omega)),
1238 c1, c3, c10, c14, c15)] # 27
1239 pr(27)
1240 conds += [And(p < q - 1, Eq(t, 0), Eq(phi, 0), s.is_positive is True, bstar.is_positive is True,
1241 cstar >= 0, cstar*pi < Abs(unbranched_argument(omega)),
1242 Abs(unbranched_argument(omega)) < (m + n - p + 1)*pi,
1243 c1, c2, c10, c14, c15)] # 28
1244 pr(28)
1245 conds += [And(
1246 p > q + 1, Eq(s, 0), Eq(phi, 0), t.is_positive is True, bstar.is_positive is True, cstar >= 0,
1247 cstar*pi < Abs(unbranched_argument(omega)),
1248 Abs(unbranched_argument(omega)) < (m + n - q + 1)*pi,
1249 c1, c3, c10, c14, c15)] # 29
1250 pr(29)
1251 conds += [And(Eq(n, 0), Eq(phi, 0), s + t > 0, m.is_positive is True, cstar.is_positive is True, bstar.is_negative is True,
1252 Abs(unbranched_argument(sigma)) < (s + t - u + 1)*pi,
1253 c1, c2, c12, c14, c15)] # 30
1254 pr(30)
1255 conds += [And(Eq(m, 0), Eq(phi, 0), s + t > v, n.is_positive is True, cstar.is_positive is True, bstar.is_negative is True,
1256 Abs(unbranched_argument(sigma)) < (s + t - v + 1)*pi,
1257 c1, c3, c12, c14, c15)] # 31
1258 pr(31)
1259 conds += [And(Eq(n, 0), Eq(phi, 0), Eq(u, v - 1), m.is_positive is True, cstar.is_positive is True,
1260 bstar >= 0, bstar*pi < Abs(unbranched_argument(sigma)),
1261 Abs(unbranched_argument(sigma)) < (bstar + 1)*pi,
1262 c1, c2, c12, c14, c15)] # 32
1263 pr(32)
1264 conds += [And(Eq(m, 0), Eq(phi, 0), Eq(u, v + 1), n.is_positive is True, cstar.is_positive is True,
1265 bstar >= 0, bstar*pi < Abs(unbranched_argument(sigma)),
1266 Abs(unbranched_argument(sigma)) < (bstar + 1)*pi,
1267 c1, c3, c12, c14, c15)] # 33
1268 pr(33)
1269 conds += [And(
1270 Eq(n, 0), Eq(phi, 0), u < v - 1, m.is_positive is True, cstar.is_positive is True, bstar >= 0,
1271 bstar*pi < Abs(unbranched_argument(sigma)),
1272 Abs(unbranched_argument(sigma)) < (s + t - u + 1)*pi,
1273 c1, c2, c12, c14, c15)] # 34
1274 pr(34)
1275 conds += [And(
1276 Eq(m, 0), Eq(phi, 0), u > v + 1, n.is_positive is True, cstar.is_positive is True, bstar >= 0,
1277 bstar*pi < Abs(unbranched_argument(sigma)),
1278 Abs(unbranched_argument(sigma)) < (s + t - v + 1)*pi,
1279 c1, c3, c12, c14, c15)] # 35
1280 pr(35)
1282 return Or(*conds)
1284 # NOTE An alternative, but as far as I can tell weaker, set of conditions
1285 # can be found in [L, section 5.6.2].
1288def _int0oo(g1, g2, x):
1289 """
1290 Express integral from zero to infinity g1*g2 using a G function,
1291 assuming the necessary conditions are fulfilled.
1293 Examples
1294 ========
1296 >>> from sympy.integrals.meijerint import _int0oo
1297 >>> from sympy.abc import s, t, m
1298 >>> from sympy import meijerg, S
1299 >>> g1 = meijerg([], [], [-S(1)/2, 0], [], s**2*t/4)
1300 >>> g2 = meijerg([], [], [m/2], [-m/2], t/4)
1301 >>> _int0oo(g1, g2, t)
1302 4*meijerg(((1/2, 0), ()), ((m/2,), (-m/2,)), s**(-2))/s**2
1303 """
1304 # See: [L, section 5.6.2, equation (1)]
1305 eta, _ = _get_coeff_exp(g1.argument, x)
1306 omega, _ = _get_coeff_exp(g2.argument, x)
1308 def neg(l):
1309 return [-x for x in l]
1310 a1 = neg(g1.bm) + list(g2.an)
1311 a2 = list(g2.aother) + neg(g1.bother)
1312 b1 = neg(g1.an) + list(g2.bm)
1313 b2 = list(g2.bother) + neg(g1.aother)
1314 return meijerg(a1, a2, b1, b2, omega/eta)/eta
1317def _rewrite_inversion(fac, po, g, x):
1318 """ Absorb ``po`` == x**s into g. """
1319 _, s = _get_coeff_exp(po, x)
1320 a, b = _get_coeff_exp(g.argument, x)
1322 def tr(l):
1323 return [t + s/b for t in l]
1324 from sympy.simplify import powdenest
1325 return (powdenest(fac/a**(s/b), polar=True),
1326 meijerg(tr(g.an), tr(g.aother), tr(g.bm), tr(g.bother), g.argument))
1329def _check_antecedents_inversion(g, x):
1330 """ Check antecedents for the laplace inversion integral. """
1331 _debug('Checking antecedents for inversion:')
1332 z = g.argument
1333 _, e = _get_coeff_exp(z, x)
1334 if e < 0:
1335 _debug(' Flipping G.')
1336 # We want to assume that argument gets large as |x| -> oo
1337 return _check_antecedents_inversion(_flip_g(g), x)
1339 def statement_half(a, b, c, z, plus):
1340 coeff, exponent = _get_coeff_exp(z, x)
1341 a *= exponent
1342 b *= coeff**c
1343 c *= exponent
1344 conds = []
1345 wp = b*exp(S.ImaginaryUnit*re(c)*pi/2)
1346 wm = b*exp(-S.ImaginaryUnit*re(c)*pi/2)
1347 if plus:
1348 w = wp
1349 else:
1350 w = wm
1351 conds += [And(Or(Eq(b, 0), re(c) <= 0), re(a) <= -1)]
1352 conds += [And(Ne(b, 0), Eq(im(c), 0), re(c) > 0, re(w) < 0)]
1353 conds += [And(Ne(b, 0), Eq(im(c), 0), re(c) > 0, re(w) <= 0,
1354 re(a) <= -1)]
1355 return Or(*conds)
1357 def statement(a, b, c, z):
1358 """ Provide a convergence statement for z**a * exp(b*z**c),
1359 c/f sphinx docs. """
1360 return And(statement_half(a, b, c, z, True),
1361 statement_half(a, b, c, z, False))
1363 # Notations from [L], section 5.7-10
1364 m, n, p, q = S([len(g.bm), len(g.an), len(g.ap), len(g.bq)])
1365 tau = m + n - p
1366 nu = q - m - n
1367 rho = (tau - nu)/2
1368 sigma = q - p
1369 if sigma == 1:
1370 epsilon = S.Half
1371 elif sigma > 1:
1372 epsilon = 1
1373 else:
1374 epsilon = S.NaN
1375 theta = ((1 - sigma)/2 + Add(*g.bq) - Add(*g.ap))/sigma
1376 delta = g.delta
1377 _debugf(' m=%s, n=%s, p=%s, q=%s, tau=%s, nu=%s, rho=%s, sigma=%s',
1378 (m, n, p, q, tau, nu, rho, sigma))
1379 _debugf(' epsilon=%s, theta=%s, delta=%s', (epsilon, theta, delta))
1381 # First check if the computation is valid.
1382 if not (g.delta >= e/2 or (p >= 1 and p >= q)):
1383 _debug(' Computation not valid for these parameters.')
1384 return False
1386 # Now check if the inversion integral exists.
1388 # Test "condition A"
1389 for a, b in itertools.product(g.an, g.bm):
1390 if (a - b).is_integer and a > b:
1391 _debug(' Not a valid G function.')
1392 return False
1394 # There are two cases. If p >= q, we can directly use a slater expansion
1395 # like [L], 5.2 (11). Note in particular that the asymptotics of such an
1396 # expansion even hold when some of the parameters differ by integers, i.e.
1397 # the formula itself would not be valid! (b/c G functions are cts. in their
1398 # parameters)
1399 # When p < q, we need to use the theorems of [L], 5.10.
1401 if p >= q:
1402 _debug(' Using asymptotic Slater expansion.')
1403 return And(*[statement(a - 1, 0, 0, z) for a in g.an])
1405 def E(z):
1406 return And(*[statement(a - 1, 0, 0, z) for a in g.an])
1408 def H(z):
1409 return statement(theta, -sigma, 1/sigma, z)
1411 def Hp(z):
1412 return statement_half(theta, -sigma, 1/sigma, z, True)
1414 def Hm(z):
1415 return statement_half(theta, -sigma, 1/sigma, z, False)
1417 # [L], section 5.10
1418 conds = []
1419 # Theorem 1 -- p < q from test above
1420 conds += [And(1 <= n, 1 <= m, rho*pi - delta >= pi/2, delta > 0,
1421 E(z*exp(S.ImaginaryUnit*pi*(nu + 1))))]
1422 # Theorem 2, statements (2) and (3)
1423 conds += [And(p + 1 <= m, m + 1 <= q, delta > 0, delta < pi/2, n == 0,
1424 (m - p + 1)*pi - delta >= pi/2,
1425 Hp(z*exp(S.ImaginaryUnit*pi*(q - m))),
1426 Hm(z*exp(-S.ImaginaryUnit*pi*(q - m))))]
1427 # Theorem 2, statement (5) -- p < q from test above
1428 conds += [And(m == q, n == 0, delta > 0,
1429 (sigma + epsilon)*pi - delta >= pi/2, H(z))]
1430 # Theorem 3, statements (6) and (7)
1431 conds += [And(Or(And(p <= q - 2, 1 <= tau, tau <= sigma/2),
1432 And(p + 1 <= m + n, m + n <= (p + q)/2)),
1433 delta > 0, delta < pi/2, (tau + 1)*pi - delta >= pi/2,
1434 Hp(z*exp(S.ImaginaryUnit*pi*nu)),
1435 Hm(z*exp(-S.ImaginaryUnit*pi*nu)))]
1436 # Theorem 4, statements (10) and (11) -- p < q from test above
1437 conds += [And(1 <= m, rho > 0, delta > 0, delta + rho*pi < pi/2,
1438 (tau + epsilon)*pi - delta >= pi/2,
1439 Hp(z*exp(S.ImaginaryUnit*pi*nu)),
1440 Hm(z*exp(-S.ImaginaryUnit*pi*nu)))]
1441 # Trivial case
1442 conds += [m == 0]
1444 # TODO
1445 # Theorem 5 is quite general
1446 # Theorem 6 contains special cases for q=p+1
1448 return Or(*conds)
1451def _int_inversion(g, x, t):
1452 """
1453 Compute the laplace inversion integral, assuming the formula applies.
1454 """
1455 b, a = _get_coeff_exp(g.argument, x)
1456 C, g = _inflate_fox_h(meijerg(g.an, g.aother, g.bm, g.bother, b/t**a), -a)
1457 return C/t*g
1460####################################################################
1461# Finally, the real meat.
1462####################################################################
1464_lookup_table = None
1467@cacheit
1468@timeit
1469def _rewrite_single(f, x, recursive=True):
1470 """
1471 Try to rewrite f as a sum of single G functions of the form
1472 C*x**s*G(a*x**b), where b is a rational number and C is independent of x.
1473 We guarantee that result.argument.as_coeff_mul(x) returns (a, (x**b,))
1474 or (a, ()).
1475 Returns a list of tuples (C, s, G) and a condition cond.
1476 Returns None on failure.
1477 """
1478 from .transforms import (mellin_transform, inverse_mellin_transform,
1479 IntegralTransformError, MellinTransformStripError)
1481 global _lookup_table
1482 if not _lookup_table:
1483 _lookup_table = {}
1484 _create_lookup_table(_lookup_table)
1486 if isinstance(f, meijerg):
1487 coeff, m = factor(f.argument, x).as_coeff_mul(x)
1488 if len(m) > 1:
1489 return None
1490 m = m[0]
1491 if m.is_Pow:
1492 if m.base != x or not m.exp.is_Rational:
1493 return None
1494 elif m != x:
1495 return None
1496 return [(1, 0, meijerg(f.an, f.aother, f.bm, f.bother, coeff*m))], True
1498 f_ = f
1499 f = f.subs(x, z)
1500 t = _mytype(f, z)
1501 if t in _lookup_table:
1502 l = _lookup_table[t]
1503 for formula, terms, cond, hint in l:
1504 subs = f.match(formula, old=True)
1505 if subs:
1506 subs_ = {}
1507 for fro, to in subs.items():
1508 subs_[fro] = unpolarify(polarify(to, lift=True),
1509 exponents_only=True)
1510 subs = subs_
1511 if not isinstance(hint, bool):
1512 hint = hint.subs(subs)
1513 if hint == False:
1514 continue
1515 if not isinstance(cond, (bool, BooleanAtom)):
1516 cond = unpolarify(cond.subs(subs))
1517 if _eval_cond(cond) == False:
1518 continue
1519 if not isinstance(terms, list):
1520 terms = terms(subs)
1521 res = []
1522 for fac, g in terms:
1523 r1 = _get_coeff_exp(unpolarify(fac.subs(subs).subs(z, x),
1524 exponents_only=True), x)
1525 try:
1526 g = g.subs(subs).subs(z, x)
1527 except ValueError:
1528 continue
1529 # NOTE these substitutions can in principle introduce oo,
1530 # zoo and other absurdities. It shouldn't matter,
1531 # but better be safe.
1532 if Tuple(*(r1 + (g,))).has(S.Infinity, S.ComplexInfinity, S.NegativeInfinity):
1533 continue
1534 g = meijerg(g.an, g.aother, g.bm, g.bother,
1535 unpolarify(g.argument, exponents_only=True))
1536 res.append(r1 + (g,))
1537 if res:
1538 return res, cond
1540 # try recursive mellin transform
1541 if not recursive:
1542 return None
1543 _debug('Trying recursive Mellin transform method.')
1545 def my_imt(F, s, x, strip):
1546 """ Calling simplify() all the time is slow and not helpful, since
1547 most of the time it only factors things in a way that has to be
1548 un-done anyway. But sometimes it can remove apparent poles. """
1549 # XXX should this be in inverse_mellin_transform?
1550 try:
1551 return inverse_mellin_transform(F, s, x, strip,
1552 as_meijerg=True, needeval=True)
1553 except MellinTransformStripError:
1554 from sympy.simplify import simplify
1555 return inverse_mellin_transform(
1556 simplify(cancel(expand(F))), s, x, strip,
1557 as_meijerg=True, needeval=True)
1558 f = f_
1559 s = _dummy('s', 'rewrite-single', f)
1560 # to avoid infinite recursion, we have to force the two g functions case
1562 def my_integrator(f, x):
1563 r = _meijerint_definite_4(f, x, only_double=True)
1564 if r is not None:
1565 from sympy.simplify import hyperexpand
1566 res, cond = r
1567 res = _my_unpolarify(hyperexpand(res, rewrite='nonrepsmall'))
1568 return Piecewise((res, cond),
1569 (Integral(f, (x, S.Zero, S.Infinity)), True))
1570 return Integral(f, (x, S.Zero, S.Infinity))
1571 try:
1572 F, strip, _ = mellin_transform(f, x, s, integrator=my_integrator,
1573 simplify=False, needeval=True)
1574 g = my_imt(F, s, x, strip)
1575 except IntegralTransformError:
1576 g = None
1577 if g is None:
1578 # We try to find an expression by analytic continuation.
1579 # (also if the dummy is already in the expression, there is no point in
1580 # putting in another one)
1581 a = _dummy_('a', 'rewrite-single')
1582 if a not in f.free_symbols and _is_analytic(f, x):
1583 try:
1584 F, strip, _ = mellin_transform(f.subs(x, a*x), x, s,
1585 integrator=my_integrator,
1586 needeval=True, simplify=False)
1587 g = my_imt(F, s, x, strip).subs(a, 1)
1588 except IntegralTransformError:
1589 g = None
1590 if g is None or g.has(S.Infinity, S.NaN, S.ComplexInfinity):
1591 _debug('Recursive Mellin transform failed.')
1592 return None
1593 args = Add.make_args(g)
1594 res = []
1595 for f in args:
1596 c, m = f.as_coeff_mul(x)
1597 if len(m) > 1:
1598 raise NotImplementedError('Unexpected form...')
1599 g = m[0]
1600 a, b = _get_coeff_exp(g.argument, x)
1601 res += [(c, 0, meijerg(g.an, g.aother, g.bm, g.bother,
1602 unpolarify(polarify(
1603 a, lift=True), exponents_only=True)
1604 *x**b))]
1605 _debug('Recursive Mellin transform worked:', g)
1606 return res, True
1609def _rewrite1(f, x, recursive=True):
1610 """
1611 Try to rewrite ``f`` using a (sum of) single G functions with argument a*x**b.
1612 Return fac, po, g such that f = fac*po*g, fac is independent of ``x``.
1613 and po = x**s.
1614 Here g is a result from _rewrite_single.
1615 Return None on failure.
1616 """
1617 fac, po, g = _split_mul(f, x)
1618 g = _rewrite_single(g, x, recursive)
1619 if g:
1620 return fac, po, g[0], g[1]
1623def _rewrite2(f, x):
1624 """
1625 Try to rewrite ``f`` as a product of two G functions of arguments a*x**b.
1626 Return fac, po, g1, g2 such that f = fac*po*g1*g2, where fac is
1627 independent of x and po is x**s.
1628 Here g1 and g2 are results of _rewrite_single.
1629 Returns None on failure.
1630 """
1631 fac, po, g = _split_mul(f, x)
1632 if any(_rewrite_single(expr, x, False) is None for expr in _mul_args(g)):
1633 return None
1634 l = _mul_as_two_parts(g)
1635 if not l:
1636 return None
1637 l = list(ordered(l, [
1638 lambda p: max(len(_exponents(p[0], x)), len(_exponents(p[1], x))),
1639 lambda p: max(len(_functions(p[0], x)), len(_functions(p[1], x))),
1640 lambda p: max(len(_find_splitting_points(p[0], x)),
1641 len(_find_splitting_points(p[1], x)))]))
1643 for recursive, (fac1, fac2) in itertools.product((False, True), l):
1644 g1 = _rewrite_single(fac1, x, recursive)
1645 g2 = _rewrite_single(fac2, x, recursive)
1646 if g1 and g2:
1647 cond = And(g1[1], g2[1])
1648 if cond != False:
1649 return fac, po, g1[0], g2[0], cond
1652def meijerint_indefinite(f, x):
1653 """
1654 Compute an indefinite integral of ``f`` by rewriting it as a G function.
1656 Examples
1657 ========
1659 >>> from sympy.integrals.meijerint import meijerint_indefinite
1660 >>> from sympy import sin
1661 >>> from sympy.abc import x
1662 >>> meijerint_indefinite(sin(x), x)
1663 -cos(x)
1664 """
1665 f = sympify(f)
1666 results = []
1667 for a in sorted(_find_splitting_points(f, x) | {S.Zero}, key=default_sort_key):
1668 res = _meijerint_indefinite_1(f.subs(x, x + a), x)
1669 if not res:
1670 continue
1671 res = res.subs(x, x - a)
1672 if _has(res, hyper, meijerg):
1673 results.append(res)
1674 else:
1675 return res
1676 if f.has(HyperbolicFunction):
1677 _debug('Try rewriting hyperbolics in terms of exp.')
1678 rv = meijerint_indefinite(
1679 _rewrite_hyperbolics_as_exp(f), x)
1680 if rv:
1681 if not isinstance(rv, list):
1682 from sympy.simplify.radsimp import collect
1683 return collect(factor_terms(rv), rv.atoms(exp))
1684 results.extend(rv)
1685 if results:
1686 return next(ordered(results))
1689def _meijerint_indefinite_1(f, x):
1690 """ Helper that does not attempt any substitution. """
1691 _debug('Trying to compute the indefinite integral of', f, 'wrt', x)
1692 from sympy.simplify import hyperexpand, powdenest
1694 gs = _rewrite1(f, x)
1695 if gs is None:
1696 # Note: the code that calls us will do expand() and try again
1697 return None
1699 fac, po, gl, cond = gs
1700 _debug(' could rewrite:', gs)
1701 res = S.Zero
1702 for C, s, g in gl:
1703 a, b = _get_coeff_exp(g.argument, x)
1704 _, c = _get_coeff_exp(po, x)
1705 c += s
1707 # we do a substitution t=a*x**b, get integrand fac*t**rho*g
1708 fac_ = fac * C / (b*a**((1 + c)/b))
1709 rho = (c + 1)/b - 1
1711 # we now use t**rho*G(params, t) = G(params + rho, t)
1712 # [L, page 150, equation (4)]
1713 # and integral G(params, t) dt = G(1, params+1, 0, t)
1714 # (or a similar expression with 1 and 0 exchanged ... pick the one
1715 # which yields a well-defined function)
1716 # [R, section 5]
1717 # (Note that this dummy will immediately go away again, so we
1718 # can safely pass S.One for ``expr``.)
1719 t = _dummy('t', 'meijerint-indefinite', S.One)
1721 def tr(p):
1722 return [a + rho + 1 for a in p]
1723 if any(b.is_integer and (b <= 0) == True for b in tr(g.bm)):
1724 r = -meijerg(
1725 tr(g.an), tr(g.aother) + [1], tr(g.bm) + [0], tr(g.bother), t)
1726 else:
1727 r = meijerg(
1728 tr(g.an) + [1], tr(g.aother), tr(g.bm), tr(g.bother) + [0], t)
1729 # The antiderivative is most often expected to be defined
1730 # in the neighborhood of x = 0.
1731 if b.is_extended_nonnegative and not f.subs(x, 0).has(S.NaN, S.ComplexInfinity):
1732 place = 0 # Assume we can expand at zero
1733 else:
1734 place = None
1735 r = hyperexpand(r.subs(t, a*x**b), place=place)
1737 # now substitute back
1738 # Note: we really do want the powers of x to combine.
1739 res += powdenest(fac_*r, polar=True)
1741 def _clean(res):
1742 """This multiplies out superfluous powers of x we created, and chops off
1743 constants:
1745 >> _clean(x*(exp(x)/x - 1/x) + 3)
1746 exp(x)
1748 cancel is used before mul_expand since it is possible for an
1749 expression to have an additive constant that does not become isolated
1750 with simple expansion. Such a situation was identified in issue 6369:
1752 Examples
1753 ========
1755 >>> from sympy import sqrt, cancel
1756 >>> from sympy.abc import x
1757 >>> a = sqrt(2*x + 1)
1758 >>> bad = (3*x*a**5 + 2*x - a**5 + 1)/a**2
1759 >>> bad.expand().as_independent(x)[0]
1760 0
1761 >>> cancel(bad).expand().as_independent(x)[0]
1762 1
1763 """
1764 res = expand_mul(cancel(res), deep=False)
1765 return Add._from_args(res.as_coeff_add(x)[1])
1767 res = piecewise_fold(res, evaluate=None)
1768 if res.is_Piecewise:
1769 newargs = []
1770 for e, c in res.args:
1771 e = _my_unpolarify(_clean(e))
1772 newargs += [(e, c)]
1773 res = Piecewise(*newargs, evaluate=False)
1774 else:
1775 res = _my_unpolarify(_clean(res))
1776 return Piecewise((res, _my_unpolarify(cond)), (Integral(f, x), True))
1779@timeit
1780def meijerint_definite(f, x, a, b):
1781 """
1782 Integrate ``f`` over the interval [``a``, ``b``], by rewriting it as a product
1783 of two G functions, or as a single G function.
1785 Return res, cond, where cond are convergence conditions.
1787 Examples
1788 ========
1790 >>> from sympy.integrals.meijerint import meijerint_definite
1791 >>> from sympy import exp, oo
1792 >>> from sympy.abc import x
1793 >>> meijerint_definite(exp(-x**2), x, -oo, oo)
1794 (sqrt(pi), True)
1796 This function is implemented as a succession of functions
1797 meijerint_definite, _meijerint_definite_2, _meijerint_definite_3,
1798 _meijerint_definite_4. Each function in the list calls the next one
1799 (presumably) several times. This means that calling meijerint_definite
1800 can be very costly.
1801 """
1802 # This consists of three steps:
1803 # 1) Change the integration limits to 0, oo
1804 # 2) Rewrite in terms of G functions
1805 # 3) Evaluate the integral
1806 #
1807 # There are usually several ways of doing this, and we want to try all.
1808 # This function does (1), calls _meijerint_definite_2 for step (2).
1809 _debugf('Integrating %s wrt %s from %s to %s.', (f, x, a, b))
1810 f = sympify(f)
1811 if f.has(DiracDelta):
1812 _debug('Integrand has DiracDelta terms - giving up.')
1813 return None
1815 if f.has(SingularityFunction):
1816 _debug('Integrand has Singularity Function terms - giving up.')
1817 return None
1819 f_, x_, a_, b_ = f, x, a, b
1821 # Let's use a dummy in case any of the boundaries has x.
1822 d = Dummy('x')
1823 f = f.subs(x, d)
1824 x = d
1826 if a == b:
1827 return (S.Zero, True)
1829 results = []
1830 if a is S.NegativeInfinity and b is not S.Infinity:
1831 return meijerint_definite(f.subs(x, -x), x, -b, -a)
1833 elif a is S.NegativeInfinity:
1834 # Integrating -oo to oo. We need to find a place to split the integral.
1835 _debug(' Integrating -oo to +oo.')
1836 innermost = _find_splitting_points(f, x)
1837 _debug(' Sensible splitting points:', innermost)
1838 for c in sorted(innermost, key=default_sort_key, reverse=True) + [S.Zero]:
1839 _debug(' Trying to split at', c)
1840 if not c.is_extended_real:
1841 _debug(' Non-real splitting point.')
1842 continue
1843 res1 = _meijerint_definite_2(f.subs(x, x + c), x)
1844 if res1 is None:
1845 _debug(' But could not compute first integral.')
1846 continue
1847 res2 = _meijerint_definite_2(f.subs(x, c - x), x)
1848 if res2 is None:
1849 _debug(' But could not compute second integral.')
1850 continue
1851 res1, cond1 = res1
1852 res2, cond2 = res2
1853 cond = _condsimp(And(cond1, cond2))
1854 if cond == False:
1855 _debug(' But combined condition is always false.')
1856 continue
1857 res = res1 + res2
1858 return res, cond
1860 elif a is S.Infinity:
1861 res = meijerint_definite(f, x, b, S.Infinity)
1862 return -res[0], res[1]
1864 elif (a, b) == (S.Zero, S.Infinity):
1865 # This is a common case - try it directly first.
1866 res = _meijerint_definite_2(f, x)
1867 if res:
1868 if _has(res[0], meijerg):
1869 results.append(res)
1870 else:
1871 return res
1873 else:
1874 if b is S.Infinity:
1875 for split in _find_splitting_points(f, x):
1876 if (a - split >= 0) == True:
1877 _debugf('Trying x -> x + %s', split)
1878 res = _meijerint_definite_2(f.subs(x, x + split)
1879 *Heaviside(x + split - a), x)
1880 if res:
1881 if _has(res[0], meijerg):
1882 results.append(res)
1883 else:
1884 return res
1886 f = f.subs(x, x + a)
1887 b = b - a
1888 a = 0
1889 if b is not S.Infinity:
1890 phi = exp(S.ImaginaryUnit*arg(b))
1891 b = Abs(b)
1892 f = f.subs(x, phi*x)
1893 f *= Heaviside(b - x)*phi
1894 b = S.Infinity
1896 _debug('Changed limits to', a, b)
1897 _debug('Changed function to', f)
1898 res = _meijerint_definite_2(f, x)
1899 if res:
1900 if _has(res[0], meijerg):
1901 results.append(res)
1902 else:
1903 return res
1904 if f_.has(HyperbolicFunction):
1905 _debug('Try rewriting hyperbolics in terms of exp.')
1906 rv = meijerint_definite(
1907 _rewrite_hyperbolics_as_exp(f_), x_, a_, b_)
1908 if rv:
1909 if not isinstance(rv, list):
1910 from sympy.simplify.radsimp import collect
1911 rv = (collect(factor_terms(rv[0]), rv[0].atoms(exp)),) + rv[1:]
1912 return rv
1913 results.extend(rv)
1914 if results:
1915 return next(ordered(results))
1918def _guess_expansion(f, x):
1919 """ Try to guess sensible rewritings for integrand f(x). """
1920 res = [(f, 'original integrand')]
1922 orig = res[-1][0]
1923 saw = {orig}
1924 expanded = expand_mul(orig)
1925 if expanded not in saw:
1926 res += [(expanded, 'expand_mul')]
1927 saw.add(expanded)
1929 expanded = expand(orig)
1930 if expanded not in saw:
1931 res += [(expanded, 'expand')]
1932 saw.add(expanded)
1934 if orig.has(TrigonometricFunction, HyperbolicFunction):
1935 expanded = expand_mul(expand_trig(orig))
1936 if expanded not in saw:
1937 res += [(expanded, 'expand_trig, expand_mul')]
1938 saw.add(expanded)
1940 if orig.has(cos, sin):
1941 from sympy.simplify.fu import sincos_to_sum
1942 reduced = sincos_to_sum(orig)
1943 if reduced not in saw:
1944 res += [(reduced, 'trig power reduction')]
1945 saw.add(reduced)
1947 return res
1950def _meijerint_definite_2(f, x):
1951 """
1952 Try to integrate f dx from zero to infinity.
1954 The body of this function computes various 'simplifications'
1955 f1, f2, ... of f (e.g. by calling expand_mul(), trigexpand()
1956 - see _guess_expansion) and calls _meijerint_definite_3 with each of
1957 these in succession.
1958 If _meijerint_definite_3 succeeds with any of the simplified functions,
1959 returns this result.
1960 """
1961 # This function does preparation for (2), calls
1962 # _meijerint_definite_3 for (2) and (3) combined.
1964 # use a positive dummy - we integrate from 0 to oo
1965 # XXX if a nonnegative symbol is used there will be test failures
1966 dummy = _dummy('x', 'meijerint-definite2', f, positive=True)
1967 f = f.subs(x, dummy)
1968 x = dummy
1970 if f == 0:
1971 return S.Zero, True
1973 for g, explanation in _guess_expansion(f, x):
1974 _debug('Trying', explanation)
1975 res = _meijerint_definite_3(g, x)
1976 if res:
1977 return res
1980def _meijerint_definite_3(f, x):
1981 """
1982 Try to integrate f dx from zero to infinity.
1984 This function calls _meijerint_definite_4 to try to compute the
1985 integral. If this fails, it tries using linearity.
1986 """
1987 res = _meijerint_definite_4(f, x)
1988 if res and res[1] != False:
1989 return res
1990 if f.is_Add:
1991 _debug('Expanding and evaluating all terms.')
1992 ress = [_meijerint_definite_4(g, x) for g in f.args]
1993 if all(r is not None for r in ress):
1994 conds = []
1995 res = S.Zero
1996 for r, c in ress:
1997 res += r
1998 conds += [c]
1999 c = And(*conds)
2000 if c != False:
2001 return res, c
2004def _my_unpolarify(f):
2005 return _eval_cond(unpolarify(f))
2008@timeit
2009def _meijerint_definite_4(f, x, only_double=False):
2010 """
2011 Try to integrate f dx from zero to infinity.
2013 Explanation
2014 ===========
2016 This function tries to apply the integration theorems found in literature,
2017 i.e. it tries to rewrite f as either one or a product of two G-functions.
2019 The parameter ``only_double`` is used internally in the recursive algorithm
2020 to disable trying to rewrite f as a single G-function.
2021 """
2022 from sympy.simplify import hyperexpand
2023 # This function does (2) and (3)
2024 _debug('Integrating', f)
2025 # Try single G function.
2026 if not only_double:
2027 gs = _rewrite1(f, x, recursive=False)
2028 if gs is not None:
2029 fac, po, g, cond = gs
2030 _debug('Could rewrite as single G function:', fac, po, g)
2031 res = S.Zero
2032 for C, s, f in g:
2033 if C == 0:
2034 continue
2035 C, f = _rewrite_saxena_1(fac*C, po*x**s, f, x)
2036 res += C*_int0oo_1(f, x)
2037 cond = And(cond, _check_antecedents_1(f, x))
2038 if cond == False:
2039 break
2040 cond = _my_unpolarify(cond)
2041 if cond == False:
2042 _debug('But cond is always False.')
2043 else:
2044 _debug('Result before branch substitutions is:', res)
2045 return _my_unpolarify(hyperexpand(res)), cond
2047 # Try two G functions.
2048 gs = _rewrite2(f, x)
2049 if gs is not None:
2050 for full_pb in [False, True]:
2051 fac, po, g1, g2, cond = gs
2052 _debug('Could rewrite as two G functions:', fac, po, g1, g2)
2053 res = S.Zero
2054 for C1, s1, f1 in g1:
2055 for C2, s2, f2 in g2:
2056 r = _rewrite_saxena(fac*C1*C2, po*x**(s1 + s2),
2057 f1, f2, x, full_pb)
2058 if r is None:
2059 _debug('Non-rational exponents.')
2060 return
2061 C, f1_, f2_ = r
2062 _debug('Saxena subst for yielded:', C, f1_, f2_)
2063 cond = And(cond, _check_antecedents(f1_, f2_, x))
2064 if cond == False:
2065 break
2066 res += C*_int0oo(f1_, f2_, x)
2067 else:
2068 continue
2069 break
2070 cond = _my_unpolarify(cond)
2071 if cond == False:
2072 _debugf('But cond is always False (full_pb=%s).', full_pb)
2073 else:
2074 _debugf('Result before branch substitutions is: %s', (res, ))
2075 if only_double:
2076 return res, cond
2077 return _my_unpolarify(hyperexpand(res)), cond
2080def meijerint_inversion(f, x, t):
2081 r"""
2082 Compute the inverse laplace transform
2083 $\int_{c+i\infty}^{c-i\infty} f(x) e^{tx}\, dx$,
2084 for real c larger than the real part of all singularities of ``f``.
2086 Note that ``t`` is always assumed real and positive.
2088 Return None if the integral does not exist or could not be evaluated.
2090 Examples
2091 ========
2093 >>> from sympy.abc import x, t
2094 >>> from sympy.integrals.meijerint import meijerint_inversion
2095 >>> meijerint_inversion(1/x, x, t)
2096 Heaviside(t)
2097 """
2098 f_ = f
2099 t_ = t
2100 t = Dummy('t', polar=True) # We don't want sqrt(t**2) = abs(t) etc
2101 f = f.subs(t_, t)
2102 _debug('Laplace-inverting', f)
2103 if not _is_analytic(f, x):
2104 _debug('But expression is not analytic.')
2105 return None
2106 # Exponentials correspond to shifts; we filter them out and then
2107 # shift the result later. If we are given an Add this will not
2108 # work, but the calling code will take care of that.
2109 shift = S.Zero
2111 if f.is_Mul:
2112 args = list(f.args)
2113 elif isinstance(f, exp):
2114 args = [f]
2115 else:
2116 args = None
2118 if args:
2119 newargs = []
2120 exponentials = []
2121 while args:
2122 arg = args.pop()
2123 if isinstance(arg, exp):
2124 arg2 = expand(arg)
2125 if arg2.is_Mul:
2126 args += arg2.args
2127 continue
2128 try:
2129 a, b = _get_coeff_exp(arg.args[0], x)
2130 except _CoeffExpValueError:
2131 b = 0
2132 if b == 1:
2133 exponentials.append(a)
2134 else:
2135 newargs.append(arg)
2136 elif arg.is_Pow:
2137 arg2 = expand(arg)
2138 if arg2.is_Mul:
2139 args += arg2.args
2140 continue
2141 if x not in arg.base.free_symbols:
2142 try:
2143 a, b = _get_coeff_exp(arg.exp, x)
2144 except _CoeffExpValueError:
2145 b = 0
2146 if b == 1:
2147 exponentials.append(a*log(arg.base))
2148 newargs.append(arg)
2149 else:
2150 newargs.append(arg)
2151 shift = Add(*exponentials)
2152 f = Mul(*newargs)
2154 if x not in f.free_symbols:
2155 _debug('Expression consists of constant and exp shift:', f, shift)
2156 cond = Eq(im(shift), 0)
2157 if cond == False:
2158 _debug('but shift is nonreal, cannot be a Laplace transform')
2159 return None
2160 res = f*DiracDelta(t + shift)
2161 _debug('Result is a delta function, possibly conditional:', res, cond)
2162 # cond is True or Eq
2163 return Piecewise((res.subs(t, t_), cond))
2165 gs = _rewrite1(f, x)
2166 if gs is not None:
2167 fac, po, g, cond = gs
2168 _debug('Could rewrite as single G function:', fac, po, g)
2169 res = S.Zero
2170 for C, s, f in g:
2171 C, f = _rewrite_inversion(fac*C, po*x**s, f, x)
2172 res += C*_int_inversion(f, x, t)
2173 cond = And(cond, _check_antecedents_inversion(f, x))
2174 if cond == False:
2175 break
2176 cond = _my_unpolarify(cond)
2177 if cond == False:
2178 _debug('But cond is always False.')
2179 else:
2180 _debug('Result before branch substitution:', res)
2181 from sympy.simplify import hyperexpand
2182 res = _my_unpolarify(hyperexpand(res))
2183 if not res.has(Heaviside):
2184 res *= Heaviside(t)
2185 res = res.subs(t, t + shift)
2186 if not isinstance(cond, bool):
2187 cond = cond.subs(t, t + shift)
2188 from .transforms import InverseLaplaceTransform
2189 return Piecewise((res.subs(t, t_), cond),
2190 (InverseLaplaceTransform(f_.subs(t, t_), x, t_, None), True))