image-match/python/py/Lib/site-packages/scipy/integrate/_quadrature.py

import numpy as np
import math
import warnings
from collections import namedtuple
from collections.abc import Callable

from scipy.special import roots_legendre
from scipy.special import gammaln, logsumexp
from scipy._lib._util import _rng_spawn
from scipy._lib._array_api import (_asarray, array_namespace, xp_result_type, xp_copy,
                                   xp_capabilities, xp_promote, xp_swapaxes, is_numpy)
import scipy._lib.array_api_extra as xpx


__all__ = ['fixed_quad', 'romb',
           'trapezoid', 'simpson',
           'cumulative_trapezoid', 'newton_cotes',
           'qmc_quad', 'cumulative_simpson']


@xp_capabilities(skip_backends=[('jax.numpy',
                                 "JAX arrays do not support item assignment")])
def trapezoid(y, x=None, dx=1.0, axis=-1):
    r"""
    Integrate along the given axis using the composite trapezoidal rule.

    If `x` is provided, the integration happens in sequence along its
    elements - they are not sorted.

    Integrate `y` (`x`) along each 1d slice on the given axis, compute
    :math:`\int y(x) dx`.
    When `x` is specified, this integrates along the parametric curve,
    computing :math:`\int_t y(t) dt =
    \int_t y(t) \left.\frac{dx}{dt}\right|_{x=x(t)} dt`.

    Parameters
    ----------
    y : array_like
        Input array to integrate.
    x : array_like, optional
        The sample points corresponding to the `y` values. If `x` is None,
        the sample points are assumed to be evenly spaced `dx` apart. The
        default is None.
    dx : scalar, optional
        The spacing between sample points when `x` is None. The default is 1.
    axis : int, optional
        The axis along which to integrate. The default is the last axis.

    Returns
    -------
    trapezoid : float or ndarray
        Definite integral of `y` = n-dimensional array as approximated along
        a single axis by the trapezoidal rule. If `y` is a 1-dimensional array,
        then the result is a float. If `n` is greater than 1, then the result
        is an `n`-1 dimensional array.

    See Also
    --------
    cumulative_trapezoid, simpson, romb

    Notes
    -----
    Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
    will be taken from `y` array, by default x-axis distances between
    points will be 1.0, alternatively they can be provided with `x` array
    or with `dx` scalar.  Return value will be equal to combined area under
    the red lines.

    References
    ----------
    .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule

    .. [2] Illustration image:
           https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png

    Examples
    --------
    Use the trapezoidal rule on evenly spaced points:

    >>> import numpy as np
    >>> from scipy import integrate
    >>> integrate.trapezoid([1, 2, 3])
    4.0

    The spacing between sample points can be selected by either the
    ``x`` or ``dx`` arguments:

    >>> integrate.trapezoid([1, 2, 3], x=[4, 6, 8])
    8.0
    >>> integrate.trapezoid([1, 2, 3], dx=2)
    8.0

    Using a decreasing ``x`` corresponds to integrating in reverse:

    >>> integrate.trapezoid([1, 2, 3], x=[8, 6, 4])
    -8.0

    More generally ``x`` is used to integrate along a parametric curve. We can
    estimate the integral :math:`\int_0^1 x^2 = 1/3` using:

    >>> x = np.linspace(0, 1, num=50)
    >>> y = x**2
    >>> integrate.trapezoid(y, x)
    0.33340274885464394

    Or estimate the area of a circle, noting we repeat the sample which closes
    the curve:

    >>> theta = np.linspace(0, 2 * np.pi, num=1000, endpoint=True)
    >>> integrate.trapezoid(np.cos(theta), x=np.sin(theta))
    3.141571941375841

    ``trapezoid`` can be applied along a specified axis to do multiple
    computations in one call:

    >>> a = np.arange(6).reshape(2, 3)
    >>> a
    array([[0, 1, 2],
           [3, 4, 5]])
    >>> integrate.trapezoid(a, axis=0)
    array([1.5, 2.5, 3.5])
    >>> integrate.trapezoid(a, axis=1)
    array([2.,  8.])
    """
    xp = array_namespace(y)
    y = _asarray(y, xp=xp, subok=True)
    # Cannot just use the broadcasted arrays that are returned
    # because trapezoid does not follow normal broadcasting rules
    # cf. https://github.com/scipy/scipy/pull/21524#issuecomment-2354105942
    result_dtype = xp_result_type(y, force_floating=True, xp=xp)
    nd = y.ndim
    slice1 = [slice(None)]*nd
    slice2 = [slice(None)]*nd
    slice1[axis] = slice(1, None)
    slice2[axis] = slice(None, -1)
    if x is None:
        d = dx
    else:
        x = _asarray(x, xp=xp, subok=True)
        if x.ndim == 1:
            d = x[1:] - x[:-1]
            # make d broadcastable to y
            slice3 = [None] * nd
            slice3[axis] = slice(None)
            d = d[tuple(slice3)]
        else:
            # if x is n-D it should be broadcastable to y
            x = xp.broadcast_to(x, y.shape)
            d = x[tuple(slice1)] - x[tuple(slice2)]
    try:
        ret = xp.sum(
            d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0,
            axis=axis, dtype=result_dtype
        )
    except ValueError:
        # Operations didn't work, cast to ndarray
        d = xp.asarray(d)
        y = xp.asarray(y)
        ret = xp.sum(
            d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0,
            axis=axis, dtype=result_dtype
        )
    return ret


def _cached_roots_legendre(n):
    """
    Cache roots_legendre results to speed up calls of the fixed_quad
    function.
    """
    if n in _cached_roots_legendre.cache:
        return _cached_roots_legendre.cache[n]

    _cached_roots_legendre.cache[n] = roots_legendre(n)
    return _cached_roots_legendre.cache[n]


_cached_roots_legendre.cache = dict()


