Computing Integrals using Meijer G-Functions

This text aims do describe in some detail the steps (and subtleties) involved in using Meijer G-functions for computing definite and indefinite integrals. We shall ignore proofs completely.

Overview

The algorithm to compute \(\int f(x) \mathrm{d}x\) or \(\int_0^\infty f(x) \mathrm{d}x\) generally consists of three steps:

1. Rewrite the integrand using Meijer G-functions (one or sometimes two).

2. Apply an integration theorem, to get the answer (usually expressed as another G-function).

3. Expand the result in named special functions.

Step (3) is implemented in the function hyperexpand (q.v.). Steps (1) and (2) are described below. Moreover, G-functions are usually branched. Thus our treatment of branched functions is described first.

Some other integrals (e.g. \(\int_{-\infty}^\infty\)) can also be computed by first recasting them into one of the above forms. There is a lot of choice involved here, and the algorithm is heuristic at best.

Polar Numbers and Branched Functions

Both Meijer G-Functions and Hypergeometric functions are typically branched (possible branchpoints being \(0\), \(\pm 1\), \(\infty\)). This is not very important when e.g. expanding a single hypergeometric function into named special functions, since sorting out the branches can be left to the human user. However this algorithm manipulates and transforms G-functions, and to do this correctly it needs at least some crude understanding of the branchings involved.

To begin, we consider the set \(\mathcal{S} = \{(r, \theta) : r > 0, \theta \in \mathbb{R}\}\). We have a map \(p: \mathcal{S}: \rightarrow \mathbb{C}-\{0\}, (r, \theta) \mapsto r e^{i \theta}\). Decreeing this to be a local biholomorphism gives \(\mathcal{S}\) both a topology and a complex structure. This Riemann Surface is usually referred to as the Riemann Surface of the logarithm, for the following reason: We can define maps \(\operatorname{Exp}: \mathbb{C} \rightarrow \mathcal{S}, (x + i y) \mapsto (\exp(x), y)\) and \(\operatorname{Log}: \mathcal{S} \rightarrow \mathbb{C}, (e^x, y) \mapsto x + iy\). These can both be shown to be holomorphic, and are indeed mutual inverses.

We also sometimes formally attach a point “zero” (\(0\)) to \(\mathcal{S}\) and denote the resulting object \(\mathcal{S}_0\). Notably there is no complex structure defined near \(0\). A fundamental system of neighbourhoods is given by \(\{\operatorname{Exp}(z) : \Re(z) < k\}\), which at least defines a topology. Elements of \(\mathcal{S}_0\) shall be called polar numbers. We further define functions \(\operatorname{Arg}: \mathcal{S} \rightarrow \mathbb{R}, (r, \theta) \mapsto \theta\) and \(|.|: \mathcal{S}_0 \rightarrow \mathbb{R}_{>0}, (r, \theta) \mapsto r\). These have evident meaning and are both continuous everywhere.

Using these maps many operations can be extended from \(\mathbb{C}\) to \(\mathcal{S}\). We define \(\operatorname{Exp}(a) \operatorname{Exp}(b) = \operatorname{Exp}(a + b)\) for \(a, b \in \mathbb{C}\), also for \(a \in \mathcal{S}\) and \(b \in \mathbb{C}\) we define \(a^b = \operatorname{Exp}(b \operatorname{Log}(a))\). It can be checked easily that using these definitions, many algebraic properties holding for positive reals (e.g. \((ab)^c = a^c b^c\)) which hold in \(\mathbb{C}\) only for some numbers (because of branch cuts) hold indeed for all polar numbers.

As one peculiarity it should be mentioned that addition of polar numbers is not usually defined. However, formal sums of polar numbers can be used to express branching behaviour. For example, consider the functions \(F(z) = \sqrt{1 + z}\) and \(G(a, b) = \sqrt{a + b}\), where \(a, b, z\) are polar numbers. The general rule is that functions of a single polar variable are defined in such a way that they are continuous on circles, and agree with the usual definition for positive reals. Thus if \(S(z)\) denotes the standard branch of the square root function on \(\mathbb{C}\), we are forced to define

\[\begin{split}F(z) = \begin{cases} S(p(z)) &: |z| < 1 \\ S(p(z)) &: -\pi < \operatorname{Arg}(z) + 4\pi n \le \pi \text{ for some } n \in \mathbb{Z} \\ -S(p(z)) &: \text{else} \end{cases}.\end{split}\]

(We are omitting \(|z| = 1\) here, this does not matter for integration.) Finally we define \(G(a, b) = \sqrt{a}F(b/a)\).

Representing Branched Functions on the Argand Plane

Suppose \(f: \mathcal{S} \to \mathbb{C}\) is a holomorphic function. We wish to define a function \(F\) on (part of) the complex numbers \(\mathbb{C}\) that represents \(f\) as closely as possible. This process is knows as “introducing branch cuts”. In our situation, there is actually a canonical way of doing this (which is adhered to in all of SymPy), as follows: Introduce the “cut complex plane” \(C = \mathbb{C} \setminus \mathbb{R}_{\le 0}\). Define a function \(l: C \to \mathcal{S}\) via \(re^{i\theta} \mapsto r \operatorname{Exp}(i\theta)\). Here \(r > 0\) and \(-\pi < \theta \le \pi\). Then \(l\) is holomorphic, and we define \(G = f \circ l\). This called “lifting to the principal branch” throughout the SymPy documentation.

Table Lookups and Inverse Mellin Transforms

Suppose we are given an integrand \(f(x)\) and are trying to rewrite it as a single G-function. To do this, we first split \(f(x)\) into the form \(x^s g(x)\) (where \(g(x)\) is supposed to be simpler than \(f(x)\)). This is because multiplicative powers can be absorbed into the G-function later. This splitting is done by _split_mul(f, x). Then we assemble a tuple of functions that occur in \(f\) (e.g. if \(f(x) = e^x \cos{x}\), we would assemble the tuple \((\cos, \exp)\)). This is done by the function _mytype(f, x). Next we index a lookup table (created using _create_lookup_table()) with this tuple. This (hopefully) yields a list of Meijer G-function formulae involving these functions, we then pattern-match all of them. If one fits, we were successful, otherwise not and we have to try something else.

Suppose now we want to rewrite as a product of two G-functions. To do this, we (try to) find all inequivalent ways of splitting \(f(x)\) into a product \(f_1(x) f_2(x)\). We could try these splittings in any order, but it is often a good idea to minimize (a) the number of powers occurring in \(f_i(x)\) and (b) the number of different functions occurring in \(f_i(x)\). Thus given e.g. \(f(x) = \sin{x}\, e^{x} \sin{2x}\) we should try \(f_1(x) = \sin{x}\, \sin{2x}\), \(f_2(x) = e^{x}\) first. All of this is done by the function _mul_as_two_parts(f).

Finally, we can try a recursive Mellin transform technique. Since the Meijer G-function is defined essentially as a certain inverse mellin transform, if we want to write a function \(f(x)\) as a G-function, we can compute its mellin transform \(F(s)\). If \(F(s)\) is in the right form, the G-function expression can be read off. This technique generalises many standard rewritings, e.g. \(e^{ax} e^{bx} = e^{(a + b) x}\).

One twist is that some functions don’t have mellin transforms, even though they can be written as G-functions. This is true for example for \(f(x) = e^x \sin{x}\) (the function grows too rapidly to have a mellin transform). However if the function is recognised to be analytic, then we can try to compute the mellin-transform of \(f(ax)\) for a parameter \(a\), and deduce the G-function expression by analytic continuation. (Checking for analyticity is easy. Since we can only deal with a certain subset of functions anyway, we only have to filter out those which are not analyitc.)

The function _rewrite_single does the table lookup and recursive mellin transform. The functions _rewrite1 and _rewrite2 respectively use above-mentioned helpers and _rewrite_single to rewrite their argument as respectively one or two G-functions.

Applying the Integral Theorems

If the integrand has been recast into G-functions, evaluating the integral is relatively easy. We first do some substitutions to reduce e.g. the exponent of the argument of the G-function to unity (see _rewrite_saxena_1 and _rewrite_saxena, respectively, for one or two G-functions). Next we go through a list of conditions under which the integral theorem applies. It can fail for basically two reasons: either the integral does not exist, or the manipulations in deriving the theorem may not be allowed (for more details, see this [BlogPost]).

