(ctx, f, *points, **kwargs)
return +v
r"""
        Computes a single, double or triple integral over a given
        1D interval, 2D rectangle, or 3D cuboid. A basic example::

            >>> from sympy.mpmath import *
            >>> mp.dps = 15; mp.pretty = True
            >>> quad(sin, [0, pi])
            2.0

        A basic 2D integral::

            >>> f = lambda x, y: cos(x+y/2)
            >>> quad(f, [-pi/2, pi/2], [0, pi])
            4.0

        **Interval format**

        The integration range for each dimension may be specified
        using a list or tuple. Arguments are interpreted as follows:

        ``quad(f, [x1, x2])`` -- calculates
        `\int_{x_1}^{x_2} f(x) \, dx`

        ``quad(f, [x1, x2], [y1, y2])`` -- calculates
        `\int_{x_1}^{x_2} \int_{y_1}^{y_2} f(x,y) \, dy \, dx`

        ``quad(f, [x1, x2], [y1, y2], [z1, z2])`` -- calculates
        `\int_{x_1}^{x_2} \int_{y_1}^{y_2} \int_{z_1}^{z_2} f(x,y,z)
        \, dz \, dy \, dx`

        Endpoints may be finite or infinite. An interval descriptor
        may also contain more than two points. In this
        case, the integration is split into subintervals, between
        each pair of consecutive points. This is useful for
        dealing with mid-interval discontinuities, or integrating
        over large intervals where the function is irregular or
        oscillates.

        **Options**

        :func:`~mpmath.quad` recognizes the following keyword arguments:

        *method*
            Chooses integration algorithm (described below).
        *error*
            If set to true, :func:`~mpmath.quad` returns `(v, e)` where `v` is the
            integral and `e` is the estimated error.
        *maxdegree*
            Maximum degree of the quadrature rule to try before
            quitting.
        *verbose*
            Print details about progress.

        **Algorithms**

        Mpmath presently implements two integration algorithms: tanh-sinh
        quadrature and Gauss-Legendre quadrature. These can be selected
        using *method='tanh-sinh'* or *method='gauss-legendre'* or by
        passing the classes *method=TanhSinh*, *method=GaussLegendre*.
        The functions :func:`~mpmath.quadts` and :func:`~mpmath.quadgl` are also available
        as shortcuts.

        Both algorithms have the property that doubling the number of
        evaluation points roughly doubles the accuracy, so both are ideal
        for high precision quadrature (hundreds or thousands of digits).

        At high precision, computing the nodes and weights for the
        integration can be expensive (more expensive than computing the
        function values). To make repeated integrations fast, nodes
        are automatically cached.

        The advantages of the tanh-sinh algorithm are that it tends to
        handle endpoint singularities well, and that the nodes are cheap
        to compute on the first run. For these reasons, it is used by
        :func:`~mpmath.quad` as the default algorithm.

        Gauss-Legendre quadrature often requires fewer function
        evaluations, and is therefore often faster for repeated use, but
        the algorithm does not handle endpoint singularities as well and
        the nodes are more expensive to compute. Gauss-Legendre quadrature
        can be a better choice if the integrand is smooth and repeated
        integrations are required (e.g. for multiple integrals).

        See the documentation for :class:`TanhSinh` and
        :class:`GaussLegendre` for additional details.

        **Examples of 1D integrals**

        Intervals may be infinite or half-infinite. The following two
        examples evaluate the limits of the inverse tangent function
        (`\int 1/(1+x^2) = \tan^{-1} x`), and the Gaussian integral
        `\int_{\infty}^{\infty} \exp(-x^2)\,dx = \sqrt{\pi}`::

            >>> mp.dps = 15
            >>> quad(lambda x: 2/(x**2+1), [0, inf])
            3.14159265358979
            >>> quad(lambda x: exp(-x**2), [-inf, inf])**2
            3.14159265358979

        Integrals can typically be resolved to high precision.
        The following computes 50 digits of `\pi` by integrating the
        area of the half-circle defined by `x^2 + y^2 \le 1`,
        `-1 \le x \le 1`, `y \ge 0`::

            >>> mp.dps = 50
            >>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1])
            3.1415926535897932384626433832795028841971693993751

        One can just as well compute 1000 digits (output truncated)::

            >>> mp.dps = 1000
            >>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1])  #doctest:+ELLIPSIS
            3.141592653589793238462643383279502884...216420198

        Complex integrals are supported. The following computes
        a residue at `z = 0` by integrating counterclockwise along the
        diamond-shaped path from `1` to `+i` to `-1` to `-i` to `1`::

            >>> mp.dps = 15
            >>> chop(quad(lambda z: 1/z, [1,j,-1,-j,1]))
            (0.0 + 6.28318530717959j)

        **Examples of 2D and 3D integrals**

        Here are several nice examples of analytically solvable
        2D integrals (taken from MathWorld [1]) that can be evaluated
        to high precision fairly rapidly by :func:`~mpmath.quad`::

            >>> mp.dps = 30
            >>> f = lambda x, y: (x-1)/((1-x*y)*log(x*y))
            >>> quad(f, [0, 1], [0, 1])
            0.577215664901532860606512090082
            >>> +euler
            0.577215664901532860606512090082

            >>> f = lambda x, y: 1/sqrt(1+x**2+y**2)
            >>> quad(f, [-1, 1], [-1, 1])
            3.17343648530607134219175646705
            >>> 4*log(2+sqrt(3))-2*pi/3
            3.17343648530607134219175646705

            >>> f = lambda x, y: 1/(1-x**2 * y**2)
            >>> quad(f, [0, 1], [0, 1])
            1.23370055013616982735431137498
            >>> pi**2 / 8
            1.23370055013616982735431137498

            >>> quad(lambda x, y: 1/(1-x*y), [0, 1], [0, 1])
            1.64493406684822643647241516665
            >>> pi**2 / 6
            1.64493406684822643647241516665

        Multiple integrals may be done over infinite ranges::

            >>> mp.dps = 15
            >>> print(quad(lambda x,y: exp(-x-y), [0, inf], [1, inf]))
            0.367879441171442
            >>> print(1/e)
            0.367879441171442

        For nonrectangular areas, one can call :func:`~mpmath.quad` recursively.
        For example, we can replicate the earlier example of calculating
        `\pi` by integrating over the unit-circle, and actually use double
        quadrature to actually measure the area circle::

            >>> f = lambda x: quad(lambda y: 1, [-sqrt(1-x**2), sqrt(1-x**2)])
            >>> quad(f, [-1, 1])
            3.14159265358979

        Here is a simple triple integral::

            >>> mp.dps = 15
            >>> f = lambda x,y,z: x*y/(1+z)
            >>> quad(f, [0,1], [0,1], [1,2], method='gauss-legendre')
            0.101366277027041
            >>> (log(3)-log(2))/4
            0.101366277027041

        **Singularities**

        Both tanh-sinh and Gauss-Legendre quadrature are designed to
        integrate smooth (infinitely differentiable) functions. Neither
        algorithm copes well with mid-interval singularities (such as
        mid-interval discontinuities in `f(x)` or `f'(x)`).
        The best solution is to split the integral into parts::

            >>> mp.dps = 15
            >>> quad(lambda x: abs(sin(x)), [0, 2*pi])   # Bad
            3.99900894176779
            >>> quad(lambda x: abs(sin(x)), [0, pi, 2*pi])  # Good
            4.0

        The tanh-sinh rule often works well for integrands having a
        singularity at one or both endpoints::

            >>> mp.dps = 15
            >>> quad(log, [0, 1], method='tanh-sinh')  # Good
            -1.0
            >>> quad(log, [0, 1], method='gauss-legendre')  # Bad
            -0.999932197413801

        However, the result may still be inaccurate for some functions::

            >>> quad(lambda x: 1/sqrt(x), [0, 1], method='tanh-sinh')
            1.99999999946942

        This problem is not due to the quadrature rule per se, but to
        numerical amplification of errors in the nodes. The problem can be
        circumvented by temporarily increasing the precision::

            >>> mp.dps = 30
            >>> a = quad(lambda x: 1/sqrt(x), [0, 1], method='tanh-sinh')
            >>> mp.dps = 15
            >>> +a
            2.0

        **Highly variable functions**

        For functions that are smooth (in the sense of being infinitely
        differentiable) but contain sharp mid-interval peaks or many
        "bumps", :func:`~mpmath.quad` may fail to provide full accuracy. For
        example, with default settings, :func:`~mpmath.quad` is able to integrate
        `\sin(x)` accurately over an interval of length 100 but not over
        length 1000::

            >>> quad(sin, [0, 100]); 1-cos(100)   # Good
            0.137681127712316
            0.137681127712316
            >>> quad(sin, [0, 1000]); 1-cos(1000)   # Bad
            -37.8587612408485
            0.437620923709297

        One solution is to break the integration into 10 intervals of
        length 100::

            >>> quad(sin, linspace(0, 1000, 10))   # Good
            0.437620923709297

        Another is to increase the degree of the quadrature::

            >>> quad(sin, [0, 1000], maxdegree=10)   # Also good
            0.437620923709297

        Whether splitting the interval or increasing the degree is
        more efficient differs from case to case. Another example is the
        function `1/(1+x^2)`, which has a sharp peak centered around
        `x = 0`::

            >>> f = lambda x: 1/(1+x**2)
            >>> quad(f, [-100, 100])   # Bad
            3.64804647105268
            >>> quad(f, [-100, 100], maxdegree=10)   # Good
            3.12159332021646
            >>> quad(f, [-100, 0, 100])   # Also good
            3.12159332021646

        **References**

        1. http://mathworld.wolfram.com/DoubleIntegral.html
r"""
        Computes a single, double or triple integral over a given
        1D interval, 2D rectangle, or 3D cuboid. A basic example::