@xp_capabilities(np_only=True)
def fixed_quad(func, a, b, args=(), n=5):
    """
    Compute a definite integral using fixed-order Gaussian quadrature.

    Integrate `func` from `a` to `b` using Gaussian quadrature of
    order `n`.

    Parameters
    ----------
    func : callable
        A Python function or method to integrate (must accept vector inputs).
        If integrating a vector-valued function, the returned array must have
        shape ``(..., len(x))``.
    a : float
        Lower limit of integration.
    b : float
        Upper limit of integration.
    args : tuple, optional
        Extra arguments to pass to function, if any.
    n : int, optional
        Order of quadrature integration. Default is 5.

    Returns
    -------
    val : float
        Gaussian quadrature approximation to the integral
    none : None
        Statically returned value of None

    See Also
    --------
    quad : adaptive quadrature using QUADPACK
    dblquad : double integrals
    tplquad : triple integrals
    romb : integrators for sampled data
    simpson : integrators for sampled data
    cumulative_trapezoid : cumulative integration for sampled data

    Examples
    --------
    >>> from scipy import integrate
    >>> import numpy as np
    >>> f = lambda x: x**8
    >>> integrate.fixed_quad(f, 0.0, 1.0, n=4)
    (0.1110884353741496, None)
    >>> integrate.fixed_quad(f, 0.0, 1.0, n=5)
    (0.11111111111111102, None)
    >>> print(1/9.0)  # analytical result
    0.1111111111111111

    >>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=4)
    (0.9999999771971152, None)
    >>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=5)
    (1.000000000039565, None)
    >>> np.sin(np.pi/2)-np.sin(0)  # analytical result
    1.0

    """
    x, w = _cached_roots_legendre(n)
    x = np.real(x)
    if np.isinf(a) or np.isinf(b):
        raise ValueError("Gaussian quadrature is only available for "
                         "finite limits.")
    y = (b-a)*(x+1)/2.0 + a
    return (b-a)/2.0 * np.sum(w*func(y, *args), axis=-1), None


def tupleset(t, i, value):
    l = list(t)
    l[i] = value
    return tuple(l)


@xp_capabilities()
def cumulative_trapezoid(y, x=None, dx=1.0, axis=-1, initial=None):
    """
    Cumulatively integrate y(x) using the composite trapezoidal rule.

    Parameters
    ----------
    y : array_like
        Values to integrate.
    x : array_like, optional
        The coordinate to integrate along. If None (default), use spacing `dx`
        between consecutive elements in `y`.
    dx : float, optional
        Spacing between elements of `y`. Only used if `x` is None.
    axis : int, optional
        Specifies the axis to cumulate. Default is -1 (last axis).
    initial : scalar, optional
        If given, insert this value at the beginning of the returned result.
        0 or None are the only values accepted. Default is None, which means
        `res` has one element less than `y` along the axis of integration.

    Returns
    -------
    res : ndarray
        The result of cumulative integration of `y` along `axis`.
        If `initial` is None, the shape is such that the axis of integration
        has one less value than `y`. If `initial` is given, the shape is equal
        to that of `y`.

    See Also
    --------
    numpy.cumsum, numpy.cumprod
    cumulative_simpson : cumulative integration using Simpson's 1/3 rule
    quad : adaptive quadrature using QUADPACK
    fixed_quad : fixed-order Gaussian quadrature
    dblquad : double integrals
    tplquad : triple integrals
    romb : integrators for sampled data

    Examples
    --------
    >>> from scipy import integrate
    >>> import numpy as np
    >>> import matplotlib.pyplot as plt

    >>> x = np.linspace(-2, 2, num=20)
    >>> y = x
    >>> y_int = integrate.cumulative_trapezoid(y, x, initial=0)
    >>> plt.plot(x, y_int, 'ro', x, y[0] + 0.5 * x**2, 'b-')
    >>> plt.show()

    """
    xp = array_namespace(y)
    y = xp.asarray(y)

    if y.shape[axis] == 0:
        raise ValueError("At least one point is required along `axis`.")

    if x is None:
        d = dx
    else:
        x = xp.asarray(x)
        if x.ndim == 1:
            d = xp.diff(x)
            # reshape to correct shape
            shape = [1] * y.ndim
            shape[axis] = -1
            d = xp.reshape(d, tuple(shape))
        elif len(x.shape) != len(y.shape):
            raise ValueError("If given, shape of x must be 1-D or the "
                             "same as y.")
        else:
            d = xp.diff(x, axis=axis)

        if d.shape[axis] != y.shape[axis] - 1:
            raise ValueError("If given, length of x along axis must be the "
                             "same as y.")

    nd = len(y.shape)
    slice1 = tupleset((slice(None),)*nd, axis, slice(1, None))
    slice2 = tupleset((slice(None),)*nd, axis, slice(None, -1))
    res = xp.cumulative_sum(d * (y[slice1] + y[slice2]) / 2.0, axis=axis)

    if initial is not None:
        if initial != 0:
            raise ValueError("`initial` must be `None` or `0`.")
        if not np.isscalar(initial):
            raise ValueError("`initial` parameter should be a scalar.")

        shape = list(res.shape)
        shape[axis] = 1
        res = xp.concat((xp.full(tuple(shape), initial, dtype=res.dtype), res),
                        axis=axis)

    return res


def _basic_simpson(y, start, stop, x, dx, axis):
    nd = len(y.shape)
    if start is None:
        start = 0
    step = 2
    slice_all = (slice(None),)*nd
    slice0 = tupleset(slice_all, axis, slice(start, stop, step))
    slice1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
    slice2 = tupleset(slice_all, axis, slice(start+2, stop+2, step))

    if x is None:  # Even-spaced Simpson's rule.
        result = np.sum(y[slice0] + 4.0*y[slice1] + y[slice2], axis=axis)
        result *= dx / 3.0
    else:
        # Account for possibly different spacings.
        #    Simpson's rule changes a bit.
        h = np.diff(x, axis=axis)
        sl0 = tupleset(slice_all, axis, slice(start, stop, step))
        sl1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
        h0 = h[sl0].astype(float, copy=False)
        h1 = h[sl1].astype(float, copy=False)
        hsum = h0 + h1
        hprod = h0 * h1
        h0divh1 = np.true_divide(h0, h1, out=np.zeros_like(h0), where=h1 != 0)
        tmp = hsum/6.0 * (y[slice0] *
                          (2.0 - np.true_divide(1.0, h0divh1,
                                                out=np.zeros_like(h0divh1),
                                                where=h0divh1 != 0)) +
                          y[slice1] * (hsum *
                                       np.true_divide(hsum, hprod,
                                                      out=np.zeros_like(hsum),
                                                      where=hprod != 0)) +
                          y[slice2] * (2.0 - h0divh1))
        result = np.sum(tmp, axis=axis)
    return result