Sometimes this can be remedied by reducing the argument of the G-functions involved. For example it is clear that the G-function representing \(e^z\) is satisfies \(G(\operatorname{Exp}(2 \pi i)z) = G(z)\) for all \(z \in \mathcal{S}\). The function meijerg.get_period() can be used to discover this, and the function principal_branch(z, period) in functions/elementary/complexes.py can be used to exploit the information. This is done transparently by the integration code.

[BlogPost]

https://nessgrh.wordpress.com/2011/07/07/tricky-branch-cuts/

The G-Function Integration Theorems

This section intends to display in detail the definite integration theorems used in the code. The following two formulae go back to Meijer (In fact he proved more general formulae; indeed in the literature formulae are usually staded in more general form. However it is very easy to deduce the general formulae from the ones we give here. It seemed best to keep the theorems as simple as possible, since they are very complicated anyway.):

1. \[\begin{split}\int_0^\infty G_{p, q}^{m, n} \left.\left(\begin{matrix} a_1, \cdots, a_p \\ b_1, \cdots, b_q \end{matrix} \right| \eta x \right) \mathrm{d}x = \frac{\prod_{j=1}^m \Gamma(b_j + 1) \prod_{j=1}^n \Gamma(-a_j)}{\eta \prod_{j=m+1}^q \Gamma(-b_j) \prod_{j=n+1}^p \Gamma(a_j + 1)}\end{split}\]
2. \[\begin{split}\int_0^\infty G_{u, v}^{s, t} \left.\left(\begin{matrix} c_1, \cdots, c_u \\ d_1, \cdots, d_v \end{matrix} \right| \sigma x \right) G_{p, q}^{m, n} \left.\left(\begin{matrix} a_1, \cdots, a_p \\ b_1, \cdots, b_q \end{matrix} \right| \omega x \right) \mathrm{d}x = G_{v+p, u+q}^{m+t, n+s} \left.\left( \begin{matrix} a_1, \cdots, a_n, -d_1, \cdots, -d_v, a_{n+1}, \cdots, a_p \\ b_1, \cdots, b_m, -c_1, \cdots, -c_u, b_{m+1}, \cdots, b_q \end{matrix} \right| \frac{\omega}{\sigma} \right)\end{split}\]

The more interesting question is under what conditions these formulae are valid. Below we detail the conditions implemented in SymPy. They are an amalgamation of conditions found in [Prudnikov1990] and [Luke1969]; please let us know if you find any errors.

Conditions of Convergence for Integral (1)

We can without loss of generality assume \(p \le q\), since the G-functions of indices \(m, n, p, q\) and of indices \(n, m, q, p\) can be related easily (see e.g. [Luke1969], section 5.3). We introduce the following notation:

\[\begin{split}\xi = m + n - p \\ \delta = m + n - \frac{p + q}{2}\end{split}\]

\[\begin{split}C_3: -\Re(b_j) < 1 \text{ for } j=1, \ldots, m \\ 0 < -\Re(a_j) \text{ for } j=1, \ldots, n\end{split}\]

\[\begin{split}C_3^*: -\Re(b_j) < 1 \text{ for } j=1, \ldots, q \\ 0 < -\Re(a_j) \text{ for } j=1, \ldots, p\end{split}\]

\[C_4: -\Re(\delta) + \frac{q + 1 - p}{2} > q - p\]

The convergence conditions will be detailed in several “cases”, numbered one to five. For later use it will be helpful to separate conditions “at infinity” from conditions “at zero”. By conditions “at infinity” we mean conditions that only depend on the behaviour of the integrand for large, positive values of \(x\), whereas by conditions “at zero” we mean conditions that only depend on the behaviour of the integrand on \((0, \epsilon)\) for any \(\epsilon > 0\). Since all our conditions are specified in terms of parameters of the G-functions, this distinction is not immediately visible. They are, however, of very distinct character mathematically; the conditions at infinity being in particular much harder to control.

In order for the integral theorem to be valid, conditions \(n\) “at zero” and “at infinity” both have to be fulfilled, for some \(n\).

These are the conditions “at infinity”:

1. \[\delta > 0 \wedge |\arg(\eta)| < \delta \pi \wedge (A \vee B \vee C),\]

   where

   \[A = 1 \le n \wedge p < q \wedge 1 \le m\]
   \[B = 1 \le p \wedge 1 \le m \wedge q = p+1 \wedge \neg (n = 0 \wedge m = p + 1 )\]
   \[C = 1 \le n \wedge q = p \wedge |\arg(\eta)| \ne (\delta - 2k)\pi \text{ for } k = 0, 1, \ldots \left\lceil \frac{\delta}{2} \right\rceil.\]
2. \[n = 0 \wedge p + 1 \le m \wedge |\arg(\eta)| < \delta \pi\]
3. \[(p < q \wedge 1 \le m \wedge \delta > 0 \wedge |\arg(\eta)| = \delta \pi) \vee (p \le q - 2 \wedge \delta = 0 \wedge \arg(\eta) = 0)\]
4. \[p = q \wedge \delta = 0 \wedge \arg(\eta) = 0 \wedge \eta \ne 0 \wedge \Re\left(\sum_{j=1}^p b_j - a_j \right) < 0\]
5. \[\delta > 0 \wedge |\arg(\eta)| < \delta \pi\]

And these are the conditions “at zero”:

1. \[\eta \ne 0 \wedge C_3\]
2. \[C_3\]
3. \[C_3 \wedge C_4\]
4. \[C_3\]
5. \[C_3\]

Conditions of Convergence for Integral (2)

We introduce the following notation:

\[b^* = s + t - \frac{u + v}{2}\]

\[c^* = m + n - \frac{p + q}{2}\]

\[\rho = \sum_{j=1}^v d_j - \sum_{j=1}^u c_j + \frac{u - v}{2} + 1\]

\[\mu = \sum_{j=1}^q b_j - \sum_{j=1}^p a_j + \frac{p - q}{2} + 1\]

\[\phi = q - p - \frac{u - v}{2} + 1\]

\[\eta = 1 - (v - u) - \mu - \rho\]

\[\psi = \frac{\pi(q - m - n) + |\arg(\omega)|}{q - p}\]

\[\theta = \frac{\pi(v - s - t) + |\arg(\sigma)|)}{v - u}\]

\[\lambda_c = (q - p)|\omega|^{1/(q - p)} \cos{\psi} + (v - u)|\sigma|^{1/(v - u)} \cos{\theta}\]

\[\lambda_{s0}(c_1, c_2) = c_1 (q - p)|\omega|^{1/(q - p)} \sin{\psi} + c_2 (v - u)|\sigma|^{1/(v - u)} \sin{\theta}\]

\[\begin{split}\lambda_s = \begin{cases} \lambda_{s0}\left(-1,-1\right) \lambda_{s0}\left(1,1\right) & \text{for}\: \arg(\omega) = 0 \wedge \arg(\sigma) = 0 \\\lambda_{s0}\left(\operatorname{sign}\left(\operatorname{\arg}\left(\omega\right)\right),-1\right) \lambda_{s0}\left(\operatorname{sign}\left(\operatorname{\arg}\left(\omega\right)\right),1\right) & \text{for}\: \arg(\omega) \ne 0 \wedge \arg(\sigma) = 0 \\\lambda_{s0}\left(-1,\operatorname{sign}\left(\operatorname{\arg}\left(\sigma\right)\right)\right) \lambda_{s0}\left(1,\operatorname{sign}\left(\operatorname{\arg}\left(\sigma\right)\right)\right) & \text{for}\: \arg(\omega) = 0 \wedge \arg(\sigma) \ne 0) \\\lambda_{s0}\left(\operatorname{sign}\left(\operatorname{\arg}\left(\omega\right)\right),\operatorname{sign}\left(\operatorname{\arg}\left(\sigma\right)\right)\right) & \text{otherwise} \end{cases}\end{split}\]

\[z_0 = \frac{\omega}{\sigma} e^{-i\pi (b^* + c^*)}\]

\[z_1 = \frac{\sigma}{\omega} e^{-i\pi (b^* + c^*)}\]

The following conditions will be helpful:

\[\begin{split}C_1: (a_i - b_j \notin \mathbb{Z}_{>0} \text{ for } i = 1, \ldots, n, j = 1, \ldots, m) \\ \wedge (c_i - d_j \notin \mathbb{Z}_{>0} \text{ for } i = 1, \ldots, t, j = 1, \ldots, s)\end{split}\]

\[C_2: \Re(1 + b_i + d_j) > 0 \text{ for } i = 1, \ldots, m, j = 1, \ldots, s\]

