Repository: scipy/scipy
File: scipy/integrate/_quadrature.py (Python)

Code before:
    >>> 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

    """
    xp = array_namespace(y, x, dx)
    y, x, dx = xp_promote(y, x, dx, xp=xp)
    nd = y.ndim
    N = y.shape[axis]
    last_dx = dx
    if x is not None:
        if x.ndim == 1:
            shapex = [1] * nd
            shapex[axis] = x.shape[0]
            x = xp.reshape(x, tuple(shapex))
        elif x.ndim != y.ndim:
            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 = 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, xp=xp)

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

            h = xp.stack([dx, dx])

Review comment (suggestion):
`returnshape` and `saveshape` were unused except for useless `x = x.reshape(saveshape)` at the end.

Code after:

    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

    """
    xp = array_namespace(y, x, dx)
    y, x, dx = xp_promote(y, x, dx, xp=xp)
    nd = y.ndim
    N = y.shape[axis]
    last_dx = dx
    if x is not None:
        if x.ndim == 1:
            shapex = [1] * nd
            shapex[axis] = x.shape[0]
            x = xp.reshape(x, tuple(shapex))
        elif x.ndim != y.ndim:
            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, xp=xp)