@xp_capabilities(np_only=True)
def simpson(y, x=None, *, dx=1.0, axis=-1):
    """
    Integrate y(x) using samples along the given axis and the composite
    Simpson's rule. If x is None, spacing of dx is assumed.

    Parameters
    ----------
    y : array_like
        Array to be integrated.
    x : array_like, optional
        If given, the points at which `y` is sampled.
    dx : float, optional
        Spacing of integration points along axis of `x`. Only used when
        `x` is None. Default is 1.
    axis : int, optional
        Axis along which to integrate. Default is the last axis.

    Returns
    -------
    float
        The estimated integral computed with the composite Simpson's rule.

    See Also
    --------
    quad : adaptive quadrature using QUADPACK
    fixed_quad : fixed-order Gaussian quadrature
    dblquad : double integrals
    tplquad : triple integrals
    romb : integrators for sampled data
    cumulative_trapezoid : cumulative integration for sampled data
    cumulative_simpson : cumulative integration using Simpson's 1/3 rule

    Notes
    -----
    For an odd number of samples that are equally spaced the result is
    exact if the function is a polynomial of order 3 or less. If
    the samples are not equally spaced, then the result is exact only
    if the function is a polynomial of order 2 or less.

    References
    ----------
    .. [1] Cartwright, Kenneth V. Simpson's Rule Cumulative Integration with
           MS Excel and Irregularly-spaced Data. Journal of Mathematical
           Sciences and Mathematics Education. 12 (2): 1-9

    Examples
    --------
    >>> from scipy import integrate
    >>> import numpy as np
    >>> x = np.arange(0, 10)
    >>> y = np.arange(0, 10)

    >>> integrate.simpson(y, x=x)
    40.5

    >>> y = np.power(x, 3)
    >>> integrate.simpson(y, x=x)
    1640.5
    >>> integrate.quad(lambda x: x**3, 0, 9)[0]
    1640.25

    """
    y = np.asarray(y)
    nd = len(y.shape)
    N = y.shape[axis]
    last_dx = dx
    returnshape = 0
    if x is not None:
        x = np.asarray(x)
        if len(x.shape) == 1:
            shapex = [1] * nd
            shapex[axis] = x.shape[0]
            saveshape = x.shape
            returnshape = 1
            x = x.reshape(tuple(shapex))
        elif len(x.shape) != len(y.shape):
            raise ValueError("If given, shape of x must be 1-D or the "
                             "same as y.")
        if x.shape[axis] != N:
            raise ValueError("If given, length of x along axis must be the "
                             "same as y.")

    if N % 2 == 0:
        val = 0.0
        result = 0.0
        slice_all = (slice(None),) * nd

        if N == 2:
            # need at least 3 points in integration axis to form parabolic
            # segment. If there are two points then any of 'avg', 'first',
            # 'last' should give the same result.
            slice1 = tupleset(slice_all, axis, -1)
            slice2 = tupleset(slice_all, axis, -2)
            if x is not None:
                last_dx = x[slice1] - x[slice2]
            val += 0.5 * last_dx * (y[slice1] + y[slice2])
        else:
            # use Simpson's rule on first intervals
            result = _basic_simpson(y, 0, N-3, x, dx, axis)

            slice1 = tupleset(slice_all, axis, -1)
            slice2 = tupleset(slice_all, axis, -2)
            slice3 = tupleset(slice_all, axis, -3)

            h = np.asarray([dx, dx], dtype=np.float64)
            if x is not None:
                # grab the last two spacings from the appropriate axis
                hm2 = tupleset(slice_all, axis, slice(-2, -1, 1))
                hm1 = tupleset(slice_all, axis, slice(-1, None, 1))

                diffs = np.float64(np.diff(x, axis=axis))
                h = [np.squeeze(diffs[hm2], axis=axis),
                     np.squeeze(diffs[hm1], axis=axis)]

            # This is the correction for the last interval according to
            # Cartwright.
            # However, I used the equations given at
            # https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule_for_irregularly_spaced_data
            # A footnote on Wikipedia says:
            # Cartwright 2017, Equation 8. The equation in Cartwright is
            # calculating the first interval whereas the equations in the
            # Wikipedia article are adjusting for the last integral. If the
            # proper algebraic substitutions are made, the equation results in
            # the values shown.
            num = 2 * h[1] ** 2 + 3 * h[0] * h[1]
            den = 6 * (h[1] + h[0])
            alpha = np.true_divide(
                num,
                den,
                out=np.zeros_like(den),
                where=den != 0
            )

            num = h[1] ** 2 + 3.0 * h[0] * h[1]
            den = 6 * h[0]
            beta = np.true_divide(
                num,
                den,
                out=np.zeros_like(den),
                where=den != 0
            )

            num = 1 * h[1] ** 3
            den = 6 * h[0] * (h[0] + h[1])
            eta = np.true_divide(
                num,
                den,
                out=np.zeros_like(den),
                where=den != 0
            )

            result += alpha*y[slice1] + beta*y[slice2] - eta*y[slice3]

        result += val
    else:
        result = _basic_simpson(y, 0, N-2, x, dx, axis)
    if returnshape:
        x = x.reshape(saveshape)
    return result


def _cumulatively_sum_simpson_integrals(
    y: np.ndarray,
    dx: np.ndarray,
    integration_func: Callable[[np.ndarray, np.ndarray], np.ndarray],
    xp
) -> np.ndarray:
    """Calculate cumulative sum of Simpson integrals.
    Takes as input the integration function to be used.
    The integration_func is assumed to return the cumulative sum using
    composite Simpson's rule. Assumes the axis of summation is -1.
    """
    sub_integrals_h1 = integration_func(y, dx)
    sub_integrals_h2 = xp.flip(
        integration_func(xp.flip(y, axis=-1), xp.flip(dx, axis=-1)),
        axis=-1
    )

    shape = list(sub_integrals_h1.shape)
    shape[-1] += 1
    sub_integrals = xp.empty(shape, dtype=xp.result_type(y, dx))
    sub_integrals[..., :-1:2] = sub_integrals_h1[..., ::2]
    sub_integrals[..., 1::2] = sub_integrals_h2[..., ::2]
    # Integral over last subinterval can only be calculated from
    # formula for h2
    sub_integrals[..., -1] = sub_integrals_h2[..., -1]
    res = xp.cumulative_sum(sub_integrals, axis=-1)
    return res