\[C_3: \Re(a_i + c_j) < 1 \text{ for } i = 1, \ldots, n, j = 1, \ldots, t\]

\[C_4: (p - q)\Re(c_i) - \Re(\mu) > -\frac{3}{2} \text{ for } i=1, \ldots, t\]

\[C_5: (p - q)\Re(1 + d_i) - \Re(\mu) > -\frac{3}{2} \text{ for } i=1, \ldots, s\]

\[C_6: (u - v)\Re(a_i) - \Re(\rho) > -\frac{3}{2} \text{ for } i=1, \ldots, n\]

\[C_7: (u - v)\Re(1 + b_i) - \Re(\rho) > -\frac{3}{2} \text{ for } i=1, \ldots, m\]

\[C_8: 0 < \lvert{\phi}\rvert + 2 \Re\left(\left(\mu -1\right) \left(- u + v\right) + \left(- p + q\right) \left(\rho -1\right) + \left(- p + q\right) \left(- u + v\right)\right)\]

\[C_9: 0 < \lvert{\phi}\rvert - 2 \Re\left(\left(\mu -1\right) \left(- u + v\right) + \left(- p + q\right) \left(\rho -1\right) + \left(- p + q\right) \left(- u + v\right)\right)\]

\[C_{10}: \lvert{\operatorname{arg}\left(\sigma\right)}\rvert < \pi b^{*}\]

\[C_{11}: \lvert{\operatorname{arg}\left(\sigma\right)}\rvert = \pi b^{*}\]

\[C_{12}: |\arg(\omega)| < c^*\pi\]

\[C_{13}: |\arg(\omega)| = c^*\pi\]

\[C_{14}^1: \left(z_0 \ne 1 \wedge |\arg(1 - z_0)| < \pi \right) \vee \left(z_0 = 1 \wedge \Re(\mu + \rho - u + v) < 1 \right)\]

\[C_{14}^2: \left(z_1 \ne 1 \wedge |\arg(1 - z_1)| < \pi \right) \vee \left(z_1 = 1 \wedge \Re(\mu + \rho - p + q) < 1 \right)\]

\[C_{14}: \phi = 0 \wedge b^* + c^* \le 1 \wedge (C_{14}^1 \vee C_{14}^2)\]

\[C_{15}: \lambda_c > 0 \vee (\lambda_c = 0 \wedge \lambda_s \ne 0 \wedge \Re(\eta) > -1) \vee (\lambda_c = 0 \wedge \lambda_s = 0 \wedge \Re(\eta) > 0)\]

\[C_{16}: \int_0^\infty G_{u, v}^{s, t}(\sigma x) \mathrm{d} x \text{ converges at infinity }\]

\[C_{17}: \int_0^\infty G_{p, q}^{m, n}(\omega x) \mathrm{d} x \text{ converges at infinity }\]

Note that \(C_{16}\) and \(C_{17}\) are the reason we split the convergence conditions for integral (1).

With this notation established, the implemented convergence conditions can be enumerated as follows:

1. \[m n s t \neq 0 \wedge 0 < b^{*} \wedge 0 < c^{*} \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{10} \wedge C_{12}\]
2. \[u = v \wedge b^{*} = 0 \wedge 0 < c^{*} \wedge 0 < \sigma \wedge \Re{\rho} < 1 \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{12}\]
3. \[p = q \wedge u = v \wedge b^{*} = 0 \wedge c^{*} = 0 \wedge 0 < \sigma \wedge 0 < \omega \wedge \Re{\mu} < 1 \wedge \Re{\rho} < 1 \wedge \sigma \neq \omega \wedge C_{1} \wedge C_{2} \wedge C_{3}\]
4. \[p = q \wedge u = v \wedge b^{*} = 0 \wedge c^{*} = 0 \wedge 0 < \sigma \wedge 0 < \omega \wedge \Re\left(\mu + \rho\right) < 1 \wedge \omega \neq \sigma \wedge C_{1} \wedge C_{2} \wedge C_{3}\]
5. \[p = q \wedge u = v \wedge b^{*} = 0 \wedge c^{*} = 0 \wedge 0 < \sigma \wedge 0 < \omega \wedge \Re\left(\mu + \rho\right) < 1 \wedge \omega \neq \sigma \wedge C_{1} \wedge C_{2} \wedge C_{3}\]
6. \[q < p \wedge 0 < s \wedge 0 < b^{*} \wedge 0 \leq c^{*} \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{5} \wedge C_{10} \wedge C_{13}\]
7. \[p < q \wedge 0 < t \wedge 0 < b^{*} \wedge 0 \leq c^{*} \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{4} \wedge C_{10} \wedge C_{13}\]
8. \[v < u \wedge 0 < m \wedge 0 < c^{*} \wedge 0 \leq b^{*} \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{7} \wedge C_{11} \wedge C_{12}\]
9. \[u < v \wedge 0 < n \wedge 0 < c^{*} \wedge 0 \leq b^{*} \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{6} \wedge C_{11} \wedge C_{12}\]
10. \[q < p \wedge u = v \wedge b^{*} = 0 \wedge 0 \leq c^{*} \wedge 0 < \sigma \wedge \Re{\rho} < 1 \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{5} \wedge C_{13}\]
11. \[p < q \wedge u = v \wedge b^{*} = 0 \wedge 0 \leq c^{*} \wedge 0 < \sigma \wedge \Re{\rho} < 1 \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{4} \wedge C_{13}\]
12. \[p = q \wedge v < u \wedge 0 \leq b^{*} \wedge c^{*} = 0 \wedge 0 < \omega \wedge \Re{\mu} < 1 \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{7} \wedge C_{11}\]
13. \[p = q \wedge u < v \wedge 0 \leq b^{*} \wedge c^{*} = 0 \wedge 0 < \omega \wedge \Re{\mu} < 1 \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{6} \wedge C_{11}\]
14. \[p < q \wedge v < u \wedge 0 \leq b^{*} \wedge 0 \leq c^{*} \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{4} \wedge C_{7} \wedge C_{11} \wedge C_{13}\]
15. \[q < p \wedge u < v \wedge 0 \leq b^{*} \wedge 0 \leq c^{*} \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{5} \wedge C_{6} \wedge C_{11} \wedge C_{13}\]
16. \[q < p \wedge v < u \wedge 0 \leq b^{*} \wedge 0 \leq c^{*} \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{5} \wedge C_{7} \wedge C_{8} \wedge C_{11} \wedge C_{13} \wedge C_{14}\]
17. \[p < q \wedge u < v \wedge 0 \leq b^{*} \wedge 0 \leq c^{*} \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{4} \wedge C_{6} \wedge C_{9} \wedge C_{11} \wedge C_{13} \wedge C_{14}\]
18. \[t = 0 \wedge 0 < s \wedge 0 < b^{*} \wedge 0 < \phi \wedge C_{1} \wedge C_{2} \wedge C_{10}\]
19. \[s = 0 \wedge 0 < t \wedge 0 < b^{*} \wedge \phi < 0 \wedge C_{1} \wedge C_{3} \wedge C_{10}\]
20. \[n = 0 \wedge 0 < m \wedge 0 < c^{*} \wedge \phi < 0 \wedge C_{1} \wedge C_{2} \wedge C_{12}\]
21. \[m = 0 \wedge 0 < n \wedge 0 < c^{*} \wedge 0 < \phi \wedge C_{1} \wedge C_{3} \wedge C_{12}\]
22. \[s t = 0 \wedge 0 < b^{*} \wedge 0 < c^{*} \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{10} \wedge C_{12}\]
23. \[m n = 0 \wedge 0 < b^{*} \wedge 0 < c^{*} \wedge C_{1} \wedge C_{2} \wedge C_{3} \wedge C_{10} \wedge C_{12}\]
24. \[\begin{split}p < m + n \wedge t = 0 \wedge \phi = 0 \wedge 0 < s \wedge 0 < b^{*} \wedge c^{*} < 0 \wedge \lvert{\operatorname{arg}\left(\omega\right)}\rvert < \\ \pi \left(m + n - p + 1\right) \wedge C_{1} \wedge C_{2} \wedge C_{10} \wedge C_{14} \wedge C_{15}\end{split}\]
25. \[\begin{split}q < m + n \wedge s = 0 \wedge \phi = 0 \wedge 0 < t \wedge 0 < b^{*} \wedge c^{*} < 0 \wedge \lvert{\operatorname{arg}\left(\omega\right)}\rvert < \\ \pi \left(m + n - q + 1\right) \wedge C_{1} \wedge C_{3} \wedge C_{10} \wedge C_{14} \wedge C_{15}\end{split}\]
26. \[p = q -1 \wedge t = 0 \wedge \phi = 0 \wedge 0 < s \wedge 0 < b^{*} \wedge 0 \leq c^{*} \wedge \pi c^{*} < \lvert{\operatorname{arg}\left(\omega\right)}\rvert \wedge C_{1} \wedge C_{2} \wedge C_{10} \wedge C_{14} \wedge C_{15}\]
27. \[p = q + 1 \wedge s = 0 \wedge \phi = 0 \wedge 0 < t \wedge 0 < b^{*} \wedge 0 \leq c^{*} \wedge \pi c^{*} < \lvert{\operatorname{arg}\left(\omega\right)}\rvert \wedge C_{1} \wedge C_{3} \wedge C_{10} \wedge C_{14} \wedge C_{15}\]
28. \[\begin{split}p < q -1 \wedge t = 0 \wedge \phi = 0 \wedge 0 < s \wedge 0 < b^{*} \wedge 0 \leq c^{*} \wedge \pi c^{*} < \lvert{\operatorname{arg}\left(\omega\right)}\rvert \wedge \lvert{\operatorname{arg}\left(\omega\right)}\rvert < \\ \pi \left(m + n - p + 1\right) \wedge C_{1} \wedge C_{2} \wedge C_{10} \wedge C_{14} \wedge C_{15}\end{split}\]
29. \[\begin{split}q + 1 < p \wedge s = 0 \wedge \phi = 0 \wedge 0 < t \wedge 0 < b^{*} \wedge 0 \leq c^{*} \wedge \pi c^{*} < \lvert{\operatorname{arg}\left(\omega\right)}\rvert \wedge \lvert{\operatorname{arg}\left(\omega\right)}\rvert < \\ \pi \left(m + n - q + 1 \right) \wedge C_{1} \wedge C_{3} \wedge C_{10} \wedge C_{14} \wedge C_{15}\end{split}\]
30. \[\begin{split}n = 0 \wedge \phi = 0 \wedge 0 < s + t \wedge 0 < m \wedge 0 < c^{*} \wedge b^{*} < 0 \wedge \lvert{\operatorname{arg}\left(\sigma\right)}\rvert < \\ \pi \left(s + t - u + 1\right) \wedge C_{1} \wedge C_{2} \wedge C_{12} \wedge C_{14} \wedge C_{15}\end{split}\]
31. \[\begin{split}m = 0 \wedge \phi = 0 \wedge v < s + t \wedge 0 < n \wedge 0 < c^{*} \wedge b^{*} < 0 \wedge \lvert{\operatorname{arg}\left(\sigma\right)}\rvert < \\ \pi \left(s + t - v + 1\right) \wedge C_{1} \wedge C_{3} \wedge C_{12} \wedge C_{14} \wedge C_{15}\end{split}\]
32. \[\begin{split}n = 0 \wedge \phi = 0 \wedge u = v -1 \wedge 0 < m \wedge 0 < c^{*} \wedge 0 \leq b^{*} \wedge \pi b^{*} < \lvert{\operatorname{arg}\left(\sigma\right)}\rvert \wedge \lvert{\operatorname{arg}\left(\sigma\right)}\rvert < \\ \pi \left(b^{*} + 1\right) \wedge C_{1} \wedge C_{2} \wedge C_{12} \wedge C_{14} \wedge C_{15}\end{split}\]
33. \[\begin{split}m = 0 \wedge \phi = 0 \wedge u = v + 1 \wedge 0 < n \wedge 0 < c^{*} \wedge 0 \leq b^{*} \wedge \pi b^{*} < \lvert{\operatorname{arg}\left(\sigma\right)}\rvert \wedge \lvert{\operatorname{arg}\left(\sigma\right)}\rvert < \\ \pi \left(b^{*} + 1\right) \wedge C_{1} \wedge C_{3} \wedge C_{12} \wedge C_{14} \wedge C_{15}\end{split}\]
34. \[\begin{split}n = 0 \wedge \phi = 0 \wedge u < v -1 \wedge 0 < m \wedge 0 < c^{*} \wedge 0 \leq b^{*} \wedge \pi b^{*} < \lvert{\operatorname{arg}\left(\sigma\right)}\rvert \wedge \lvert{\operatorname{arg}\left(\sigma\right)}\rvert < \\ \pi \left(s + t - u + 1\right) \wedge C_{1} \wedge C_{2} \wedge C_{12} \wedge C_{14} \wedge C_{15}\end{split}\]
35. \[\begin{split}m = 0 \wedge \phi = 0 \wedge v + 1 < u \wedge 0 < n \wedge 0 < c^{*} \wedge 0 \leq b^{*} \wedge \pi b^{*} < \lvert{\operatorname{arg}\left(\sigma\right)}\rvert \wedge \lvert{\operatorname{arg}\left(\sigma\right)}\rvert < \\ \pi \left(s + t - v + 1 \right) \wedge C_{1} \wedge C_{3} \wedge C_{12} \wedge C_{14} \wedge C_{15}\end{split}\]
36. \[C_{17} \wedge t = 0 \wedge u < s \wedge 0 < b^{*} \wedge C_{10} \wedge C_{1} \wedge C_{2} \wedge C_{3}\]
37. \[C_{17} \wedge s = 0 \wedge v < t \wedge 0 < b^{*} \wedge C_{10} \wedge C_{1} \wedge C_{2} \wedge C_{3}\]
38. \[C_{16} \wedge n = 0 \wedge p < m \wedge 0 < c^{*} \wedge C_{12} \wedge C_{1} \wedge C_{2} \wedge C_{3}\]
39. \[C_{16} \wedge m = 0 \wedge q < n \wedge 0 < c^{*} \wedge C_{12} \wedge C_{1} \wedge C_{2} \wedge C_{3}\]

The Inverse Laplace Transform of a G-function

The inverse laplace transform of a Meijer G-function can be expressed as another G-function. This is a fairly versatile method for computing this transform. However, I could not find the details in the literature, so I work them out here. In [Luke1969], section 5.6.3, there is a formula for the inverse Laplace transform of a G-function of argument \(bz\), and convergence conditions are also given. However, we need a formula for argument \(bz^a\) for rational \(a\).

We are asked to compute

\[f(t) = \frac{1}{2\pi i} \int_{c-i\infty}^{c+i\infty} e^{zt} G(bz^a) \mathrm{d}z,\]

for positive real \(t\). Three questions arise:

1. When does this integral converge?

2. How can we compute the integral?

3. When is our computation valid?

How to compute the integral

We shall work formally for now. Denote by \(\Delta(s)\) the product of gamma functions appearing in the definition of \(G\), so that

\[G(z) = \frac{1}{2\pi i} \int_L \Delta(s) z^s \mathrm{d}s.\]

Thus

\[f(t) = \frac{1}{(2\pi i)^2} \int_{c - i\infty}^{c + i\infty} \int_L e^{zt} \Delta(s) b^s z^{as} \mathrm{d}s \mathrm{d}z.\]

We interchange the order of integration to get

\[f(t) = \frac{1}{2\pi i} \int_L b^s \Delta(s) \int_{c-i\infty}^{c+i\infty} e^{zt} z^{as} \frac{\mathrm{d}z}{2\pi i} \mathrm{d}s.\]

The inner integral is easily seen to be \(\frac{1}{\Gamma(-as)} \frac{1}{t^{1+as}}\). (Using Cauchy’s theorem and Jordan’s lemma deform the contour to run from \(-\infty\) to \(-\infty\), encircling \(0\) once in the negative sense. For \(as\) real and greater than one, this contour can be pushed onto the negative real axis and the integral is recognised as a product of a sine and a gamma function. The formula is then proved using the functional equation of the gamma function, and extended to the entire domain of convergence of the original integral by appealing to analytic continuation.) Hence we find

\[f(t) = \frac{1}{t} \frac{1}{2\pi i} \int_L \Delta(s) \frac{1}{\Gamma(-as)} \left(\frac{b}{t^a}\right)^s \mathrm{d}s,\]

which is a so-called Fox H function (of argument \(\frac{b}{t^a}\)). For rational \(a\), this can be expressed as a Meijer G-function using the gamma function multiplication theorem.

When this computation is valid