r Computes a single double or triple integral over a given 1D interval 2D rectangle or 3D cuboid . A basic example ::
def quad(ctx, f, *points, **kwargs):
        r"""
        Computes a single, double or triple integral over a given
        1D interval, 2D rectangle, or 3D cuboid. A basic example::

            >>> from sympy.mpmath import *
            >>> mp.dps = 15; mp.pretty = True
            >>> quad(sin, [0, pi])
            2.0

        A basic 2D integral::

            >>> f = lambda x, y: cos(x+y/2)
            >>> quad(f, [-pi/2, pi/2], [0, pi])
            4.0

        **Interval format**

        The integration range for each dimension may be specified
        using a list or tuple. Arguments are interpreted as follows:

        ``quad(f, [x1, x2])`` -- calculates
        `\int_{x_1}^{x_2} f(x) \, dx`

        ``quad(f, [x1, x2], [y1, y2])`` -- calculates
        `\int_{x_1}^{x_2} \int_{y_1}^{y_2} f(x,y) \, dy \, dx`

        ``quad(f, [x1, x2], [y1, y2], [z1, z2])`` -- calculates
        `\int_{x_1}^{x_2} \int_{y_1}^{y_2} \int_{z_1}^{z_2} f(x,y,z)
        \, dz \, dy \, dx`

        Endpoints may be finite or infinite. An interval descriptor
        may also contain more than two points. In this
        case, the integration is split into subintervals, between
        each pair of consecutive points. This is useful for
        dealing with mid-interval discontinuities, or integrating
        over large intervals where the function is irregular or
        oscillates.

        **Options**

        :func:`~mpmath.quad` recognizes the following keyword arguments:

        *method*
            Chooses integration algorithm (described below).
        *error*
            If set to true, :func:`~mpmath.quad` returns `(v, e)` where `v` is the
            integral and `e` is the estimated error.
        *maxdegree*
            Maximum degree of the quadrature rule to try before
            quitting.
        *verbose*
            Print details about progress.

        **Algorithms**

        Mpmath presently implements two integration algorithms: tanh-sinh
        quadrature and Gauss-Legendre quadrature. These can be selected
        using *method='tanh-sinh'* or *method='gauss-legendre'* or by
        passing the classes *method=TanhSinh*, *method=GaussLegendre*.
        The functions :func:`~mpmath.quadts` and :func:`~mpmath.quadgl` are also available
        as shortcuts.

        Both algorithms have the property that doubling the number of
        evaluation points roughly doubles the accuracy, so both are ideal
        for high precision quadrature (hundreds or thousands of digits).

        At high precision, computing the nodes and weights for the
        integration can be expensive (more expensive than computing the
        function values). To make repeated integrations fast, nodes
        are automatically cached.

        The advantages of the tanh-sinh algorithm are that it tends to
        handle endpoint singularities well, and that the nodes are cheap
        to compute on the first run. For these reasons, it is used by
        :func:`~mpmath.quad` as the default algorithm.

        Gauss-Legendre quadrature often requires fewer function
        evaluations, and is therefore often faster for repeated use, but
        the algorithm does not handle endpoint singularities as well and
        the nodes are more expensive to compute. Gauss-Legendre quadrature
        can be a better choice if the integrand is smooth and repeated
        integrations are required (e.g. for multiple integrals).

        See the documentation for :class:`TanhSinh` and
        :class:`GaussLegendre` for additional details.

        **Examples of 1D integrals**

        Intervals may be infinite or half-infinite. The following two
        examples evaluate the limits of the inverse tangent function
        (`\int 1/(1+x^2) = \tan^{-1} x`), and the Gaussian integral
        `\int_{\infty}^{\infty} \exp(-x^2)\,dx = \sqrt{\pi}`::

            >>> mp.dps = 15
            >>> quad(lambda x: 2/(x**2+1), [0, inf])
            3.14159265358979
            >>> quad(lambda x: exp(-x**2), [-inf, inf])**2
            3.14159265358979

        Integrals can typically be resolved to high precision.
        The following computes 50 digits of `\pi` by integrating the
        area of the half-circle defined by `x^2 + y^2 \le 1`,
        `-1 \le x \le 1`, `y \ge 0`::

            >>> mp.dps = 50
            >>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1])
            3.1415926535897932384626433832795028841971693993751

        One can just as well compute 1000 digits (output truncated)::

            >>> mp.dps = 1000
            >>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1])  #doctest:+ELLIPSIS
            3.141592653589793238462643383279502884...216420198

        Complex integrals are supported. The following computes
        a residue at `z = 0` by integrating counterclockwise along the
        diamond-shaped path from `1` to `+i` to `-1` to `-i` to `1`::

            >>> mp.dps = 15
            >>> chop(quad(lambda z: 1/z, [1,j,-1,-j,1]))
            (0.0 + 6.28318530717959j)

        **Examples of 2D and 3D integrals**

        Here are several nice examples of analytically solvable
        2D integrals (taken from MathWorld [1]) that can be evaluated
        to high precision fairly rapidly by :func:`~mpmath.quad`::

            >>> mp.dps = 30
            >>> f = lambda x, y: (x-1)/((1-x*y)*log(x*y))
            >>> quad(f, [0, 1], [0, 1])
            0.577215664901532860606512090082
            >>> +euler
            0.577215664901532860606512090082

            >>> f = lambda x, y: 1/sqrt(1+x**2+y**2)
            >>> quad(f, [-1, 1], [-1, 1])
            3.17343648530607134219175646705
            >>> 4*log(2+sqrt(3))-2*pi/3
            3.17343648530607134219175646705

            >>> f = lambda x, y: 1/(1-x**2 * y**2)
            >>> quad(f, [0, 1], [0, 1])
            1.23370055013616982735431137498
            >>> pi**2 / 8
            1.23370055013616982735431137498

            >>> quad(lambda x, y: 1/(1-x*y), [0, 1], [0, 1])
            1.64493406684822643647241516665
            >>> pi**2 / 6
            1.64493406684822643647241516665

        Multiple integrals may be done over infinite ranges::

            >>> mp.dps = 15
            >>> print(quad(lambda x,y: exp(-x-y), [0, inf], [1, inf]))
            0.367879441171442
            >>> print(1/e)
            0.367879441171442

        For nonrectangular areas, one can call :func:`~mpmath.quad` recursively.
        For example, we can replicate the earlier example of calculating
        `\pi` by integrating over the unit-circle, and actually use double
        quadrature to actually measure the area circle::

            >>> f = lambda x: quad(lambda y: 1, [-sqrt(1-x**2), sqrt(1-x**2)])
            >>> quad(f, [-1, 1])
            3.14159265358979

        Here is a simple triple integral::

            >>> mp.dps = 15
            >>> f = lambda x,y,z: x*y/(1+z)
            >>> quad(f, [0,1], [0,1], [1,2], method='gauss-legendre')
            0.101366277027041
            >>> (log(3)-log(2))/4
            0.101366277027041

        **Singularities**

        Both tanh-sinh and Gauss-Legendre quadrature are designed to
        integrate smooth (infinitely differentiable) functions. Neither
        algorithm copes well with mid-interval singularities (such as
        mid-interval discontinuities in `f(x)` or `f'(x)`).
        The best solution is to split the integral into parts::

            >>> mp.dps = 15
            >>> quad(lambda x: abs(sin(x)), [0, 2*pi])   # Bad
            3.99900894176779
            >>> quad(lambda x: abs(sin(x)), [0, pi, 2*pi])  # Good
            4.0

        The tanh-sinh rule often works well for integrands having a
        singularity at one or both endpoints::

            >>> mp.dps = 15
            >>> quad(log, [0, 1], method='tanh-sinh')  # Good
            -1.0
            >>> quad(log, [0, 1], method='gauss-legendre')  # Bad
            -0.999932197413801

        However, the result may still be inaccurate for some functions::

            >>> quad(lambda x: 1/sqrt(x), [0, 1], method='tanh-sinh')
            1.99999999946942

        This problem is not due to the quadrature rule per se, but to
        numerical amplification of errors in the nodes. The problem can be
        circumvented by temporarily increasing the precision::

            >>> mp.dps = 30
            >>> a = quad(lambda x: 1/sqrt(x), [0, 1], method='tanh-sinh')
            >>> mp.dps = 15
            >>> +a
            2.0

        **Highly variable functions**

        For functions that are smooth (in the sense of being infinitely
        differentiable) but contain sharp mid-interval peaks or many
        "bumps", :func:`~mpmath.quad` may fail to provide full accuracy. For
        example, with default settings, :func:`~mpmath.quad` is able to integrate
        `\sin(x)` accurately over an interval of length 100 but not over
        length 1000::

            >>> quad(sin, [0, 100]); 1-cos(100)   # Good
            0.137681127712316
            0.137681127712316
            >>> quad(sin, [0, 1000]); 1-cos(1000)   # Bad
            -37.8587612408485
            0.437620923709297

        One solution is to break the integration into 10 intervals of
        length 100::

            >>> quad(sin, linspace(0, 1000, 10))   # Good
            0.437620923709297

        Another is to increase the degree of the quadrature::

            >>> quad(sin, [0, 1000], maxdegree=10)   # Also good
            0.437620923709297

        Whether splitting the interval or increasing the degree is
        more efficient differs from case to case. Another example is the
        function `1/(1+x^2)`, which has a sharp peak centered around
        `x = 0`::

            >>> f = lambda x: 1/(1+x**2)
            >>> quad(f, [-100, 100])   # Bad
            3.64804647105268
            >>> quad(f, [-100, 100], maxdegree=10)   # Good
            3.12159332021646
            >>> quad(f, [-100, 0, 100])   # Also good
            3.12159332021646

        **References**

        1. http://mathworld.wolfram.com/DoubleIntegral.html

        """
        rule = kwargs.get('method', 'tanh-sinh')
        if type(rule) is str:
            if rule == 'tanh-sinh':
                rule = ctx._tanh_sinh
            elif rule == 'gauss-legendre':
                rule = ctx._gauss_legendre
            else:
                raise ValueError("unknown quadrature rule: %s" % rule)
        else:
            rule = rule(ctx)
        verbose = kwargs.get('verbose')
        dim = len(points)
        orig = prec = ctx.prec
        epsilon = ctx.eps/8
        m = kwargs.get('maxdegree') or rule.guess_degree(prec)
        points = [ctx._as_points(p) for p in points]
        try:
            ctx.prec += 20
            if dim == 1:
                v, err = rule.summation(f, points[0], prec, epsilon, m, verbose)
            elif dim == 2:
                v, err = rule.summation(lambda x: \
                        rule.summation(lambda y: f(x,y), \
                        points[1], prec, epsilon, m)[0],
                    points[0], prec, epsilon, m, verbose)
            elif dim == 3:
                v, err = rule.summation(lambda x: \
                        rule.summation(lambda y: \
                            rule.summation(lambda z: f(x,y,z), \
                            points[2], prec, epsilon, m)[0],
                        points[1], prec, epsilon, m)[0],
                    points[0], prec, epsilon, m, verbose)
            else:
                raise NotImplementedError("quadrature must have dim 1, 2 or 3")
        finally:
            ctx.prec = orig
        if kwargs.get("error"):
            return +v, err
        return +v