def _cumulative_simpson_equal_intervals(y: np.ndarray, dx: np.ndarray) -> np.ndarray:
    """Calculate the Simpson integrals for all h1 intervals assuming equal interval
    widths. The function can also be used to calculate the integral for all
    h2 intervals by reversing the inputs, `y` and `dx`.
    """
    d = dx[..., :-1]
    f1 = y[..., :-2]
    f2 = y[..., 1:-1]
    f3 = y[..., 2:]

    # Calculate integral over the subintervals (eqn (10) of Reference [2])
    return d / 3 * (5 * f1 / 4 + 2 * f2 - f3 / 4)


def _cumulative_simpson_unequal_intervals(y: np.ndarray, dx: np.ndarray) -> np.ndarray:
    """Calculate the Simpson integrals for all h1 intervals assuming unequal interval
    widths. The function can also be used to calculate the integral for all
    h2 intervals by reversing the inputs, `y` and `dx`.
    """
    x21 = dx[..., :-1]
    x32 = dx[..., 1:]
    f1 = y[..., :-2]
    f2 = y[..., 1:-1]
    f3 = y[..., 2:]

    x31 = x21 + x32
    x21_x31 = x21/x31
    x21_x32 = x21/x32
    x21x21_x31x32 = x21_x31 * x21_x32

    # Calculate integral over the subintervals (eqn (8) of Reference [2])
    coeff1 = 3 - x21_x31
    coeff2 = 3 + x21x21_x31x32 + x21_x31
    coeff3 = -x21x21_x31x32

    return x21/6 * (coeff1*f1 + coeff2*f2 + coeff3*f3)


@xp_capabilities(allow_dask_compute=1,
                 skip_backends=[("jax.numpy", "item assignment")])
def cumulative_simpson(y, *, x=None, dx=1.0, axis=-1, initial=None):
    r"""
    Cumulatively integrate y(x) using the composite Simpson's 1/3 rule.
    The integral of the samples at every point is calculated by assuming a
    quadratic relationship between each point and the two adjacent points.

    Parameters
    ----------
    y : array_like
        Values to integrate. Requires at least one point along `axis`. If two or fewer
        points are provided along `axis`, Simpson's integration is not possible and the
        result is calculated with `cumulative_trapezoid`.
    x : array_like, optional
        The coordinate to integrate along. Must have the same shape as `y` or
        must be 1D with the same length as `y` along `axis`. `x` must also be
        strictly increasing along `axis`.
        If `x` is None (default), integration is performed using spacing `dx`
        between consecutive elements in `y`.
    dx : scalar or array_like, optional
        Spacing between elements of `y`. Only used if `x` is None. Can either
        be a float, or an array with the same shape as `y`, but of length one along
        `axis`. Default is 1.0.
    axis : int, optional
        Specifies the axis to integrate along. Default is -1 (last axis).
    initial : scalar or array_like, optional
        If given, insert this value at the beginning of the returned result,
        and add it to the rest of the result. Default is None, which means no
        value at ``x[0]`` is returned and `res` has one element less than `y`
        along the axis of integration. Can either be a float, or an array with
        the same shape as `y`, but of length one along `axis`.

    Returns
    -------
    res : ndarray
        The result of cumulative integration of `y` along `axis`.
        If `initial` is None, the shape is such that the axis of integration
        has one less value than `y`. If `initial` is given, the shape is equal
        to that of `y`.

    See Also
    --------
    numpy.cumsum
    cumulative_trapezoid : cumulative integration using the composite
        trapezoidal rule
    simpson : integrator for sampled data using the Composite Simpson's Rule

    Notes
    -----

    .. versionadded:: 1.12.0

    The composite Simpson's 1/3 method can be used to approximate the definite
    integral of a sampled input function :math:`y(x)` [1]_. The method assumes
    a quadratic relationship over the interval containing any three consecutive
    sampled points.

    Consider three consecutive points:
    :math:`(x_1, y_1), (x_2, y_2), (x_3, y_3)`.

    Assuming a quadratic relationship over the three points, the integral over
    the subinterval between :math:`x_1` and :math:`x_2` is given by formula
    (8) of [2]_:

    .. math::
        \int_{x_1}^{x_2} y(x) dx\ &= \frac{x_2-x_1}{6}\left[\
        \left\{3-\frac{x_2-x_1}{x_3-x_1}\right\} y_1 + \
        \left\{3 + \frac{(x_2-x_1)^2}{(x_3-x_2)(x_3-x_1)} + \
        \frac{x_2-x_1}{x_3-x_1}\right\} y_2\\
        - \frac{(x_2-x_1)^2}{(x_3-x_2)(x_3-x_1)} y_3\right]

    The integral between :math:`x_2` and :math:`x_3` is given by swapping
    appearances of :math:`x_1` and :math:`x_3`. The integral is estimated
    separately for each subinterval and then cumulatively summed to obtain
    the final result.

    For samples that are equally spaced, the result is exact if the function
    is a polynomial of order three or less [1]_ and the number of subintervals
    is even. Otherwise, the integral is exact for polynomials of order two or
    less.

    References
    ----------
    .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Simpson's_rule
    .. [2] Cartwright, Kenneth V. Simpson's Rule Cumulative Integration with
            MS Excel and Irregularly-spaced Data. Journal of Mathematical
            Sciences and Mathematics Education. 12 (2): 1-9

    Examples
    --------
    >>> from scipy import integrate
    >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> x = np.linspace(-2, 2, num=20)
    >>> y = x**2
    >>> y_int = integrate.cumulative_simpson(y, x=x, initial=0)
    >>> fig, ax = plt.subplots()
    >>> ax.plot(x, y_int, 'ro', x, x**3/3 - (x[0])**3/3, 'b-')
    >>> ax.grid()
    >>> plt.show()

    The output of `cumulative_simpson` is similar to that of iteratively
    calling `simpson` with successively higher upper limits of integration, but
    not identical.

    >>> def cumulative_simpson_reference(y, x):
    ...     return np.asarray([integrate.simpson(y[:i], x=x[:i])
    ...                        for i in range(2, len(y) + 1)])
    >>>
    >>> rng = np.random.default_rng(354673834679465)
    >>> x, y = rng.random(size=(2, 10))
    >>> x.sort()
    >>>
    >>> res = integrate.cumulative_simpson(y, x=x)
    >>> ref = cumulative_simpson_reference(y, x)
    >>> equal = np.abs(res - ref) < 1e-15
    >>> equal  # not equal when `simpson` has even number of subintervals
    array([False,  True, False,  True, False,  True, False,  True,  True])

    This is expected: because `cumulative_simpson` has access to more
    information than `simpson`, it can typically produce more accurate
    estimates of the underlying integral over subintervals.

    """
    xp = array_namespace(y)
    y = xp_promote(y, force_floating=True, xp=xp)

    # validate `axis` and standardize to work along the last axis
    original_y = y
    original_shape = y.shape
    try:
        y = xp_swapaxes(y, axis, -1, xp)
    except IndexError as e:
        message = f"`axis={axis}` is not valid for `y` with `y.ndim={y.ndim}`."
        raise ValueError(message) from e
    if y.shape[-1] < 3:
        res = cumulative_trapezoid(original_y, x, dx=dx, axis=axis, initial=None)
        res = xp_swapaxes(res, axis, -1, xp)

    elif x is not None:
        x = xp_promote(x, force_floating=True, xp=xp)
        message = ("If given, shape of `x` must be the same as `y` or 1-D with "
                   "the same length as `y` along `axis`.")
        if not (x.shape == original_shape
                or (x.ndim == 1 and x.shape[0] == original_shape[axis])):
            raise ValueError(message)

        x = xp.broadcast_to(x, y.shape) if x.ndim == 1 else xp_swapaxes(x, axis, -1, xp)
        dx = xp.diff(x, axis=-1)
        if xp.any(dx <= 0):
            raise ValueError("Input x must be strictly increasing.")
        res = _cumulatively_sum_simpson_integrals(
            y, dx, _cumulative_simpson_unequal_intervals, xp
        )

    else:
        dx = xp_promote(xp.asarray(dx), force_floating=True, xp=xp)
        final_dx_shape = tupleset(original_shape, axis, original_shape[axis] - 1)
        alt_input_dx_shape = tupleset(original_shape, axis, 1)
        message = ("If provided, `dx` must either be a scalar or have the same "
                   "shape as `y` but with only 1 point along `axis`.")
        if not (dx.ndim == 0 or dx.shape == alt_input_dx_shape):
            raise ValueError(message)
        dx = xp.broadcast_to(dx, final_dx_shape)
        dx = xp_swapaxes(dx, axis, -1, xp)
        res = _cumulatively_sum_simpson_integrals(
            y, dx, _cumulative_simpson_equal_intervals, xp
        )

    if initial is not None:
        initial = xp_promote(initial, force_floating=True, xp=xp)
        alt_initial_input_shape = tupleset(original_shape, axis, 1)
        message = ("If provided, `initial` must either be a scalar or have the "
                   "same shape as `y` but with only 1 point along `axis`.")
        if not (initial.ndim == 0 or initial.shape == alt_initial_input_shape):
            raise ValueError(message)
        initial = xp.broadcast_to(initial, alt_initial_input_shape)
        initial = xp_swapaxes(initial, axis, -1, xp)

        res += initial
        res = xp.concat((initial, res), axis=-1)

    res = xp_swapaxes(res, -1, axis, xp)
    return res


