text: from sympy.concrete.expr_with_limits import AddWithLimits
  from sympy.core.add import Add
  from sympy.core.basic import Basic
  from sympy.core.compatibility import is_sequence
  from sympy.core.containers import Tuple
  from sympy.core.expr import Expr
  from sympy.core.function import diff
  from sympy.core.logic import fuzzy_bool
  from sympy.core.mul import Mul
  from sympy.core.numbers import oo, pi
  from sympy.core.relational import Ne
  from sympy.core.singleton import S
  from sympy.core.symbol import (Dummy, Symbol, Wild)
  from sympy.core.sympify import sympify
  from sympy.functions import Piecewise, sqrt, piecewise_fold, tan, cot, atan
  from sympy.functions.elementary.exponential import log
  from sympy.functions.elementary.integers import floor
  from sympy.functions.elementary.complexes import Abs, sign
  from sympy.functions.elementary.miscellaneous import Min, Max
  from sympy.integrals.manualintegrate import manualintegrate
  from sympy.integrals.trigonometry import trigintegrate
  from sympy.integrals.meijerint import meijerint_definite, meijerint_indefinite
  from sympy.matrices import MatrixBase
  from sympy.polys import Poly, PolynomialError
  from sympy.series import limit
  from sympy.series.order import Order
  from sympy.series.formal import FormalPowerSeries
  from sympy.simplify.fu import sincos_to_sum
  from sympy.utilities.misc import filldedent
  from sympy.utilities.exceptions import SymPyDeprecationWarning
  
  
  class Integral(AddWithLimits):
      """Represents unevaluated integral."""
  
      __slots__ = ('is_commutative',)
  
      def __new__(cls, function, *symbols, **assumptions):
          """Create an unevaluated integral.
  
          Arguments are an integrand followed by one or more limits.
  
          If no limits are given and there is only one free symbol in the
          expression, that symbol will be used, otherwise an error will be
          raised.
  
          >>> from sympy import Integral
          >>> from sympy.abc import x, y
          >>> Integral(x)
          Integral(x, x)
          >>> Integral(y)
          Integral(y, y)
  
          When limits are provided, they are interpreted as follows (using
          ``x`` as though it were the variable of integration):
  
              (x,) or x - indefinite integral
              (x, a) - "evaluate at" integral is an abstract antiderivative
              (x, a, b) - definite integral
  
          The ``as_dummy`` method can be used to see which symbols cannot be
          targeted by subs: those with a prepended underscore cannot be
          changed with ``subs``. (Also, the integration variables themselves --
          the first element of a limit -- can never be changed by subs.)
  
          >>> i = Integral(x, x)
          >>> at = Integral(x, (x, x))
          >>> i.as_dummy()
          Integral(x, x)
          >>> at.as_dummy()
          Integral(_0, (_0, x))
  
          """
  
          #This will help other classes define their own definitions
          #of behaviour with Integral.
          if hasattr(function, '_eval_Integral'):
              return function._eval_Integral(*symbols, **assumptions)
  
          if isinstance(function, Poly):
              SymPyDeprecationWarning(
                  feature="Using integrate/Integral with Poly",
                  issue=18613,
                  deprecated_since_version="1.6",
                  useinstead="the as_expr or integrate methods of Poly").warn()
  
          obj = AddWithLimits.__new__(cls, function, *symbols, **assumptions)
          return obj
  
      def __getnewargs__(self):
          return (self.function,) + tuple([tuple(xab) for xab in self.limits])
  
      @property
      def free_symbols(self):
          """
          This method returns the symbols that will exist when the
          integral is evaluated. This is useful if one is trying to
          determine whether an integral depends on a certain
          symbol or not.
  
          Examples
          ========
  
          >>> from sympy import Integral
          >>> from sympy.abc import x, y
          >>> Integral(x, (x, y, 1)).free_symbols
          {y}
  
          See Also
          ========
  
          sympy.concrete.expr_with_limits.ExprWithLimits.function
          sympy.concrete.expr_with_limits.ExprWithLimits.limits
          sympy.concrete.expr_with_limits.ExprWithLimits.variables
          """
          return AddWithLimits.free_symbols.fget(self)
  
      def _eval_is_zero(self):
          # This is a very naive and quick test, not intended to do the integral to
          # answer whether it is zero or not, e.g. Integral(sin(x), (x, 0, 2*pi))
          # is zero but this routine should return None for that case. But, like
          # Mul, there are trivial situations for which the integral will be
          # zero so we check for those.
          if self.function.is_zero:
              return True
          got_none = False
          for l in self.limits:
              if len(l) == 3:
                  z = (l[1] == l[2]) or (l[1] - l[2]).is_zero
                  if z:
                      return True
                  elif z is None:
                      got_none = True
          free = self.function.free_symbols
          for xab in self.limits:
              if len(xab) == 1:
                  free.add(xab[0])
                  continue
              if len(xab) == 2 and xab[0] not in free:
                  if xab[1].is_zero:
                      return True
                  elif xab[1].is_zero is None:
                      got_none = True
              # take integration symbol out of free since it will be replaced
              # with the free symbols in the limits
              free.discard(xab[0])
              # add in the new symbols
              for i in xab[1:]:
                  free.update(i.free_symbols)
          if self.function.is_zero is False and got_none is False:
              return False
  
      def transform(self, x, u):
          r"""
          Performs a change of variables from `x` to `u` using the relationship
          given by `x` and `u` which will define the transformations `f` and `F`
          (which are inverses of each other) as follows:
  
          1) If `x` is a Symbol (which is a variable of integration) then `u`
             will be interpreted as some function, f(u), with inverse F(u).
             This, in effect, just makes the substitution of x with f(x).
  
          2) If `u` is a Symbol then `x` will be interpreted as some function,
             F(x), with inverse f(u). This is commonly referred to as
             u-substitution.
  
          Once f and F have been identified, the transformation is made as
          follows:
  
          .. math:: \int_a^b x \mathrm{d}x \rightarrow \int_{F(a)}^{F(b)} f(x)
                    \frac{\mathrm{d}}{\mathrm{d}x}
  
          where `F(x)` is the inverse of `f(x)` and the limits and integrand have
          been corrected so as to retain the same value after integration.
  
          Notes
          =====
  
          The mappings, F(x) or f(u), must lead to a unique integral. Linear
          or rational linear expression, `2*x`, `1/x` and `sqrt(x)`, will
          always work; quadratic expressions like `x**2 - 1` are acceptable
          as long as the resulting integrand does not depend on the sign of
          the solutions (see examples).
  
          The integral will be returned unchanged if `x` is not a variable of
          integration.
  
          `x` must be (or contain) only one of of the integration variables. If
          `u` has more than one free symbol then it should be sent as a tuple
          (`u`, `uvar`) where `uvar` identifies which variable is replacing
          the integration variable.
          XXX can it contain another integration variable?
  
          Examples
          ========
  
          >>> from sympy.abc import a, x, u
          >>> from sympy import Integral, cos, sqrt
  
          >>> i = Integral(x*cos(x**2 - 1), (x, 0, 1))
  
          transform can change the variable of integration
  
          >>> i.transform(x, u)
          Integral(u*cos(u**2 - 1), (u, 0, 1))
  
          transform can perform u-substitution as long as a unique
          integrand is obtained:
  
          >>> i.transform(x**2 - 1, u)
          Integral(cos(u)/2, (u, -1, 0))
  
          This attempt fails because x = +/-sqrt(u + 1) and the
          sign does not cancel out of the integrand:
  
          >>> Integral(cos(x**2 - 1), (x, 0, 1)).transform(x**2 - 1, u)
          Traceback (most recent call last):
          ...
          ValueError:
          The mapping between F(x) and f(u) did not give a unique integrand.
  
          transform can do a substitution. Here, the previous
          result is transformed back into the original expression
          using "u-substitution":
  
          >>> ui = _
          >>> _.transform(sqrt(u + 1), x) == i
          True
  
          We can accomplish the same with a regular substitution:
  
          >>> ui.transform(u, x**2 - 1) == i
          True
  
          If the `x` does not contain a symbol of integration then
          the integral will be returned unchanged. Integral `i` does
          not have an integration variable `a` so no change is made:
  
          >>> i.transform(a, x) == i
          True
  
          When `u` has more than one free symbol the symbol that is
          replacing `x` must be identified by passing `u` as a tuple:
  
          >>> Integral(x, (x, 0, 1)).transform(x, (u + a, u))
          Integral(a + u, (u, -a, 1 - a))
          >>> Integral(x, (x, 0, 1)).transform(x, (u + a, a))
          Integral(a + u, (a, -u, 1 - u))
  
          See Also
          ========
  
          sympy.concrete.expr_with_limits.ExprWithLimits.variables : Lists the integration variables
          as_dummy : Replace integration variables with dummy ones
          """
          from sympy.solvers.solvers import solve, posify
          d = Dummy('d')
  
          xfree = x.free_symbols.intersection(self.variables)
          if len(xfree) > 1:
              raise ValueError(
                  'F(x) can only contain one of: %s' % self.variables)
          xvar = xfree.pop() if xfree else d
  
          if xvar not in self.variables:
              return self
  
          u = sympify(u)
          if isinstance(u, Expr):
              ufree = u.free_symbols
              if len(ufree) == 0:
                  raise ValueError(filldedent('''
                  f(u) cannot be a constant'''))
              if len(ufree) > 1:
                  raise ValueError(filldedent('''
                  When f(u) has more than one free symbol, the one replacing x
                  must be identified: pass f(u) as (f(u), u)'''))
              uvar = ufree.pop()
          else:
              u, uvar = u
              if uvar not in u.free_symbols:
                  raise ValueError(filldedent('''
                  Expecting a tuple (expr, symbol) where symbol identified
                  a free symbol in expr, but symbol is not in expr's free
                  symbols.'''))
              if not isinstance(uvar, Symbol):
                  # This probably never evaluates to True
                  raise ValueError(filldedent('''
                  Expecting a tuple (expr, symbol) but didn't get
                  a symbol; got %s''' % uvar))
  
          if x.is_Symbol and u.is_Symbol:
              return self.xreplace({x: u})
  
          if not x.is_Symbol and not u.is_Symbol:
              raise ValueError('either x or u must be a symbol')
  
          if uvar == xvar:
              return self.transform(x, (u.subs(uvar, d), d)).xreplace({d: uvar})
  
          if uvar in self.limits:
              raise ValueError(filldedent('''
              u must contain the same variable as in x
              or a variable that is not already an integration variable'''))
  
          if not x.is_Symbol:
              F = [x.subs(xvar, d)]
              soln = solve(u - x, xvar, check=False)
              if not soln:
                  raise ValueError('no solution for solve(F(x) - f(u), x)')
              f = [fi.subs(uvar, d) for fi in soln]
          else:
              f = [u.subs(uvar, d)]
              pdiff, reps = posify(u - x)
              puvar = uvar.subs([(v, k) for k, v in reps.items()])
              soln = [s.subs(reps) for s in solve(pdiff, puvar)]
              if not soln:
                  raise ValueError('no solution for solve(F(x) - f(u), u)')
              F = [fi.subs(xvar, d) for fi in soln]
  
          newfuncs = {(self.function.subs(xvar, fi)*fi.diff(d)
                          ).subs(d, uvar) for fi in f}
          if len(newfuncs) > 1:
              raise ValueError(filldedent('''
              The mapping between F(x) and f(u) did not give
              a unique integrand.'''))
          newfunc = newfuncs.pop()
  
          def _calc_limit_1(F, a, b):
              """
              replace d with a, using subs if possible, otherwise limit
              where sign of b is considered
              """
              wok = F.subs(d, a)
              if wok is S.NaN or wok.is_finite is False and a.is_finite:
                  return limit(sign(b)*F, d, a)
              return wok
  
          def _calc_limit(a, b):
              """
              replace d with a, using subs if possible, otherwise limit
              where sign of b is considered
              """
              avals = list({_calc_limit_1(Fi, a, b) for Fi in F})
              if len(avals) > 1:
                  raise ValueError(filldedent('''
                  The mapping between F(x) and f(u) did not
                  give a unique limit.'''))
              return avals[0]
  
          newlimits = []
          for xab in self.limits:
              sym = xab[0]
              if sym == xvar:
                  if len(xab) == 3:
                      a, b = xab[1:]
                      a, b = _calc_limit(a, b), _calc_limit(b, a)
                      if fuzzy_bool(a - b > 0):
                          a, b = b, a
                          newfunc = -newfunc
                      newlimits.append((uvar, a, b))
                  elif len(xab) == 2:
                      a = _calc_limit(xab[1], 1)
                      newlimits.append((uvar, a))
                  else:
                      newlimits.append(uvar)
              else:
                  newlimits.append(xab)
  
          return self.func(newfunc, *newlimits)
  
      def doit(self, **hints):
          """
          Perform the integration using any hints given.
  
          Examples
          ========
  
          >>> from sympy import Piecewise, S
          >>> from sympy.abc import x, t
          >>> p = x**2 + Piecewise((0, x/t < 0), (1, True))
          >>> p.integrate((t, S(4)/5, 1), (x, -1, 1))
          1/3
  
          See Also
          ========
  
          sympy.integrals.trigonometry.trigintegrate
          sympy.integrals.heurisch.heurisch
          sympy.integrals.rationaltools.ratint
          as_sum : Approximate the integral using a sum
          """
          from sympy.concrete.summations import Sum
          if not hints.get('integrals', True):
              return self
  
          deep = hints.get('deep', True)
          meijerg = hints.get('meijerg', None)
          conds = hints.get('conds', 'piecewise')
          risch = hints.get('risch', None)
          heurisch = hints.get('heurisch', None)
          manual = hints.get('manual', None)
          if len(list(filter(None, (manual, meijerg, risch, heurisch)))) > 1:
              raise ValueError("At most one of manual, meijerg, risch, heurisch can be True")
          elif manual:
              meijerg = risch = heurisch = False
          elif meijerg:
              manual = risch = heurisch = False
          elif risch:
              manual = meijerg = heurisch = False
          elif heurisch:
              manual = meijerg = risch = False
          eval_kwargs = dict(meijerg=meijerg, risch=risch, manual=manual, heurisch=heurisch,
              conds=conds)
  
          if conds not in ['separate', 'piecewise', 'none']:
              raise ValueError('conds must be one of "separate", "piecewise", '
                               '"none", got: %s' % conds)
  
          if risch and any(len(xab) > 1 for xab in self.limits):
              raise ValueError('risch=True is only allowed for indefinite integrals.')
  
          # check for the trivial zero
          if self.is_zero:
              return S.Zero
  
          # hacks to handle integrals of
          # nested summations
          if isinstance(self.function, Sum):
              if any(v in self.function.limits[0] for v in self.variables):
                  raise ValueError('Limit of the sum cannot be an integration variable.')
              if any(l.is_infinite for l in self.function.limits[0][1:]):
                  return self
              _i = self
              _sum = self.function
              return _sum.func(_i.func(_sum.function, *_i.limits).doit(), *_sum.limits).doit()
  
          # now compute and check the function
          function = self.function
          if deep:
              function = function.doit(**hints)
          if function.is_zero:
              return S.Zero
  
          # hacks to handle special cases
          if isinstance(function, MatrixBase):
              return function.applyfunc(
                  lambda f: self.func(f, self.limits).doit(**hints))
  
          if isinstance(function, FormalPowerSeries):
              if len(self.limits) > 1:
                  raise NotImplementedError
              xab = self.limits[0]
              if len(xab) > 1:
                  return function.integrate(xab, **eval_kwargs)
              else:
                  return function.integrate(xab[0], **eval_kwargs)
  
          # There is no trivial answer and special handling
          # is done so continue
  
          # first make sure any definite limits have integration
          # variables with matching assumptions
          reps = {}
          for xab in self.limits:
              if len(xab) != 3:
                  continue
              x, a, b = xab
              l = (a, b)
              if all(i.is_nonnegative for i in l) and not x.is_nonnegative:
                  d = Dummy(positive=True)
              elif all(i.is_nonpositive for i in l) and not x.is_nonpositive:
                  d = Dummy(negative=True)
              elif all(i.is_real for i in l) and not x.is_real:
                  d = Dummy(real=True)
              else:
                  d = None
              if d:
                  reps[x] = d
          if reps:
              undo = {v: k for k, v in reps.items()}
              did = self.xreplace(reps).doit(**hints)
              if type(did) is tuple:  # when separate=True
                  did = tuple([i.xreplace(undo) for i in did])
              else:
                  did = did.xreplace(undo)
              return did
  
          # continue with existing assumptions
          undone_limits = []
          # ulj = free symbols of any undone limits' upper and lower limits
          ulj = set()
          for xab in self.limits:
              # compute uli, the free symbols in the
              # Upper and Lower limits of limit I
              if len(xab) == 1:
                  uli = set(xab[:1])
              elif len(xab) == 2:
                  uli = xab[1].free_symbols
              elif len(xab) == 3:
                  uli = xab[1].free_symbols.union(xab[2].free_symbols)
              # this integral can be done as long as there is no blocking
              # limit that has been undone. An undone limit is blocking if
              # it contains an integration variable that is in this limit's
              # upper or lower free symbols or vice versa
              if xab[0] in ulj or any(v[0] in uli for v in undone_limits):
                  undone_limits.append(xab)
                  ulj.update(uli)
                  function = self.func(*([function] + [xab]))
                  factored_function = function.factor()
                  if not isinstance(factored_function, Integral):
                      function = factored_function
                  continue
  
              if function.has(Abs, sign) and (
                  (len(xab) < 3 and all(x.is_extended_real for x in xab)) or
                  (len(xab) == 3 and all(x.is_extended_real and not x.is_infinite for
                   x in xab[1:]))):
                      # some improper integrals are better off with Abs
                      xr = Dummy("xr", real=True)
                      function = (function.xreplace({xab[0]: xr})
                          .rewrite(Piecewise).xreplace({xr: xab[0]}))
              elif function.has(Min, Max):
                  function = function.rewrite(Piecewise)
              if (function.has(Piecewise) and
                  not isinstance(function, Piecewise)):
                      function = piecewise_fold(function)
              if isinstance(function, Piecewise):
                  if len(xab) == 1:
                      antideriv = function._eval_integral(xab[0],
                          **eval_kwargs)
                  else:
                      antideriv = self._eval_integral(
                          function, xab[0], **eval_kwargs)
              else:
                  # There are a number of tradeoffs in using the
                  # Meijer G method. It can sometimes be a lot faster
                  # than other methods, and sometimes slower. And
                  # there are certain types of integrals for which it
                  # is more likely to work than others. These
                  # heuristics are incorporated in deciding what
                  # integration methods to try, in what order. See the
                  # integrate() docstring for details.
                  def try_meijerg(function, xab):
                      ret = None
                      if len(xab) == 3 and meijerg is not False:
                          x, a, b = xab
                          try:
                              res = meijerint_definite(function, x, a, b)
                          except NotImplementedError:
                              from sympy.integrals.meijerint import _debug
                              _debug('NotImplementedError '
                                  'from meijerint_definite')
                              res = None
                          if res is not None:
                              f, cond = res
                              if conds == 'piecewise':
                                  ret = Piecewise(
                                      (f, cond),
                                      (self.func(
                                      function, (x, a, b)), True))
                              elif conds == 'separate':
                                  if len(self.limits) != 1:
                                      raise ValueError(filldedent('''
                                          conds=separate not supported in
                                          multiple integrals'''))
                                  ret = f, cond
                              else:
                                  ret = f
                      return ret
  
                  meijerg1 = meijerg
                  if (meijerg is not False and
                          len(xab) == 3 and xab[1].is_extended_real and xab[2].is_extended_real
                          and not function.is_Poly and
                          (xab[1].has(oo, -oo) or xab[2].has(oo, -oo))):
                      ret = try_meijerg(function, xab)
                      if ret is not None:
                          function = ret
                          continue
                      meijerg1 = False
                  # If the special meijerg code did not succeed in
                  # finding a definite integral, then the code using
                  # meijerint_indefinite will not either (it might
                  # find an antiderivative, but the answer is likely
                  # to be nonsensical). Thus if we are requested to
                  # only use Meijer G-function methods, we give up at
                  # this stage. Otherwise we just disable G-function
                  # methods.
                  if meijerg1 is False and meijerg is True:
                      antideriv = None
                  else:
                      antideriv = self._eval_integral(
                          function, xab[0], **eval_kwargs)
                      if antideriv is None and meijerg is True:
                          ret = try_meijerg(function, xab)
                          if ret is not None:
                              function = ret
                              continue
  
              if not isinstance(antideriv, Integral) and antideriv is not None:
                  for atan_term in antideriv.atoms(atan):
                      atan_arg = atan_term.args[0]
                      # Checking `atan_arg` to be linear combination of `tan` or `cot`
                      for tan_part in atan_arg.atoms(tan):
                          x1 = Dummy('x1')
                          tan_exp1 = atan_arg.subs(tan_part, x1)
                          # The coefficient of `tan` should be constant
                          coeff = tan_exp1.diff(x1)
                          if x1 not in coeff.free_symbols:
                              a = tan_part.args[0]
                              antideriv = antideriv.subs(atan_term, Add(atan_term,
                                  sign(coeff)*pi*floor((a-pi/2)/pi)))
                      for cot_part in atan_arg.atoms(cot):
                          x1 = Dummy('x1')
                          cot_exp1 = atan_arg.subs(cot_part, x1)
                          # The coefficient of `cot` should be constant
                          coeff = cot_exp1.diff(x1)
                          if x1 not in coeff.free_symbols:
                              a = cot_part.args[0]
                              antideriv = antideriv.subs(atan_term, Add(atan_term,
                                  sign(coeff)*pi*floor((a)/pi)))
  
              if antideriv is None:
                  undone_limits.append(xab)
                  function = self.func(*([function] + [xab])).factor()
                  factored_function = function.factor()
                  if not isinstance(factored_function, Integral):
                      function = factored_function
                  continue
              else:
                  if len(xab) == 1:
                      function = antideriv
                  else:
                      if len(xab) == 3:
                          x, a, b = xab
                      elif len(xab) == 2:
                          x, b = xab
                          a = None
                      else:
                          raise NotImplementedError
  
                      if deep:
                          if isinstance(a, Basic):
                              a = a.doit(**hints)
                          if isinstance(b, Basic):
                              b = b.doit(**hints)
  
                      if antideriv.is_Poly:
                          gens = list(antideriv.gens)
                          gens.remove(x)
  
                          antideriv = antideriv.as_expr()
  
                          function = antideriv._eval_interval(x, a, b)
                          function = Poly(function, *gens)
                      else:
                          def is_indef_int(g, x):
                              return (isinstance(g, Integral) and
                                      any(i == (x,) for i in g.limits))
  
                          def eval_factored(f, x, a, b):
                              # _eval_interval for integrals with
                              # (constant) factors
                              # a single indefinite integral is assumed
                              args = []
                              for g in Mul.make_args(f):
                                  if is_indef_int(g, x):
                                      args.append(g._eval_interval(x, a, b))
                                  else:
                                      args.append(g)
                              return Mul(*args)
  
                          integrals, others, piecewises = [], [], []
                          for f in Add.make_args(antideriv):
                              if any(is_indef_int(g, x)
                                     for g in Mul.make_args(f)):
                                  integrals.append(f)
                              elif any(isinstance(g, Piecewise)
                                       for g in Mul.make_args(f)):
                                  piecewises.append(piecewise_fold(f))
                              else:
                                  others.append(f)
                          uneval = Add(*[eval_factored(f, x, a, b)
                                         for f in integrals])
                          try:
                              evalued = Add(*others)._eval_interval(x, a, b)
                              evalued_pw = piecewise_fold(Add(*piecewises))._eval_interval(x, a, b)
                              function = uneval + evalued + evalued_pw
                          except NotImplementedError:
                              # This can happen if _eval_interval depends in a
                              # complicated way on limits that cannot be computed
                              undone_limits.append(xab)
                              function = self.func(*([function] + [xab]))
                              factored_function = function.factor()
                              if not isinstance(factored_function, Integral):
                                  function = factored_function
          return function
  
      def _eval_derivative(self, sym):
          """Evaluate the derivative of the current Integral object by
          differentiating under the integral sign [1], using the Fundamental
          Theorem of Calculus [2] when possible.
  
          Whenever an Integral is encountered that is equivalent to zero or
          has an integrand that is independent of the variable of integration
          those integrals are performed. All others are returned as Integral
          instances which can be resolved with doit() (provided they are integrable).
  
          References:
             [1] https://en.wikipedia.org/wiki/Differentiation_under_the_integral_sign
             [2] https://en.wikipedia.org/wiki/Fundamental_theorem_of_calculus
  
          Examples
          ========
  
          >>> from sympy import Integral
          >>> from sympy.abc import x, y
          >>> i = Integral(x + y, y, (y, 1, x))
          >>> i.diff(x)
          Integral(x + y, (y, x)) + Integral(1, y, (y, 1, x))
          >>> i.doit().diff(x) == i.diff(x).doit()
          True
          >>> i.diff(y)
          0
  
          The previous must be true since there is no y in the evaluated integral:
  
          >>> i.free_symbols
          {x}
          >>> i.doit()
          2*x**3/3 - x/2 - 1/6
  
          """
  
          # differentiate under the integral sign; we do not
          # check for regularity conditions (TODO), see issue 4215
  
          # get limits and the function
          f, limits = self.function, list(self.limits)
  
          # the order matters if variables of integration appear in the limits
          # so work our way in from the outside to the inside.
          limit = limits.pop(-1)
          if len(limit) == 3:
              x, a, b = limit
          elif len(limit) == 2:
              x, b = limit
              a = None
          else:
              a = b = None
              x = limit[0]
  
          if limits:  # f is the argument to an integral
              f = self.func(f, *tuple(limits))
  
          # assemble the pieces
          def _do(f, ab):
              dab_dsym = diff(ab, sym)
              if not dab_dsym:
                  return S.Zero
              if isinstance(f, Integral):
                  limits = [(x, x) if (len(l) == 1 and l[0] == x) else l
                            for l in f.limits]
                  f = self.func(f.function, *limits)
              return f.subs(x, ab)*dab_dsym
  
          rv = S.Zero
          if b is not None:
              rv += _do(f, b)
          if a is not None:
              rv -= _do(f, a)
          if len(limit) == 1 and sym == x:
              # the dummy variable *is* also the real-world variable
              arg = f
              rv += arg
          else:
              # the dummy variable might match sym but it's
              # only a dummy and the actual variable is determined
              # by the limits, so mask off the variable of integration
              # while differentiating
              u = Dummy('u')
              arg = f.subs(x, u).diff(sym).subs(u, x)
              if arg:
                  rv += self.func(arg, Tuple(x, a, b))
          return rv
  
      def _eval_integral(self, f, x, meijerg=None, risch=None, manual=None,
                         heurisch=None, conds='piecewise'):
          """
          Calculate the anti-derivative to the function f(x).
  
          The following algorithms are applied (roughly in this order):
  
          1. Simple heuristics (based on pattern matching and integral table):
  
             - most frequently used functions (e.g. polynomials, products of
               trig functions)
  
          2. Integration of rational functions:
  
             - A complete algorithm for integrating rational functions is
               implemented (the Lazard-Rioboo-Trager algorithm).  The algorithm
               also uses the partial fraction decomposition algorithm
               implemented in apart() as a preprocessor to make this process
               faster.  Note that the integral of a rational function is always
               elementary, but in general, it may include a RootSum.
  
          3. Full Risch algorithm:
  
             - The Risch algorithm is a complete decision
               procedure for integrating elementary functions, which means that
               given any elementary function, it will either compute an
               elementary antiderivative, or else prove that none exists.
               Currently, part of transcendental case is implemented, meaning
               elementary integrals containing exponentials, logarithms, and
               (soon!) trigonometric functions can be computed.  The algebraic
               case, e.g., functions containing roots, is much more difficult
               and is not implemented yet.
  
             - If the routine fails (because the integrand is not elementary, or
               because a case is not implemented yet), it continues on to the
               next algorithms below.  If the routine proves that the integrals
               is nonelementary, it still moves on to the algorithms below,
               because we might be able to find a closed-form solution in terms
               of special functions.  If risch=True, however, it will stop here.
  
          4. The Meijer G-Function algorithm:
  
             - This algorithm works by first rewriting the integrand in terms of
               very general Meijer G-Function (meijerg in SymPy), integrating
               it, and then rewriting the result back, if possible.  This
               algorithm is particularly powerful for definite integrals (which
               is actually part of a different method of Integral), since it can
               compute closed-form solutions of definite integrals even when no
               closed-form indefinite integral exists.  But it also is capable
               of computing many indefinite integrals as well.
  
             - Another advantage of this method is that it can use some results
               about the Meijer G-Function to give a result in terms of a
               Piecewise expression, which allows to express conditionally
               convergent integrals.
  
             - Setting meijerg=True will cause integrate() to use only this
               method.
  
          5. The "manual integration" algorithm:
  
             - This algorithm tries to mimic how a person would find an
               antiderivative by hand, for example by looking for a
               substitution or applying integration by parts. This algorithm
               does not handle as many integrands but can return results in a
               more familiar form.
  
             - Sometimes this algorithm can evaluate parts of an integral; in
               this case integrate() will try to evaluate the rest of the
               integrand using the other methods here.
  
             - Setting manual=True will cause integrate() to use only this
               method.
  
          6. The Heuristic Risch algorithm:
  
             - This is a heuristic version of the Risch algorithm, meaning that
               it is not deterministic.  This is tried as a last resort because
               it can be very slow.  It is still used because not enough of the
               full Risch algorithm is implemented, so that there are still some
               integrals that can only be computed using this method.  The goal
               is to implement enough of the Risch and Meijer G-function methods
               so that this can be deleted.
  
               Setting heurisch=True will cause integrate() to use only this
               method. Set heurisch=False to not use it.
  
          """
          from sympy.integrals.deltafunctions import deltaintegrate
          from sympy.integrals.singularityfunctions import singularityintegrate
          from sympy.integrals.heurisch import heurisch as heurisch_, heurisch_wrapper
          from sympy.integrals.rationaltools import ratint
          from sympy.integrals.risch import risch_integrate
  
          if risch:
              try:
                  return risch_integrate(f, x, conds=conds)
              except NotImplementedError:
                  return None
  
          if manual:
              try:
                  result = manualintegrate(f, x)
                  if result is not None and result.func != Integral:
                      return result
              except (ValueError, PolynomialError):
                  pass
  
          eval_kwargs = dict(meijerg=meijerg, risch=risch, manual=manual,
              heurisch=heurisch, conds=conds)
  
          # if it is a poly(x) then let the polynomial integrate itself (fast)
          #
          # It is important to make this check first, otherwise the other code
          # will return a sympy expression instead of a Polynomial.
          #
          # see Polynomial for details.
          if isinstance(f, Poly) and not (manual or meijerg or risch):
              SymPyDeprecationWarning(
                  feature="Using integrate/Integral with Poly",
                  issue=18613,
                  deprecated_since_version="1.6",
                  useinstead="the as_expr or integrate methods of Poly").warn()
              return f.integrate(x)
  
          # Piecewise antiderivatives need to call special integrate.
          if isinstance(f, Piecewise):
              return f.piecewise_integrate(x, **eval_kwargs)
  
          # let's cut it short if `f` does not depend on `x`; if
          # x is only a dummy, that will be handled below
          if not f.has(x):
              return f*x
  
          # try to convert to poly(x) and then integrate if successful (fast)
          poly = f.as_poly(x)
          if poly is not None and not (manual or meijerg or risch):
              return poly.integrate().as_expr()
  
          if risch is not False:
              try:
                  result, i = risch_integrate(f, x, separate_integral=True,
                      conds=conds)
              except NotImplementedError:
                  pass
              else:
                  if i:
                      # There was a nonelementary integral. Try integrating it.
  
                      # if no part of the NonElementaryIntegral is integrated by
                      # the Risch algorithm, then use the original function to
                      # integrate, instead of re-written one
                      if result == 0:
                          from sympy.integrals.risch import NonElementaryIntegral
                          return NonElementaryIntegral(f, x).doit(risch=False)
                      else:
                          return result + i.doit(risch=False)
                  else:
                      return result
  
          # since Integral(f=g1+g2+...) == Integral(g1) + Integral(g2) + ...
          # we are going to handle Add terms separately,
          # if `f` is not Add -- we only have one term
  
          # Note that in general, this is a bad idea, because Integral(g1) +
          # Integral(g2) might not be computable, even if Integral(g1 + g2) is.
          # For example, Integral(x**x + x**x*log(x)).  But many heuristics only
          # work term-wise.  So we compute this step last, after trying
          # risch_integrate.  We also try risch_integrate again in this loop,
          # because maybe the integral is a sum of an elementary part and a
          # nonelementary part (like erf(x) + exp(x)).  risch_integrate() is
          # quite fast, so this is acceptable.
          parts = []
          args = Add.make_args(f)
          for g in args:
              coeff, g = g.as_independent(x)
  
              # g(x) = const
              if g is S.One and not meijerg:
                  parts.append(coeff*x)
                  continue
  
              # g(x) = expr + O(x**n)
              order_term = g.getO()
  
              if order_term is not None:
                  h = self._eval_integral(g.removeO(), x, **eval_kwargs)
  
                  if h is not None:
                      h_order_expr = self._eval_integral(order_term.expr, x, **eval_kwargs)
  
                      if h_order_expr is not None:
                          h_order_term = order_term.func(
                              h_order_expr, *order_term.variables)
                          parts.append(coeff*(h + h_order_term))
                          continue
  
                  # NOTE: if there is O(x**n) and we fail to integrate then
                  # there is no point in trying other methods because they
                  # will fail, too.
                  return None
  
              #               c
              # g(x) = (a*x+b)
              if g.is_Pow and not g.exp.has(x) and not meijerg:
                  a = Wild('a', exclude=[x])
                  b = Wild('b', exclude=[x])
  
                  M = g.base.match(a*x + b)
  
                  if M is not None:
                      if g.exp == -1:
                          h = log(g.base)
                      elif conds != 'piecewise':
                          h = g.base**(g.exp + 1) / (g.exp + 1)
                      else:
                          h1 = log(g.base)
                          h2 = g.base**(g.exp + 1) / (g.exp + 1)
                          h = Piecewise((h2, Ne(g.exp, -1)), (h1, True))
  
                      parts.append(coeff * h / M[a])
                      continue
  
              #        poly(x)
              # g(x) = -------
              #        poly(x)
              if g.is_rational_function(x) and not (manual or meijerg or risch):
                  parts.append(coeff * ratint(g, x))
                  continue
  
              if not (manual or meijerg or risch):
                  # g(x) = Mul(trig)
                  h = trigintegrate(g, x, conds=conds)
                  if h is not None:
                      parts.append(coeff * h)
                      continue
  
                  # g(x) has at least a DiracDelta term
                  h = deltaintegrate(g, x)
                  if h is not None:
                      parts.append(coeff * h)
                      continue
  
                  # g(x) has at least a Singularity Function term
                  h = singularityintegrate(g, x)
                  if h is not None:
                      parts.append(coeff * h)
                      continue
  
                  # Try risch again.
                  if risch is not False:
                      try:
                          h, i = risch_integrate(g, x,
                              separate_integral=True, conds=conds)
                      except NotImplementedError:
                          h = None
                      else:
                          if i:
                              h = h + i.doit(risch=False)
  
                          parts.append(coeff*h)
                          continue
  
                  # fall back to heurisch
                  if heurisch is not False:
                      try:
                          if conds == 'piecewise':
                              h = heurisch_wrapper(g, x, hints=[])
                          else:
                              h = heurisch_(g, x, hints=[])
                      except PolynomialError:
                          # XXX: this exception means there is a bug in the
                          # implementation of heuristic Risch integration
                          # algorithm.
                          h = None
              else:
                  h = None
  
              if meijerg is not False and h is None:
                  # rewrite using G functions
                  try:
                      h = meijerint_indefinite(g, x)
                  except NotImplementedError:
                      from sympy.integrals.meijerint import _debug
                      _debug('NotImplementedError from meijerint_definite')
                  if h is not None:
                      parts.append(coeff * h)
                      continue
  
              if h is None and manual is not False:
                  try:
                      result = manualintegrate(g, x)
                      if result is not None and not isinstance(result, Integral):
                          if result.has(Integral) and not manual:
                              # Try to have other algorithms do the integrals
                              # manualintegrate can't handle,
                              # unless we were asked to use manual only.
                              # Keep the rest of eval_kwargs in case another
                              # method was set to False already
                              new_eval_kwargs = eval_kwargs
                              new_eval_kwargs["manual"] = False
                              result = result.func(*[
                                  arg.doit(**new_eval_kwargs) if
                                  arg.has(Integral) else arg
                                  for arg in result.args
                              ]).expand(multinomial=False,
                                        log=False,
                                        power_exp=False,
                                        power_base=False)
                          if not result.has(Integral):
                              parts.append(coeff * result)
                              continue
                  except (ValueError, PolynomialError):
                      # can't handle some SymPy expressions
                      pass
  
              # if we failed maybe it was because we had
              # a product that could have been expanded,
              # so let's try an expansion of the whole
              # thing before giving up; we don't try this
              # at the outset because there are things
              # that cannot be solved unless they are
              # NOT expanded e.g., x**x*(1+log(x)). There
              # should probably be a checker somewhere in this
              # routine to look for such cases and try to do
              # collection on the expressions if they are already
              # in an expanded form
              if not h and len(args) == 1:
                  f = sincos_to_sum(f).expand(mul=True, deep=False)
                  if f.is_Add:
                      # Note: risch will be identical on the expanded
                      # expression, but maybe it will be able to pick out parts,
                      # like x*(exp(x) + erf(x)).
                      return self._eval_integral(f, x, **eval_kwargs)
  
              if h is not None:
                  parts.append(coeff * h)
              else:
                  return None
  
          return Add(*parts)
  
      def _eval_lseries(self, x, logx):
          expr = self.as_dummy()
          symb = x
          for l in expr.limits:
              if x in l[1:]:
                  symb = l[0]
                  break
          for term in expr.function.lseries(symb, logx):
              yield integrate(term, *expr.limits)
  
      def _eval_nseries(self, x, n, logx):
          expr = self.as_dummy()
          symb = x
          for l in expr.limits:
              if x in l[1:]:
                  symb = l[0]
                  break
          terms, order = expr.function.nseries(
              x=symb, n=n, logx=logx).as_coeff_add(Order)
          order = [o.subs(symb, x) for o in order]
          return integrate(terms, *expr.limits) + Add(*order)*x
  
      def _eval_as_leading_term(self, x):
          series_gen = self.args[0].lseries(x)
          for leading_term in series_gen:
              if leading_term != 0:
                  break
          return integrate(leading_term, *self.args[1:])
  
      def _eval_simplify(self, **kwargs):
          from sympy.core.exprtools import factor_terms
          from sympy.simplify.simplify import simplify
  
          expr = factor_terms(self)
          if isinstance(expr, Integral):
              return expr.func(*[simplify(i, **kwargs) for i in expr.args])
          return expr.simplify(**kwargs)
  
      def as_sum(self, n=None, method="midpoint", evaluate=True):
          """
          Approximates a definite integral by a sum.
  
          Arguments
          ---------
          n
              The number of subintervals to use, optional.
          method
              One of: 'left', 'right', 'midpoint', 'trapezoid'.
          evaluate
              If False, returns an unevaluated Sum expression. The default
              is True, evaluate the sum.
  
          These methods of approximate integration are described in [1].
  
          [1] https://en.wikipedia.org/wiki/Riemann_sum#Methods
  
          Examples
          ========
  
          >>> from sympy import sin, sqrt
          >>> from sympy.abc import x, n
          >>> from sympy.integrals import Integral
          >>> e = Integral(sin(x), (x, 3, 7))
          >>> e
          Integral(sin(x), (x, 3, 7))
  
          For demonstration purposes, this interval will only be split into 2
          regions, bounded by [3, 5] and [5, 7].
  
          The left-hand rule uses function evaluations at the left of each
          interval:
  
          >>> e.as_sum(2, 'left')
          2*sin(5) + 2*sin(3)
  
          The midpoint rule uses evaluations at the center of each interval:
  
          >>> e.as_sum(2, 'midpoint')
          2*sin(4) + 2*sin(6)
  
          The right-hand rule uses function evaluations at the right of each
          interval:
  
          >>> e.as_sum(2, 'right')
          2*sin(5) + 2*sin(7)
  
          The trapezoid rule uses function evaluations on both sides of the
          intervals. This is equivalent to taking the average of the left and
          right hand rule results:
  
          >>> e.as_sum(2, 'trapezoid')
          2*sin(5) + sin(3) + sin(7)
          >>> (e.as_sum(2, 'left') + e.as_sum(2, 'right'))/2 == _
          True
  
          Here, the discontinuity at x = 0 can be avoided by using the
          midpoint or right-hand method:
  
          >>> e = Integral(1/sqrt(x), (x, 0, 1))
          >>> e.as_sum(5).n(4)
          1.730
          >>> e.as_sum(10).n(4)
          1.809
          >>> e.doit().n(4)  # the actual value is 2
          2.000
  
          The left- or trapezoid method will encounter the discontinuity and
          return infinity:
  
          >>> e.as_sum(5, 'left')
          zoo
  
          The number of intervals can be symbolic. If omitted, a dummy symbol
          will be used for it.
  
          >>> e = Integral(x**2, (x, 0, 2))
          >>> e.as_sum(n, 'right').expand()
          8/3 + 4/n + 4/(3*n**2)
  
          This shows that the midpoint rule is more accurate, as its error
          term decays as the square of n:
  
          >>> e.as_sum(method='midpoint').expand()
          8/3 - 2/(3*_n**2)
  
          A symbolic sum is returned with evaluate=False:
  
          >>> e.as_sum(n, 'midpoint', evaluate=False)
          2*Sum((2*_k/n - 1/n)**2, (_k, 1, n))/n
  
          See Also
          ========
  
          Integral.doit : Perform the integration using any hints
          """
  
          from sympy.concrete.summations import Sum
          limits = self.limits
          if len(limits) > 1:
              raise NotImplementedError(
                  "Multidimensional midpoint rule not implemented yet")
          else:
              limit = limits[0]
              if (len(limit) != 3 or limit[1].is_finite is False or
                  limit[2].is_finite is False):
                  raise ValueError("Expecting a definite integral over "
                                    "a finite interval.")
          if n is None:
              n = Dummy('n', integer=True, positive=True)
          else:
              n = sympify(n)
          if (n.is_positive is False or n.is_integer is False or
              n.is_finite is False):
              raise ValueError("n must be a positive integer, got %s" % n)
          x, a, b = limit
          dx = (b - a)/n
          k = Dummy('k', integer=True, positive=True)
          f = self.function
  
          if method == "left":
              result = dx*Sum(f.subs(x, a + (k-1)*dx), (k, 1, n))
          elif method == "right":
              result = dx*Sum(f.subs(x, a + k*dx), (k, 1, n))
          elif method == "midpoint":
              result = dx*Sum(f.subs(x, a + k*dx - dx/2), (k, 1, n))
          elif method == "trapezoid":
              result = dx*((f.subs(x, a) + f.subs(x, b))/2 +
                  Sum(f.subs(x, a + k*dx), (k, 1, n - 1)))
          else:
              raise ValueError("Unknown method %s" % method)
          return result.doit() if evaluate else result
  
      def _sage_(self):
          import sage.all as sage
          f, limits = self.function._sage_(), list(self.limits)
          for limit_ in limits:
              if len(limit_) == 1:
                  x = limit_[0]
                  f = sage.integral(f,
                                      x._sage_(),
                                      hold=True)
              elif len(limit_) == 2:
                  x, b = limit_
                  f = sage.integral(f,
                                      x._sage_(),
                                      b._sage_(),
                                      hold=True)
              else:
                  x, a, b = limit_
                  f = sage.integral(f,
                                    (x._sage_(),
                                      a._sage_(),
                                      b._sage_()),
                                      hold=True)
          return f
  
      def principal_value(self, **kwargs):
          """
          Compute the Cauchy Principal Value of the definite integral of a real function in the given interval
          on the real axis.
          In mathematics, the Cauchy principal value, is a method for assigning values to certain improper
          integrals which would otherwise be undefined.
  
          Examples
          ========
  
          >>> from sympy import oo
          >>> from sympy.integrals.integrals import Integral
          >>> from sympy.abc import x
          >>> Integral(x+1, (x, -oo, oo)).principal_value()
          oo
          >>> f = 1 / (x**3)
          >>> Integral(f, (x, -oo, oo)).principal_value()
          0
          >>> Integral(f, (x, -10, 10)).principal_value()
          0
          >>> Integral(f, (x, -10, oo)).principal_value() + Integral(f, (x, -oo, 10)).principal_value()
          0
  
          References
          ==========
          .. [1] https://en.wikipedia.org/wiki/Cauchy_principal_value
          .. [2] http://mathworld.wolfram.com/CauchyPrincipalValue.html
          """
          from sympy.calculus import singularities
          if len(self.limits) != 1 or len(list(self.limits[0])) != 3:
              raise ValueError("You need to insert a variable, lower_limit, and upper_limit correctly to calculate "
                               "cauchy's principal value")
          x, a, b = self.limits[0]
          if not (a.is_comparable and b.is_comparable and a <= b):
              raise ValueError("The lower_limit must be smaller than or equal to the upper_limit to calculate "
                               "cauchy's principal value. Also, a and b need to be comparable.")
          if a == b:
              return 0
          r = Dummy('r')
          f = self.function
          singularities_list = [s for s in singularities(f, x) if s.is_comparable and a <= s <= b]
          for i in singularities_list:
              if (i == b) or (i == a):
                  raise ValueError(
                      'The principal value is not defined in the given interval due to singularity at %d.' % (i))
          F = integrate(f, x, **kwargs)
          if F.has(Integral):
              return self
          if a is -oo and b is oo:
              I = limit(F - F.subs(x, -x), x, oo)
          else:
              I = limit(F, x, b, '-') - limit(F, x, a, '+')
          for s in singularities_list:
              I += limit(((F.subs(x, s - r)) - F.subs(x, s + r)), r, 0, '+')
          return I
  
  
  
  def integrate(*args, **kwargs):
      """integrate(f, var, ...)
  
      Compute definite or indefinite integral of one or more variables
      using Risch-Norman algorithm and table lookup. This procedure is
      able to handle elementary algebraic and transcendental functions
      and also a huge class of special functions, including Airy,
      Bessel, Whittaker and Lambert.
  
      var can be:
  
      - a symbol                   -- indefinite integration
      - a tuple (symbol, a)        -- indefinite integration with result
                                      given with `a` replacing `symbol`
      - a tuple (symbol, a, b)     -- definite integration
  
      Several variables can be specified, in which case the result is
      multiple integration. (If var is omitted and the integrand is
      univariate, the indefinite integral in that variable will be performed.)
  
      Indefinite integrals are returned without terms that are independent
      of the integration variables. (see examples)
  
      Definite improper integrals often entail delicate convergence
      conditions. Pass conds='piecewise', 'separate' or 'none' to have
      these returned, respectively, as a Piecewise function, as a separate
      result (i.e. result will be a tuple), or not at all (default is
      'piecewise').
  
      **Strategy**
  
      SymPy uses various approaches to definite integration. One method is to
      find an antiderivative for the integrand, and then use the fundamental
      theorem of calculus. Various functions are implemented to integrate
      polynomial, rational and trigonometric functions, and integrands
      containing DiracDelta terms.
  
      SymPy also implements the part of the Risch algorithm, which is a decision
      procedure for integrating elementary functions, i.e., the algorithm can
      either find an elementary antiderivative, or prove that one does not
      exist.  There is also a (very successful, albeit somewhat slow) general
      implementation of the heuristic Risch algorithm.  This algorithm will
      eventually be phased out as more of the full Risch algorithm is
      implemented. See the docstring of Integral._eval_integral() for more
      details on computing the antiderivative using algebraic methods.
  
      The option risch=True can be used to use only the (full) Risch algorithm.
      This is useful if you want to know if an elementary function has an
      elementary antiderivative.  If the indefinite Integral returned by this
      function is an instance of NonElementaryIntegral, that means that the
      Risch algorithm has proven that integral to be non-elementary.  Note that
      by default, additional methods (such as the Meijer G method outlined
      below) are tried on these integrals, as they may be expressible in terms
      of special functions, so if you only care about elementary answers, use
      risch=True.  Also note that an unevaluated Integral returned by this
      function is not necessarily a NonElementaryIntegral, even with risch=True,
      as it may just be an indication that the particular part of the Risch
      algorithm needed to integrate that function is not yet implemented.
  
      Another family of strategies comes from re-writing the integrand in
      terms of so-called Meijer G-functions. Indefinite integrals of a
      single G-function can always be computed, and the definite integral
      of a product of two G-functions can be computed from zero to
      infinity. Various strategies are implemented to rewrite integrands
      as G-functions, and use this information to compute integrals (see
      the ``meijerint`` module).
  
      The option manual=True can be used to use only an algorithm that tries
      to mimic integration by hand. This algorithm does not handle as many
      integrands as the other algorithms implemented but may return results in
      a more familiar form. The ``manualintegrate`` module has functions that
      return the steps used (see the module docstring for more information).
  
      In general, the algebraic methods work best for computing
      antiderivatives of (possibly complicated) combinations of elementary
      functions. The G-function methods work best for computing definite
      integrals from zero to infinity of moderately complicated
      combinations of special functions, or indefinite integrals of very
      simple combinations of special functions.
  
      The strategy employed by the integration code is as follows:
  
      - If computing a definite integral, and both limits are real,
        and at least one limit is +- oo, try the G-function method of
        definite integration first.
  
      - Try to find an antiderivative, using all available methods, ordered
        by performance (that is try fastest method first, slowest last; in
        particular polynomial integration is tried first, Meijer
        G-functions second to last, and heuristic Risch last).
  
      - If still not successful, try G-functions irrespective of the
        limits.
  
      The option meijerg=True, False, None can be used to, respectively:
      always use G-function methods and no others, never use G-function
      methods, or use all available methods (in order as described above).
      It defaults to None.
  
      Examples
      ========
  
      >>> from sympy import integrate, log, exp, oo
      >>> from sympy.abc import a, x, y
  
      >>> integrate(x*y, x)
      x**2*y/2
  
      >>> integrate(log(x), x)
      x*log(x) - x
  
      >>> integrate(log(x), (x, 1, a))
      a*log(a) - a + 1
  
      >>> integrate(x)
      x**2/2
  
      Terms that are independent of x are dropped by indefinite integration:
  
      >>> from sympy import sqrt
      >>> integrate(sqrt(1 + x), (x, 0, x))
      2*(x + 1)**(3/2)/3 - 2/3
      >>> integrate(sqrt(1 + x), x)
      2*(x + 1)**(3/2)/3
  
      >>> integrate(x*y)
      Traceback (most recent call last):
      ...
      ValueError: specify integration variables to integrate x*y
  
      Note that ``integrate(x)`` syntax is meant only for convenience
      in interactive sessions and should be avoided in library code.
  
      >>> integrate(x**a*exp(-x), (x, 0, oo)) # same as conds='piecewise'
      Piecewise((gamma(a + 1), re(a) > -1),
          (Integral(x**a*exp(-x), (x, 0, oo)), True))
  
      >>> integrate(x**a*exp(-x), (x, 0, oo), conds='none')
      gamma(a + 1)
  
      >>> integrate(x**a*exp(-x), (x, 0, oo), conds='separate')
      (gamma(a + 1), -re(a) < 1)
  
      See Also
      ========
  
      Integral, Integral.doit
  
      """
      doit_flags = {
          'deep': False,
          'meijerg': kwargs.pop('meijerg', None),
          'conds': kwargs.pop('conds', 'piecewise'),
          'risch': kwargs.pop('risch', None),
          'heurisch': kwargs.pop('heurisch', None),
          'manual': kwargs.pop('manual', None)
          }
      integral = Integral(*args, **kwargs)
  
      if isinstance(integral, Integral):
          return integral.doit(**doit_flags)
      else:
          new_args = [a.doit(**doit_flags) if isinstance(a, Integral) else a
              for a in integral.args]
          return integral.func(*new_args)
  
  
  def line_integrate(field, curve, vars):
      """line_integrate(field, Curve, variables)
  
      Compute the line integral.
  
      Examples
      ========
  
      >>> from sympy import Curve, line_integrate, E, ln
      >>> from sympy.abc import x, y, t
      >>> C = Curve([E**t + 1, E**t - 1], (t, 0, ln(2)))
      >>> line_integrate(x + y, C, [x, y])
      3*sqrt(2)
  
      See Also
      ========
  
      sympy.integrals.integrals.integrate, Integral
      """
      from sympy.geometry import Curve
      F = sympify(field)
      if not F:
          raise ValueError(
              "Expecting function specifying field as first argument.")
      if not isinstance(curve, Curve):
          raise ValueError("Expecting Curve entity as second argument.")
      if not is_sequence(vars):
          raise ValueError("Expecting ordered iterable for variables.")
      if len(curve.functions) != len(vars):
          raise ValueError("Field variable size does not match curve dimension.")
  
      if curve.parameter in vars:
          raise ValueError("Curve parameter clashes with field parameters.")
  
      # Calculate derivatives for line parameter functions
      # F(r) -> F(r(t)) and finally F(r(t)*r'(t))
      Ft = F
      dldt = 0
      for i, var in enumerate(vars):
          _f = curve.functions[i]
          _dn = diff(_f, curve.parameter)
          # ...arc length
          dldt = dldt + (_dn * _dn)
          Ft = Ft.subs(var, _f)
      Ft = Ft * sqrt(dldt)
  
      integral = Integral(Ft, curve.limits).doit(deep=False)
      return integral
  