There are a number of obstacles in this computation. Interchange of integrals is only valid if all integrals involved are absolutely convergent. In particular the inner integral has to converge. Also, for our identification of the final integral as a Fox H / Meijer G-function to be correct, the poles of the newly obtained gamma function must be separated properly.

It is easy to check that the inner integral converges absolutely for \(\Re(as) < -1\). Thus the contour \(L\) has to run left of the line \(\Re(as) = -1\). Under this condition, the poles of the newly-introduced gamma function are separated properly.

It remains to observe that the Meijer G-function is an analytic, unbranched function of its parameters, and of the coefficient \(b\). Hence so is \(f(t)\). Thus the final computation remains valid as long as the initial integral converges, and if there exists a changed set of parameters where the computation is valid. If we assume w.l.o.g. that \(a > 0\), then the latter condition is fulfilled if \(G\) converges along contours (2) or (3) of [Luke1969], section 5.2, i.e. either \(\delta \ge \frac{a}{2}\) or \(p \ge 1, p \ge q\).

When the integral exists

Using [Luke1969], section 5.10, for any given meijer G-function we can find a dominant term of the form \(z^a e^{bz^c}\) (although this expression might not be the best possible, because of cancellation).

We must thus investigate

\[\lim_{T \to \infty} \int_{c-iT}^{c+iT} e^{zt} z^a e^{bz^c} \mathrm{d}z.\]

(This principal value integral is the exact statement used in the Laplace inversion theorem.) We write \(z = c + i \tau\). Then \(arg(z) \to \pm \frac{\pi}{2}\), and so \(e^{zt} \sim e^{it \tau}\) (where \(\sim\) shall always mean “asymptotically equivalent up to a positive real multiplicative constant”). Also \(z^{x + iy} \sim |\tau|^x e^{i y \log{|\tau|}} e^{\pm x i \frac{\pi}{2}}.\)

Set \(\omega_{\pm} = b e^{\pm i \Re(c) \frac{\pi}{2}}\). We have three cases:

1. \(b=0\) or \(\Re(c) \le 0\). In this case the integral converges if \(\Re(a) \le -1\).

2. \(b \ne 0\), \(\Im(c) = 0\), \(\Re(c) > 0\). In this case the integral converges if \(\Re(\omega_{\pm}) < 0\).

3. \(b \ne 0\), \(\Im(c) = 0\), \(\Re(c) > 0\), \(\Re(\omega_{\pm}) \le 0\), and at least one of \(\Re(\omega_{\pm}) = 0\). Here the same condition as in (1) applies.

Implemented G-Function Formulae

An important part of the algorithm is a table expressing various functions as Meijer G-functions. This is essentially a table of Mellin Transforms in disguise. The following automatically generated table shows the formulae currently implemented in SymPy. An entry “generated” means that the corresponding G-function has a variable number of parameters. This table is intended to shrink in future, when the algorithm’s capabilities of deriving new formulae improve. Of course it has to grow whenever a new class of special functions is to be dealt with.

Elementary functions:

\[\begin{split}a = a {G_{1, 1}^{1, 0}\left(\begin{matrix} & 1 \\0 & \end{matrix} \middle| {z} \right)} + a {G_{1, 1}^{0, 1}\left(\begin{matrix} 1 & \\ & 0 \end{matrix} \middle| {z} \right)}\end{split}\]

\[\begin{split}\left(z^{q} p + b\right)^{- a} = \frac{b^{- a} {G_{1, 1}^{1, 1}\left(\begin{matrix} 1 - a & \\0 & \end{matrix} \middle| {\frac{z^{q} p}{b}} \right)}}{\Gamma\left(a\right)}\end{split}\]

\[\begin{split}\frac{- b^{a} + \left(z^{q} p\right)^{a}}{z^{q} p - b} = \frac{b^{a - 1} {G_{2, 2}^{2, 2}\left(\begin{matrix} 0, a & \\0, a & \end{matrix} \middle| {\frac{z^{q} p}{b}} \right)} \sin{\left(\pi a \right)}}{\pi}\end{split}\]

\[\begin{split}\left(a + \sqrt{z^{q} p + a^{2}}\right)^{b} = - \frac{a^{b} b {G_{2, 2}^{1, 2}\left(\begin{matrix} \frac{b}{2} + \frac{1}{2}, \frac{b}{2} + 1 & \\0 & b \end{matrix} \middle| {\frac{z^{q} p}{a^{2}}} \right)}}{2 \sqrt{\pi}}\end{split}\]

\[\begin{split}\left(- a + \sqrt{z^{q} p + a^{2}}\right)^{b} = \frac{a^{b} b {G_{2, 2}^{1, 2}\left(\begin{matrix} \frac{b}{2} + \frac{1}{2}, \frac{b}{2} + 1 & \\b & 0 \end{matrix} \middle| {\frac{z^{q} p}{a^{2}}} \right)}}{2 \sqrt{\pi}}\end{split}\]

\[\begin{split}\frac{\left(a + \sqrt{z^{q} p + a^{2}}\right)^{b}}{\sqrt{z^{q} p + a^{2}}} = \frac{a^{b - 1} {G_{2, 2}^{1, 2}\left(\begin{matrix} \frac{b}{2}, \frac{b}{2} + \frac{1}{2} & \\0 & b \end{matrix} \middle| {\frac{z^{q} p}{a^{2}}} \right)}}{\sqrt{\pi}}\end{split}\]

\[\begin{split}\frac{\left(- a + \sqrt{z^{q} p + a^{2}}\right)^{b}}{\sqrt{z^{q} p + a^{2}}} = \frac{a^{b - 1} {G_{2, 2}^{1, 2}\left(\begin{matrix} \frac{b}{2}, \frac{b}{2} + \frac{1}{2} & \\b & 0 \end{matrix} \middle| {\frac{z^{q} p}{a^{2}}} \right)}}{\sqrt{\pi}}\end{split}\]

\[\begin{split}\left(z^{\frac{q}{2}} \sqrt{p} + \sqrt{z^{q} p + a}\right)^{b} = - \frac{a^{\frac{b}{2}} b {G_{2, 2}^{2, 1}\left(\begin{matrix} \frac{b}{2} + 1 & 1 - \frac{b}{2} \\0, \frac{1}{2} & \end{matrix} \middle| {\frac{z^{q} p}{a}} \right)}}{2 \sqrt{\pi}}\end{split}\]

\[\begin{split}\left(- z^{\frac{q}{2}} \sqrt{p} + \sqrt{z^{q} p + a}\right)^{b} = \frac{a^{\frac{b}{2}} b {G_{2, 2}^{2, 1}\left(\begin{matrix} 1 - \frac{b}{2} & \frac{b}{2} + 1 \\0, \frac{1}{2} & \end{matrix} \middle| {\frac{z^{q} p}{a}} \right)}}{2 \sqrt{\pi}}\end{split}\]

\[\begin{split}\frac{\left(z^{\frac{q}{2}} \sqrt{p} + \sqrt{z^{q} p + a}\right)^{b}}{\sqrt{z^{q} p + a}} = \frac{a^{\frac{b}{2} - \frac{1}{2}} {G_{2, 2}^{2, 1}\left(\begin{matrix} \frac{b}{2} + \frac{1}{2} & \frac{1}{2} - \frac{b}{2} \\0, \frac{1}{2} & \end{matrix} \middle| {\frac{z^{q} p}{a}} \right)}}{\sqrt{\pi}}\end{split}\]

\[\begin{split}\frac{\left(- z^{\frac{q}{2}} \sqrt{p} + \sqrt{z^{q} p + a}\right)^{b}}{\sqrt{z^{q} p + a}} = \frac{a^{\frac{b}{2} - \frac{1}{2}} {G_{2, 2}^{2, 1}\left(\begin{matrix} \frac{1}{2} - \frac{b}{2} & \frac{b}{2} + \frac{1}{2} \\0, \frac{1}{2} & \end{matrix} \middle| {\frac{z^{q} p}{a}} \right)}}{\sqrt{\pi}}\end{split}\]

Functions involving \(\theta\left(z^{q} p - b\right)\):

\[\begin{split}\left(z^{q} p - b\right)^{a - 1} \theta\left(z^{q} p - b\right) = b^{a - 1} {G_{1, 1}^{0, 1}\left(\begin{matrix} a & \\ & 0 \end{matrix} \middle| {\frac{z^{q} p}{b}} \right)} \Gamma\left(a\right),\text{ if } b > 0\end{split}\]