@xp_capabilities()
def romb(y, dx=1.0, axis=-1, show=False):
    """
    Romberg integration using samples of a function.

    Parameters
    ----------
    y : array_like
        A vector of ``2**k + 1`` equally-spaced samples of a function.
    dx : float, optional
        The sample spacing. Default is 1.
    axis : int, optional
        The axis along which to integrate. Default is -1 (last axis).
    show : bool, optional
        When `y` is a single 1-D array, then if this argument is True
        print the table showing Richardson extrapolation from the
        samples. Default is False.

    Returns
    -------
    romb : ndarray
        The integrated result for `axis`.

    See Also
    --------
    quad : adaptive quadrature using QUADPACK
    fixed_quad : fixed-order Gaussian quadrature
    dblquad : double integrals
    tplquad : triple integrals
    simpson : integrators for sampled data
    cumulative_trapezoid : cumulative integration for sampled data

    Examples
    --------
    >>> from scipy import integrate
    >>> import numpy as np
    >>> x = np.arange(10, 14.25, 0.25)
    >>> y = np.arange(3, 12)

    >>> integrate.romb(y)
    56.0

    >>> y = np.sin(np.power(x, 2.5))
    >>> integrate.romb(y)
    -0.742561336672229

    >>> integrate.romb(y, show=True)
    Richardson Extrapolation Table for Romberg Integration
    ======================================================
    -0.81576
     4.63862  6.45674
    -1.10581 -3.02062 -3.65245
    -2.57379 -3.06311 -3.06595 -3.05664
    -1.34093 -0.92997 -0.78776 -0.75160 -0.74256
    ======================================================
    -0.742561336672229  # may vary

    >>> integrate.romb([[1, 2, 3], [4, 5, 6]], show=True)
    *** Printing table only supported for integrals of a single data set.
    array([ 4., 10.])
    """
    xp = array_namespace(y)
    y = xp.asarray(y)
    nd = len(y.shape)
    Nsamps = y.shape[axis]
    Ninterv = Nsamps-1
    n = 1
    k = 0
    while n < Ninterv:
        n <<= 1
        k += 1
    if n != Ninterv:
        raise ValueError("Number of samples must be one plus a "
                         "non-negative power of 2.")

    R = {}
    slice_all = (slice(None),) * nd
    slice0 = tupleset(slice_all, axis, 0)
    slicem1 = tupleset(slice_all, axis, -1)
    h = Ninterv * xp.asarray(dx, dtype=xp.float64)
    R[(0, 0)] = (y[slice0] + y[slicem1])/2.0*h
    slice_R = slice_all
    start = stop = step = Ninterv
    for i in range(1, k+1):
        start >>= 1
        slice_R = tupleset(slice_R, axis, slice(start, stop, step))
        step >>= 1
        R[(i, 0)] = 0.5*(R[(i-1, 0)] + h*xp.sum(y[slice_R], axis=axis))
        for j in range(1, i+1):
            prev = R[(i, j-1)]
            R[(i, j)] = prev + (prev-R[(i-1, j-1)]) / ((1 << (2*j))-1)
        h /= 2.0

    if show:
        if not np.isscalar(R[(0, 0)]):
            print("*** Printing table only supported for integrals" +
                  " of a single data set.")
        else:
            try:
                precis = show[0]
            except (TypeError, IndexError):
                precis = 5
            try:
                width = show[1]
            except (TypeError, IndexError):
                width = 8
            formstr = f"%{width}.{precis}f"

            title = "Richardson Extrapolation Table for Romberg Integration"
            print(title, "=" * len(title), sep="\n", end="\n")
            for i in range(k+1):
                for j in range(i+1):
                    print(formstr % R[(i, j)], end=" ")
                print()
            print("=" * len(title))

    return R[(k, k)]