def quad ( ctx , f , * points , * * kwargs ) : rule = kwargs . get ( 'method' , 'tanh-sinh' ) if type ( rule ) is str : if rule == 'tanh-sinh' : rule = ctx . _tanh_sinh elif rule == 'gauss-legendre' : rule = ctx . _gauss_legendre else : raise ValueError ( "unknown quadrature rule: %s" % rule ) else : rule = rule ( ctx ) verbose = kwargs . get ( 'verbose' ) dim = len ( points ) orig = prec = ctx . prec epsilon = ctx . eps / 8 m = kwargs . get ( 'maxdegree' ) or rule . guess_degree ( prec ) points = [ ctx . _as_points ( p ) for p in points ] try : ctx . prec += 20 if dim == 1 : v , err = rule . summation ( f , points [ 0 ] , prec , epsilon , m , verbose ) elif dim == 2 : v , err = rule . summation ( lambda x : rule . summation ( lambda y : f ( x , y ) , points [ 1 ] , prec , epsilon , m ) [ 0 ] , points [ 0 ] , prec , epsilon , m , verbose ) elif dim == 3 : v , err = rule . summation ( lambda x : rule . summation ( lambda y : rule . summation ( lambda z : f ( x , y , z ) , points [ 2 ] , prec , epsilon , m ) [ 0 ] , points [ 1 ] , prec , epsilon , m ) [ 0 ] , points [ 0 ] , prec , epsilon , m , verbose ) else : raise NotImplementedError ( "quadrature must have dim 1, 2 or 3" ) finally : ctx . prec = orig if kwargs . get ( "error" ) : return + v , err return + v