\[\begin{split}\left(- z^{q} p + b\right)^{a - 1} \theta\left(- z^{q} p + b\right) = b^{a - 1} {G_{1, 1}^{1, 0}\left(\begin{matrix} & a \\0 & \end{matrix} \middle| {\frac{z^{q} p}{b}} \right)} \Gamma\left(a\right),\text{ if } b > 0\end{split}\]

\[\begin{split}\left(z^{q} p - b\right)^{a - 1} \theta\left(z - \left(\frac{b}{p}\right)^{\frac{1}{q}}\right) = b^{a - 1} {G_{1, 1}^{0, 1}\left(\begin{matrix} a & \\ & 0 \end{matrix} \middle| {\frac{z^{q} p}{b}} \right)} \Gamma\left(a\right),\text{ if } b > 0\end{split}\]

\[\begin{split}\left(- z^{q} p + b\right)^{a - 1} \theta\left(- z + \left(\frac{b}{p}\right)^{\frac{1}{q}}\right) = b^{a - 1} {G_{1, 1}^{1, 0}\left(\begin{matrix} & a \\0 & \end{matrix} \middle| {\frac{z^{q} p}{b}} \right)} \Gamma\left(a\right),\text{ if } b > 0\end{split}\]

Functions involving \(\left|{z^{q} p - b}\right|\):

\[\begin{split}\left|{z^{q} p - b}\right|^{- a} = 2 {G_{2, 2}^{1, 1}\left(\begin{matrix} 1 - a & \frac{1}{2} - \frac{a}{2} \\0 & \frac{1}{2} - \frac{a}{2} \end{matrix} \middle| {\frac{z^{q} p}{b}} \right)} \sin{\left(\frac{\pi a}{2} \right)} \left|{b}\right|^{- a} \Gamma\left(1 - a\right),\text{ if } \operatorname{re}{\left(a\right)} < 1\end{split}\]

Functions involving \(e^{z^{q} p e^{i \pi}}\):

\[\begin{split}e^{z^{q} p e^{i \pi}} = {G_{0, 1}^{1, 0}\left(\begin{matrix} & \\0 & \end{matrix} \middle| {z^{q} p} \right)}\end{split}\]

Functions involving \(\sinh{\left(z^{q} p \right)}\):