# Coefficients for Newton-Cotes quadrature
#
# These are the points being used
#  to construct the local interpolating polynomial
#  a are the weights for Newton-Cotes integration
#  B is the error coefficient.
#  error in these coefficients grows as N gets larger.
#  or as samples are closer and closer together

# You can use maxima to find these rational coefficients
#  for equally spaced data using the commands
#  a(i,N) := (integrate(product(r-j,j,0,i-1) * product(r-j,j,i+1,N),r,0,N)
#             / ((N-i)! * i!) * (-1)^(N-i));
#  Be(N) := N^(N+2)/(N+2)! * (N/(N+3) - sum((i/N)^(N+2)*a(i,N),i,0,N));
#  Bo(N) := N^(N+1)/(N+1)! * (N/(N+2) - sum((i/N)^(N+1)*a(i,N),i,0,N));
#  B(N) := (if (mod(N,2)=0) then Be(N) else Bo(N));
#
# pre-computed for equally-spaced weights
#
# num_a, den_a, int_a, num_B, den_B = _builtincoeffs[N]
#
#  a = num_a*array(int_a)/den_a
#  B = num_B*1.0 / den_B
#
#  integrate(f(x),x,x_0,x_N) = dx*sum(a*f(x_i)) + B*(dx)^(2k+3) f^(2k+2)(x*)
#    where k = N // 2
#
_builtincoeffs = {
    1: (1,2,[1,1],-1,12),
    2: (1,3,[1,4,1],-1,90),
    3: (3,8,[1,3,3,1],-3,80),
    4: (2,45,[7,32,12,32,7],-8,945),
    5: (5,288,[19,75,50,50,75,19],-275,12096),
    6: (1,140,[41,216,27,272,27,216,41],-9,1400),
    7: (7,17280,[751,3577,1323,2989,2989,1323,3577,751],-8183,518400),
    8: (4,14175,[989,5888,-928,10496,-4540,10496,-928,5888,989],
        -2368,467775),
    9: (9,89600,[2857,15741,1080,19344,5778,5778,19344,1080,
                 15741,2857], -4671, 394240),
    10: (5,299376,[16067,106300,-48525,272400,-260550,427368,
                   -260550,272400,-48525,106300,16067],
         -673175, 163459296),
    11: (11,87091200,[2171465,13486539,-3237113, 25226685,-9595542,
                      15493566,15493566,-9595542,25226685,-3237113,
                      13486539,2171465], -2224234463, 237758976000),
    12: (1, 5255250, [1364651,9903168,-7587864,35725120,-51491295,
                      87516288,-87797136,87516288,-51491295,35725120,
                      -7587864,9903168,1364651], -3012, 875875),
    13: (13, 402361344000,[8181904909, 56280729661, -31268252574,
                           156074417954,-151659573325,206683437987,
                           -43111992612,-43111992612,206683437987,
                           -151659573325,156074417954,-31268252574,
                           56280729661,8181904909], -2639651053,
         344881152000),
    14: (7, 2501928000, [90241897,710986864,-770720657,3501442784,
                         -6625093363,12630121616,-16802270373,19534438464,
                         -16802270373,12630121616,-6625093363,3501442784,
                         -770720657,710986864,90241897], -3740727473,
         1275983280000)
    }


@xp_capabilities(np_only=True)
def newton_cotes(rn, equal=0):
    r"""
    Return weights and error coefficient for Newton-Cotes integration.

    Suppose we have :math:`(N+1)` samples of :math:`f` at the positions
    :math:`x_0, x_1, ..., x_N`. Then an :math:`N`-point Newton-Cotes formula
    for the integral between :math:`x_0` and :math:`x_N` is:

    :math:`\int_{x_0}^{x_N} f(x)dx = \Delta x \sum_{i=0}^{N} a_i f(x_i)
    + B_N (\Delta x)^{N+2} f^{N+1} (\xi)`

    where :math:`\xi \in [x_0,x_N]`
    and :math:`\Delta x = \frac{x_N-x_0}{N}` is the average samples spacing.

    If the samples are equally-spaced and N is even, then the error
    term is :math:`B_N (\Delta x)^{N+3} f^{N+2}(\xi)`.

    Parameters
    ----------
    rn : int
        The integer order for equally-spaced data or the relative positions of
        the samples with the first sample at 0 and the last at N, where N+1 is
        the length of `rn`. N is the order of the Newton-Cotes integration.
    equal : int, optional
        Set to 1 to enforce equally spaced data.

    Returns
    -------
    an : ndarray
        1-D array of weights to apply to the function at the provided sample
        positions.
    B : float
        Error coefficient.

    Notes
    -----
    Normally, the Newton-Cotes rules are used on smaller integration
    regions and a composite rule is used to return the total integral.

    Examples
    --------
    Compute the integral of sin(x) in [0, :math:`\pi`]:

    >>> from scipy.integrate import newton_cotes
    >>> import numpy as np
    >>> def f(x):
    ...     return np.sin(x)
    >>> a = 0
    >>> b = np.pi
    >>> exact = 2
    >>> for N in [2, 4, 6, 8, 10]:
    ...     x = np.linspace(a, b, N + 1)
    ...     an, B = newton_cotes(N, 1)
    ...     dx = (b - a) / N
    ...     quad = dx * np.sum(an * f(x))
    ...     error = abs(quad - exact)
    ...     print('{:2d}  {:10.9f}  {:.5e}'.format(N, quad, error))
    ...
     2   2.094395102   9.43951e-02
     4   1.998570732   1.42927e-03
     6   2.000017814   1.78136e-05
     8   1.999999835   1.64725e-07
    10   2.000000001   1.14677e-09

    """
    try:
        N = len(rn)-1
        if equal:
            rn = np.arange(N+1)
        elif np.all(np.diff(rn) == 1):
            equal = 1
    except Exception:
        N = rn
        rn = np.arange(N+1)
        equal = 1

    if equal and N in _builtincoeffs:
        na, da, vi, nb, db = _builtincoeffs[N]
        an = na * np.array(vi, dtype=float) / da
        return an, float(nb)/db

    if (rn[0] != 0) or (rn[-1] != N):
        raise ValueError("The sample positions must start at 0"
                         " and end at N")
    yi = rn / float(N)
    ti = 2 * yi - 1
    nvec = np.arange(N+1)
    C = ti ** nvec[:, np.newaxis]
    Cinv = np.linalg.inv(C)
    # improve precision of result
    for i in range(2):
        Cinv = 2*Cinv - Cinv.dot(C).dot(Cinv)
    vec = 2.0 / (nvec[::2]+1)
    ai = Cinv[:, ::2].dot(vec) * (N / 2.)

    if (N % 2 == 0) and equal:
        BN = N/(N+3.)
        power = N+2
    else:
        BN = N/(N+2.)
        power = N+1

    BN = BN - np.dot(yi**power, ai)
    p1 = power+1
    fac = power*math.log(N) - gammaln(p1)
    fac = math.exp(fac)
    return ai, BN*fac


