docstr: Computes a single, double or triple integral over a given
  1D interval, 2D rectangle, or 3D cuboid. A basic example::
  
      >>> from 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...216420199
  
  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
docstr_tok: ["Computes","a","single",",","double","or","triple","integral","over","a","given","1D","interval",",","2D","rectangle",",","or","3D","cuboid",".","A","basic","example",":",":",">",">",">","from","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","...","216420199","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"]
code: 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 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...216420199
  
      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
  
imports: import math
  from ..libmp.backend import xrange
label_desc: Use image node_id 2 for calling the QuadratureMethods obj's underlying member method code with example usage: obj.quad(ctx, f) and returns: unknown, unknown,  err