\[\begin{split}\sinh{\left(z^{q} p \right)} = \pi^{\frac{3}{2}} {G_{1, 3}^{1, 0}\left(\begin{matrix} & 1 \\\frac{1}{2} & 0, 1 \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}\end{split}\]

Functions involving \(\cosh{\left(z^{q} p \right)}\):

\[\begin{split}\cosh{\left(z^{q} p \right)} = \pi^{\frac{3}{2}} {G_{1, 3}^{1, 0}\left(\begin{matrix} & \frac{1}{2} \\0 & \frac{1}{2}, \frac{1}{2} \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}\end{split}\]

Functions involving \(\sin{\left(z^{q} p \right)}\):

\[\begin{split}\sin{\left(z^{q} p \right)} = \sqrt{\pi} {G_{0, 2}^{1, 0}\left(\begin{matrix} & \\\frac{1}{2} & 0 \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}\end{split}\]

Functions involving \(\cos{\left(z^{q} p \right)}\):

\[\begin{split}\cos{\left(z^{q} p \right)} = \sqrt{\pi} {G_{0, 2}^{1, 0}\left(\begin{matrix} & \\0 & \frac{1}{2} \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}\end{split}\]

Functions involving \(\operatorname{sinc}{\left(z^{q} p \right)}\):

\[\begin{split}\operatorname{sinc}{\left(z^{q} p \right)} = \frac{\sqrt{\pi} {G_{0, 2}^{1, 0}\left(\begin{matrix} & \\0 & - \frac{1}{2} \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}}{2}\end{split}\]

Functions involving \(\log{\left(z^{q} p \right)}\), \(\theta\left(- z^{q} p + 1\right)\):

\[\log{\left(z^{q} p \right)}^{n} \theta\left(- z^{q} p + 1\right) = \text{generated}\]

\[\log{\left(z^{q} p \right)}^{n} \theta\left(z^{q} p - 1\right) = \text{generated}\]

Functions involving \(\log{\left(z^{q} p \right)}\):

\[\log{\left(z^{q} p \right)}^{n} = \text{generated}\]

\[\begin{split}\log{\left(z^{q} p + a \right)} = {G_{1, 1}^{1, 0}\left(\begin{matrix} & 1 \\0 & \end{matrix} \middle| {z} \right)} \log{\left(a \right)} + {G_{1, 1}^{0, 1}\left(\begin{matrix} 1 & \\ & 0 \end{matrix} \middle| {z} \right)} \log{\left(a \right)} + {G_{2, 2}^{1, 2}\left(\begin{matrix} 1, 1 & \\1 & 0 \end{matrix} \middle| {\frac{z^{q} p}{a}} \right)}\end{split}\]

\[\begin{split}\log{\left(\left|{z^{q} p - a}\right| \right)} = {G_{1, 1}^{1, 0}\left(\begin{matrix} & 1 \\0 & \end{matrix} \middle| {z} \right)} \log{\left(\left|{a}\right| \right)} + {G_{1, 1}^{0, 1}\left(\begin{matrix} 1 & \\ & 0 \end{matrix} \middle| {z} \right)} \log{\left(\left|{a}\right| \right)} + \pi {G_{3, 3}^{1, 2}\left(\begin{matrix} 1, 1 & \frac{1}{2} \\1 & 0, \frac{1}{2} \end{matrix} \middle| {\frac{z^{q} p}{a}} \right)}\end{split}\]

Functions involving \(\operatorname{Ei}{\left(z^{q} p \right)}\):

\[\begin{split}\operatorname{Ei}{\left(z^{q} p \right)} = - i \pi {G_{1, 1}^{1, 0}\left(\begin{matrix} & 1 \\0 & \end{matrix} \middle| {z} \right)} - {G_{1, 2}^{2, 0}\left(\begin{matrix} & 1 \\0, 0 & \end{matrix} \middle| {z^{q} p e^{i \pi}} \right)} - i \pi {G_{1, 1}^{0, 1}\left(\begin{matrix} 1 & \\ & 0 \end{matrix} \middle| {z} \right)}\end{split}\]

Functions involving \(\operatorname{Si}{\left(z^{q} p \right)}\):

\[\begin{split}\operatorname{Si}{\left(z^{q} p \right)} = \frac{\sqrt{\pi} {G_{1, 3}^{1, 1}\left(\begin{matrix} 1 & \\\frac{1}{2} & 0, 0 \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}}{2}\end{split}\]

Functions involving \(\operatorname{Ci}{\left(z^{q} p \right)}\):

\[\begin{split}\operatorname{Ci}{\left(z^{q} p \right)} = - \frac{\sqrt{\pi} {G_{1, 3}^{2, 0}\left(\begin{matrix} & 1 \\0, 0 & \frac{1}{2} \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}}{2}\end{split}\]

Functions involving \(\operatorname{Shi}{\left(z^{q} p \right)}\):

\[\begin{split}\operatorname{Shi}{\left(z^{q} p \right)} = \frac{z^{q} \sqrt{\pi} p {G_{1, 3}^{1, 1}\left(\begin{matrix} \frac{1}{2} & \\0 & - \frac{1}{2}, - \frac{1}{2} \end{matrix} \middle| {\frac{z^{2 q} p^{2} e^{i \pi}}{4}} \right)}}{4}\end{split}\]

Functions involving \(\operatorname{Chi}\left(z^{q} p\right)\):

\[\begin{split}\operatorname{Chi}\left(z^{q} p\right) = - \frac{\pi^{\frac{3}{2}} {G_{2, 4}^{2, 0}\left(\begin{matrix} & \frac{1}{2}, 1 \\0, 0 & \frac{1}{2}, \frac{1}{2} \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}}{2}\end{split}\]

Functions involving \(\operatorname{E}_{a}\left(z^{q} p\right)\):

\[\begin{split}\operatorname{E}_{a}\left(z^{q} p\right) = {G_{1, 2}^{2, 0}\left(\begin{matrix} & a \\0, a - 1 & \end{matrix} \middle| {z^{q} p} \right)}\end{split}\]

Functions involving \(\operatorname{erf}{\left(z^{q} p \right)}\):

\[\begin{split}\operatorname{erf}{\left(z^{q} p \right)} = \frac{{G_{1, 2}^{1, 1}\left(\begin{matrix} 1 & \\\frac{1}{2} & 0 \end{matrix} \middle| {z^{2 q} p^{2}} \right)}}{\sqrt{\pi}}\end{split}\]

Functions involving \(\operatorname{erfc}{\left(z^{q} p \right)}\):

\[\begin{split}\operatorname{erfc}{\left(z^{q} p \right)} = \frac{{G_{1, 2}^{2, 0}\left(\begin{matrix} & 1 \\0, \frac{1}{2} & \end{matrix} \middle| {z^{2 q} p^{2}} \right)}}{\sqrt{\pi}}\end{split}\]

Functions involving \(\operatorname{erfi}{\left(z^{q} p \right)}\):

\[\begin{split}\operatorname{erfi}{\left(z^{q} p \right)} = \frac{z^{q} p {G_{1, 2}^{1, 1}\left(\begin{matrix} \frac{1}{2} & \\0 & - \frac{1}{2} \end{matrix} \middle| {- z^{2 q} p^{2}} \right)}}{\sqrt{\pi}}\end{split}\]

Functions involving \(S\left(z^{q} p\right)\):

\[\begin{split}S\left(z^{q} p\right) = \frac{{G_{1, 3}^{1, 1}\left(\begin{matrix} 1 & \\\frac{3}{4} & 0, \frac{1}{4} \end{matrix} \middle| {\frac{z^{4 q} \pi^{2} p^{4}}{16}} \right)}}{2}\end{split}\]

Functions involving \(C\left(z^{q} p\right)\):

\[\begin{split}C\left(z^{q} p\right) = \frac{{G_{1, 3}^{1, 1}\left(\begin{matrix} 1 & \\\frac{1}{4} & 0, \frac{3}{4} \end{matrix} \middle| {\frac{z^{4 q} \pi^{2} p^{4}}{16}} \right)}}{2}\end{split}\]

Functions involving \(J_{a}\left(z^{q} p\right)\):

\[\begin{split}J_{a}\left(z^{q} p\right) = {G_{0, 2}^{1, 0}\left(\begin{matrix} & \\\frac{a}{2} & - \frac{a}{2} \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}\end{split}\]

Functions involving \(Y_{a}\left(z^{q} p\right)\):

\[\begin{split}Y_{a}\left(z^{q} p\right) = {G_{1, 3}^{2, 0}\left(\begin{matrix} & - \frac{a}{2} - \frac{1}{2} \\- \frac{a}{2}, \frac{a}{2} & - \frac{a}{2} - \frac{1}{2} \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}\end{split}\]

Functions involving \(I_{a}\left(z^{q} p\right)\):

\[\begin{split}I_{a}\left(z^{q} p\right) = \pi {G_{1, 3}^{1, 0}\left(\begin{matrix} & \frac{a}{2} + \frac{1}{2} \\\frac{a}{2} & - \frac{a}{2}, \frac{a}{2} + \frac{1}{2} \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}\end{split}\]

Functions involving \(K_{a}\left(z^{q} p\right)\):

\[\begin{split}K_{a}\left(z^{q} p\right) = \frac{{G_{0, 2}^{2, 0}\left(\begin{matrix} & \\- \frac{a}{2}, \frac{a}{2} & \end{matrix} \middle| {\frac{z^{2 q} p^{2}}{4}} \right)}}{2}\end{split}\]

Functions involving \(K\left(z^{q} p\right)\):

\[\begin{split}K\left(z^{q} p\right) = \frac{{G_{2, 2}^{1, 2}\left(\begin{matrix} \frac{1}{2}, \frac{1}{2} & \\0 & 0 \end{matrix} \middle| {- z^{q} p} \right)}}{2}\end{split}\]

Functions involving \(E\left(z^{q} p\right)\):

\[\begin{split}E\left(z^{q} p\right) = - \frac{{G_{2, 2}^{1, 2}\left(\begin{matrix} \frac{1}{2}, \frac{3}{2} & \\0 & 0 \end{matrix} \middle| {- z^{q} p} \right)}}{4}\end{split}\]

Internal API Reference

Integrate functions by rewriting them as Meijer G-functions.

There are three user-visible functions that can be used by other parts of the sympy library to solve various integration problems:

• meijerint_indefinite

• meijerint_definite

• meijerint_inversion

They can be used to compute, respectively, indefinite integrals, definite integrals over intervals of the real line, and inverse laplace-type integrals (from c-I*oo to c+I*oo). See the respective docstrings for details.

The main references for this are:

[L] Luke, Y. L. (1969), The Special Functions and Their Approximations,

Volume 1

[R] Kelly B. Roach. Meijer G Function Representations.

In: Proceedings of the 1997 International Symposium on Symbolic and Algebraic Computation, pages 205-211, New York, 1997. ACM.

[P] A. P. Prudnikov, Yu. A. Brychkov and O. I. Marichev (1990).

Integrals and Series: More Special Functions, Vol. 3,. Gordon and Breach Science Publisher

exception sympy.integrals.meijerint._CoeffExpValueError

Exception raised by _get_coeff_exp, for internal use only.

sympy.integrals.meijerint._check_antecedents(g1, g2, x)

Return a condition under which the integral theorem applies.

sympy.integrals.meijerint._check_antecedents_1(g, x, helper=False)

Return a condition under which the mellin transform of g exists. Any power of x has already been absorbed into the G function, so this is just \(\int_0^\infty g\, dx\).

See [L, section 5.6.1]. (Note that s=1.)

If helper is True, only check if the MT exists at infinity, i.e. if \(\int_1^\infty g\, dx\) exists.

sympy.integrals.meijerint._check_antecedents_inversion(g, x)

Check antecedents for the laplace inversion integral.

sympy.integrals.meijerint._condsimp(cond, first=True)

Do naive simplifications on cond.

Explanation

Note that this routine is completely ad-hoc, simplification rules being added as need arises rather than following any logical pattern.

Examples

>>> from sympy.integrals.meijerint import _condsimp as simp
>>> from sympy import Or, Eq
>>> from sympy.abc import x, y
>>> simp(Or(x < y, Eq(x, y)))
x <= y

sympy.integrals.meijerint._create_lookup_table(table)

Add formulae for the function -> meijerg lookup table.

sympy.integrals.meijerint._dummy(name, token, expr, **kwargs)

Return a dummy. This will return the same dummy if the same token+name is requested more than once, and it is not already in expr. This is for being cache-friendly.

sympy.integrals.meijerint._dummy_(name, token, **kwargs)

Return a dummy associated to name and token. Same effect as declaring it globally.

sympy.integrals.meijerint._eval_cond(cond)

Re-evaluate the conditions.

sympy.integrals.meijerint._exponents(expr, x)

Find the exponents of x (not including zero) in expr.

Examples

>>> from sympy.integrals.meijerint import _exponents
>>> from sympy.abc import x, y
>>> from sympy import sin
>>> _exponents(x, x)
{1}
>>> _exponents(x**2, x)
{2}
>>> _exponents(x**2 + x, x)
{1, 2}
>>> _exponents(x**3*sin(x + x**y) + 1/x, x)
{-1, 1, 3, y}

sympy.integrals.meijerint._find_splitting_points(expr, x)

Find numbers a such that a linear substitution x -> x + a would (hopefully) simplify expr.

Examples

>>> from sympy.integrals.meijerint import _find_splitting_points as fsp
>>> from sympy import sin
>>> from sympy.abc import x
>>> fsp(x, x)
{0}
>>> fsp((x-1)**3, x)
{1}
>>> fsp(sin(x+3)*x, x)
{-3, 0}

sympy.integrals.meijerint._flip_g(g)

Turn the G function into one of inverse argument (i.e. G(1/x) -> G’(x))

sympy.integrals.meijerint._functions(expr, x)

Find the types of functions in expr, to estimate the complexity.

sympy.integrals.meijerint._get_coeff_exp(expr, x)

When expr is known to be of the form c*x**b, with c and/or b possibly 1, return c, b.

Examples

>>> from sympy.abc import x, a, b
>>> from sympy.integrals.meijerint import _get_coeff_exp
>>> _get_coeff_exp(a*x**b, x)
(a, b)
>>> _get_coeff_exp(x, x)
(1, 1)
>>> _get_coeff_exp(2*x, x)
(2, 1)
>>> _get_coeff_exp(x**3, x)
(1, 3)

sympy.integrals.meijerint._guess_expansion(f, x)

Try to guess sensible rewritings for integrand f(x).

sympy.integrals.meijerint._inflate_fox_h(g, a)

Let d denote the integrand in the definition of the G function g. Consider the function H which is defined in the same way, but with integrand d/Gamma(a*s) (contour conventions as usual).

If a is rational, the function H can be written as C*G, for a constant C and a G-function G.

This function returns C, G.

sympy.integrals.meijerint._inflate_g(g, n)

Return C, h such that h is a G function of argument z**n and g = C*h.

sympy.integrals.meijerint._int0oo(g1, g2, x)

Express integral from zero to infinity g1*g2 using a G function, assuming the necessary conditions are fulfilled.

Examples

>>> from sympy.integrals.meijerint import _int0oo
>>> from sympy.abc import s, t, m
>>> from sympy import meijerg, S
>>> g1 = meijerg([], [], [-S(1)/2, 0], [], s**2*t/4)
>>> g2 = meijerg([], [], [m/2], [-m/2], t/4)
>>> _int0oo(g1, g2, t)
4*meijerg(((0, 1/2), ()), ((m/2,), (-m/2,)), s**(-2))/s**2

sympy.integrals.meijerint._int0oo_1(g, x)

Evaluate \(\int_0^\infty g\, dx\) using G functions, assuming the necessary conditions are fulfilled.

Examples

>>> from sympy.abc import a, b, c, d, x, y
>>> from sympy import meijerg
>>> from sympy.integrals.meijerint import _int0oo_1
>>> _int0oo_1(meijerg([a], [b], [c], [d], x*y), x)
gamma(-a)*gamma(c + 1)/(y*gamma(-d)*gamma(b + 1))

sympy.integrals.meijerint._int_inversion(g, x, t)

Compute the laplace inversion integral, assuming the formula applies.

sympy.integrals.meijerint._is_analytic(f, x)

Check if f(x), when expressed using G functions on the positive reals, will in fact agree with the G functions almost everywhere

sympy.integrals.meijerint._meijerint_definite_2(f, x)

Try to integrate f dx from zero to infinity.

The body of this function computes various ‘simplifications’ f1, f2, … of f (e.g. by calling expand_mul(), trigexpand() - see _guess_expansion) and calls _meijerint_definite_3 with each of these in succession. If _meijerint_definite_3 succeeds with any of the simplified functions, returns this result.

sympy.integrals.meijerint._meijerint_definite_3(f, x)

Try to integrate f dx from zero to infinity.

This function calls _meijerint_definite_4 to try to compute the integral. If this fails, it tries using linearity.

sympy.integrals.meijerint._meijerint_definite_4(f, x, only_double=False)

Try to integrate f dx from zero to infinity.

Explanation

This function tries to apply the integration theorems found in literature, i.e. it tries to rewrite f as either one or a product of two G-functions.

The parameter only_double is used internally in the recursive algorithm to disable trying to rewrite f as a single G-function.

sympy.integrals.meijerint._meijerint_indefinite_1(f, x)

Helper that does not attempt any substitution.

sympy.integrals.meijerint._mul_args(f)

Return a list L such that Mul(*L) == f.

If f is not a Mul or Pow, L=[f]. If f=g**n for an integer n, L=[g]*n. If f is a Mul, L comes from applying _mul_args to all factors of f.

sympy.integrals.meijerint._mul_as_two_parts(f)

Find all the ways to split f into a product of two terms. Return None on failure.

Explanation

Although the order is canonical from multiset_partitions, this is not necessarily the best order to process the terms. For example, if the case of len(gs) == 2 is removed and multiset is allowed to sort the terms, some tests fail.

Examples

>>> from sympy.integrals.meijerint import _mul_as_two_parts
>>> from sympy import sin, exp, ordered
>>> from sympy.abc import x
>>> list(ordered(_mul_as_two_parts(x*sin(x)*exp(x))))
[(x, exp(x)*sin(x)), (x*exp(x), sin(x)), (x*sin(x), exp(x))]

sympy.integrals.meijerint._my_principal_branch(expr, period, full_pb=False)

Bring expr nearer to its principal branch by removing superfluous factors. This function does not guarantee to yield the principal branch, to avoid introducing opaque principal_branch() objects, unless full_pb=True.

sympy.integrals.meijerint._mytype(

f: Basic,

x: Symbol,

) → tuple[type[Basic], ...]

Create a hashable entity describing the type of f.

sympy.integrals.meijerint._rewrite1(f, x, recursive=True)

Try to rewrite f using a (sum of) single G functions with argument a*x**b. Return fac, po, g such that f = fac*po*g, fac is independent of x. and po = x**s. Here g is a result from _rewrite_single. Return None on failure.

sympy.integrals.meijerint._rewrite2(f, x)

Try to rewrite f as a product of two G functions of arguments a*x**b. Return fac, po, g1, g2 such that f = fac*po*g1*g2, where fac is independent of x and po is x**s. Here g1 and g2 are results of _rewrite_single. Returns None on failure.

sympy.integrals.meijerint._rewrite_inversion(fac, po, g, x)

Absorb po == x**s into g.

sympy.integrals.meijerint._rewrite_saxena(fac, po, g1, g2, x, full_pb=False)

Rewrite the integral fac*po*g1*g2 from 0 to oo in terms of G functions with argument c*x.

Explanation

Return C, f1, f2 such that integral C f1 f2 from 0 to infinity equals integral fac po, g1, g2 from 0 to infinity.

Examples

>>> from sympy.integrals.meijerint import _rewrite_saxena
>>> from sympy.abc import s, t, m
>>> from sympy import meijerg
>>> g1 = meijerg([], [], [0], [], s*t)
>>> g2 = meijerg([], [], [m/2], [-m/2], t**2/4)
>>> r = _rewrite_saxena(1, t**0, g1, g2, t)
>>> r[0]
s/(4*sqrt(pi))
>>> r[1]
meijerg(((), ()), ((-1/2, 0), ()), s**2*t/4)
>>> r[2]
meijerg(((), ()), ((m/2,), (-m/2,)), t/4)

sympy.integrals.meijerint._rewrite_saxena_1(fac, po, g, x)

Rewrite the integral fac*po*g dx, from zero to infinity, as integral fac*G, where G has argument a*x. Note po=x**s. Return fac, G.

sympy.integrals.meijerint._rewrite_single(f, x, recursive=True)

Try to rewrite f as a sum of single G functions of the form C*x**s*G(a*x**b), where b is a rational number and C is independent of x. We guarantee that result.argument.as_coeff_mul(x) returns (a, (x**b,)) or (a, ()). Returns a list of tuples (C, s, G) and a condition cond. Returns None on failure.

sympy.integrals.meijerint._split_mul(f, x)

Split expression f into fac, po, g, where fac is a constant factor, po = x**s for some s independent of s, and g is “the rest”.

Examples

>>> from sympy.integrals.meijerint import _split_mul
>>> from sympy import sin
>>> from sympy.abc import s, x
>>> _split_mul((3*x)**s*sin(x**2)*x, x)
(3**s, x*x**s, sin(x**2))

sympy.integrals.meijerint.meijerint_definite(f, x, a, b)

Integrate f over the interval [a, b], by rewriting it as a product of two G functions, or as a single G function.

Return res, cond, where cond are convergence conditions.

Examples

>>> from sympy.integrals.meijerint import meijerint_definite
>>> from sympy import exp, oo
>>> from sympy.abc import x
>>> meijerint_definite(exp(-x**2), x, -oo, oo)
(sqrt(pi), True)

This function is implemented as a succession of functions meijerint_definite, _meijerint_definite_2, _meijerint_definite_3, _meijerint_definite_4. Each function in the list calls the next one (presumably) several times. This means that calling meijerint_definite can be very costly.

sympy.integrals.meijerint.meijerint_indefinite(f, x)

Compute an indefinite integral of f by rewriting it as a G function.

Examples

>>> from sympy.integrals.meijerint import meijerint_indefinite
>>> from sympy import sin
>>> from sympy.abc import x
>>> meijerint_indefinite(sin(x), x)
-cos(x)

sympy.integrals.meijerint.meijerint_inversion(f, x, t)

Compute the inverse laplace transform \(\int_{c+i\infty}^{c-i\infty} f(x) e^{tx}\, dx\), for real c larger than the real part of all singularities of f.

Note that t is always assumed real and positive.

Return None if the integral does not exist or could not be evaluated.

Examples

>>> from sympy.abc import x, t
>>> from sympy.integrals.meijerint import meijerint_inversion
>>> meijerint_inversion(1/x, x, t)
Heaviside(t)