def _qmc_quad_iv(func, a, b, n_points, n_estimates, qrng, log, xp):
    # lazy import to avoid issues with partially-initialized submodule
    if not hasattr(qmc_quad, 'qmc'):
        from scipy import stats
        qmc_quad.stats = stats
    else:
        stats = qmc_quad.stats

    if not callable(func):
        message = "`func` must be callable."
        raise TypeError(message)

    # a, b will be modified, so copy. Oh well if it's copied twice.
    a, b = xp_promote(a, b, broadcast=True, force_floating=True, xp=xp)
    a = xpx.atleast_nd(a, ndim=1, xp=xp)
    b = xpx.atleast_nd(b, ndim=1, xp=xp)
    a, b = xp.broadcast_arrays(a, b)
    a, b = xp_copy(a), xp_copy(b)
    dim = a.shape[0]

    try:
        func((a + b) / 2)
    except Exception as e:
        message = ("`func` must evaluate the integrand at points within "
                   "the integration range; e.g. `func( (a + b) / 2)` "
                   "must return the integrand at the centroid of the "
                   "integration volume.")
        raise ValueError(message) from e

    try:
        func(xp.stack([a, b]).T)
        vfunc = func
    except Exception as e:
        if not is_numpy(xp):
            message = ("Exception encountered when attempting vectorized call to "
                       f"`func`: {e}. When using array library {xp}, `func` must "
                       "accept two-dimensional array `x` with shape `(a.shape[0], "
                       "n_points)` and return an array of the integrand value at "
                       "each of the `n_points`.")
            raise ValueError(message)

        message = ("Exception encountered when attempting vectorized call to "
                   f"`func`: {e}. For better performance, `func` should "
                   "accept two-dimensional array `x` with shape `(len(a), "
                   "n_points)` and return an array of the integrand value at "
                   "each of the `n_points`.")
        warnings.warn(message, stacklevel=3)

        def vfunc(x):
            return np.apply_along_axis(func, axis=-1, arr=x)

    n_points_int = int(n_points)
    if n_points != n_points_int:
        message = "`n_points` must be an integer."
        raise TypeError(message)

    n_estimates_int = int(n_estimates)
    if n_estimates != n_estimates_int:
        message = "`n_estimates` must be an integer."
        raise TypeError(message)

    if qrng is None:
        qrng = stats.qmc.Halton(dim)
    elif not isinstance(qrng, stats.qmc.QMCEngine):
        message = "`qrng` must be an instance of scipy.stats.qmc.QMCEngine."
        raise TypeError(message)

    if qrng.d != a.shape[0]:
        message = ("`qrng` must be initialized with dimensionality equal to "
                   "the number of variables in `a`, i.e., "
                   "`qrng.random().shape[-1]` must equal `a.shape[0]`.")
        raise ValueError(message)

    rng_seed = getattr(qrng, 'rng_seed', None)
    rng = stats._qmc.check_random_state(rng_seed)

    if log not in {True, False}:
        message = "`log` must be boolean (`True` or `False`)."
        raise TypeError(message)

    return (vfunc, a, b, n_points_int, n_estimates_int, qrng, rng, log, stats)


QMCQuadResult = namedtuple('QMCQuadResult', ['integral', 'standard_error'])


@xp_capabilities(skip_backends=[("dask.array",
                                 "Dask arrays are confused about their shape")],
                 jax_jit=False)
