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

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

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

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.

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

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.

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.

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.

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.):

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.

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:

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\).

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\).

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

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 >= \frac{a}{2}\) or \(p \ge 1, p \ge q\).

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).

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.