def qmc_quad(func, a, b, *, n_estimates=8, n_points=1024, qrng=None,
             log=False):
    """
    Compute an integral in N-dimensions using Quasi-Monte Carlo quadrature.

    Parameters
    ----------
    func : callable
        The integrand. Must accept a single argument ``x``, an array which
        specifies the point(s) at which to evaluate the scalar-valued
        integrand, and return the value(s) of the integrand.
        For efficiency, the function should be vectorized to accept an array of
        shape ``(d, n_points)``, where ``d`` is the number of variables (i.e.
        the dimensionality of the function domain) and `n_points` is the number
        of quadrature points, and return an array of shape ``(n_points,)``,
        the integrand at each quadrature point.
    a, b : array-like
        One-dimensional arrays specifying the lower and upper integration
        limits, respectively, of each of the ``d`` variables.
    n_estimates, n_points : int, optional
        `n_estimates` (default: 8) statistically independent QMC samples, each
        of `n_points` (default: 1024) points, will be generated by `qrng`.
        The total number of points at which the integrand `func` will be
        evaluated is ``n_points * n_estimates``. See Notes for details.
    qrng : `~scipy.stats.qmc.QMCEngine`, optional
        An instance of the QMCEngine from which to sample QMC points.
        The QMCEngine must be initialized to a number of dimensions ``d``
        corresponding with the number of variables ``x1, ..., xd`` passed to
        `func`.
        The provided QMCEngine is used to produce the first integral estimate.
        If `n_estimates` is greater than one, additional QMCEngines are
        spawned from the first (with scrambling enabled, if it is an option.)
        If a QMCEngine is not provided, the default `scipy.stats.qmc.Halton`
        will be initialized with the number of dimensions determine from
        the length of `a`.
    log : boolean, default: False
        When set to True, `func` returns the log of the integrand, and
        the result object contains the log of the integral.

    Returns
    -------
    result : object
        A result object with attributes:

        integral : float
            The estimate of the integral.
        standard_error :
            The error estimate. See Notes for interpretation.

    Notes
    -----
    Values of the integrand at each of the `n_points` points of a QMC sample
    are used to produce an estimate of the integral. This estimate is drawn
    from a population of possible estimates of the integral, the value of
    which we obtain depends on the particular points at which the integral
    was evaluated. We perform this process `n_estimates` times, each time
    evaluating the integrand at different scrambled QMC points, effectively
    drawing i.i.d. random samples from the population of integral estimates.
    The sample mean :math:`m` of these integral estimates is an
    unbiased estimator of the true value of the integral, and the standard
    error of the mean :math:`s` of these estimates may be used to generate
    confidence intervals using the t distribution with ``n_estimates - 1``
    degrees of freedom. Perhaps counter-intuitively, increasing `n_points`
    while keeping the total number of function evaluation points
    ``n_points * n_estimates`` fixed tends to reduce the actual error, whereas
    increasing `n_estimates` tends to decrease the error estimate.

    Examples
    --------
    QMC quadrature is particularly useful for computing integrals in higher
    dimensions. An example integrand is the probability density function
    of a multivariate normal distribution.

    >>> import numpy as np
    >>> from scipy import stats
    >>> dim = 8
    >>> mean = np.zeros(dim)
    >>> cov = np.eye(dim)
    >>> def func(x):
    ...     # `multivariate_normal` expects the _last_ axis to correspond with
    ...     # the dimensionality of the space, so `x` must be transposed
    ...     return stats.multivariate_normal.pdf(x.T, mean, cov)

    To compute the integral over the unit hypercube:

    >>> from scipy.integrate import qmc_quad
    >>> a = np.zeros(dim)
    >>> b = np.ones(dim)
    >>> rng = np.random.default_rng()
    >>> qrng = stats.qmc.Halton(d=dim, seed=rng)
    >>> n_estimates = 8
    >>> res = qmc_quad(func, a, b, n_estimates=n_estimates, qrng=qrng)
    >>> res.integral, res.standard_error
    (0.00018429555666024108, 1.0389431116001344e-07)

    A two-sided, 99% confidence interval for the integral may be estimated
    as:

    >>> t = stats.t(df=n_estimates-1, loc=res.integral,
    ...             scale=res.standard_error)
    >>> t.interval(0.99)
    (0.0001839319802536469, 0.00018465913306683527)

    Indeed, the value reported by `scipy.stats.multivariate_normal` is
    within this range.

    >>> stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
    0.00018430867675187443

    """
    xp = array_namespace(a, b)
    args = _qmc_quad_iv(func, a, b, n_points, n_estimates, qrng, log, xp)
    func, a, b, n_points, n_estimates, qrng, rng, log, stats = args

    def sum_product(integrands, dA, log=False):
        if log:
            return logsumexp(integrands) + math.log(dA)
        else:
            return xp.sum(integrands * dA)

    def mean(estimates, log=False):
        if log:
            return logsumexp(estimates) - math.log(n_estimates)
        else:
            return xp.mean(estimates)

    def std(estimates, m=None, ddof=0, log=False):
        m = m or mean(estimates, log)
        if log:
            estimates, m = xp.broadcast_arrays(estimates, m)
            temp = xp.stack((estimates, m + xp.pi * 1j))
            diff = logsumexp(temp, axis=0)
            return xp.real(0.5 * (logsumexp(2 * diff)
                                  - math.log(n_estimates - ddof)))
        else:
            return xp.std(estimates, correction=ddof)

    def sem(estimates, m=None, s=None, log=False):
        m = m or mean(estimates, log)
        s = s or std(estimates, m, ddof=1, log=log)
        if log:
            return s - 0.5*math.log(n_estimates)
        else:
            return s / math.sqrt(n_estimates)

    # The sign of the integral depends on the order of the limits. Fix this by
    # ensuring that lower bounds are indeed lower and setting sign of resulting
    # integral manually
    if xp.any(a == b):
        message = ("A lower limit was equal to an upper limit, so the value "
                   "of the integral is zero by definition.")
        warnings.warn(message, stacklevel=2)
        zero = xp.asarray(-xp.inf if log else 0, dtype=a.dtype)
        return QMCQuadResult(zero, xp.asarray(0., dtype=a.dtype))

    i_swap = b < a
    sign = (-1)**(xp.count_nonzero(i_swap, axis=-1))  # odd # of swaps -> negative
    sign = xp.astype(sign, a.dtype)
    # a[i_swap], b[i_swap] = b[i_swap], a[i_swap]
    a_iswap = a[i_swap]
    b_iswap = b[i_swap]
    a = xpx.at(a)[i_swap].set(b_iswap)
    b = xpx.at(b)[i_swap].set(a_iswap)

    A = xp.prod(b - a)
    dA = A / n_points

    estimates = xp.zeros(n_estimates, dtype=a.dtype)
    rngs = _rng_spawn(qrng.rng, n_estimates)
    for i in range(n_estimates):
        # Generate integral estimate
        sample = xp.asarray(qrng.random(n_points), dtype=a.dtype)
        # The rationale for transposing is that this allows users to easily
        # unpack `x` into separate variables, if desired. This is consistent
        # with the `xx` array passed into the `scipy.integrate.nquad` `func`.
        x = (sample * (b - a) + a).T  # (n_dim, n_points)
        integrands = func(x)
        estimates = xpx.at(estimates)[i].set(sum_product(integrands, dA, log))

        # Get a new, independently-scrambled QRNG for next time
        qrng = type(qrng)(seed=rngs[i], **qrng._init_quad)

    integral = mean(estimates, log)
    standard_error = sem(estimates, m=integral, log=log)
    integral = integral + xp.pi*1j if (log and sign < 0) else integral*sign
    return QMCQuadResult(integral, standard_error)