Integrate function numerically using Simpson's Rule

integrate_simpson.py

def IntegrateSimpson(x, y, *, dtype = 'float64'):
    # Integrates function f represented by (x, y) samples using Simpson's Rule.
    # See https://en.wikipedia.org/wiki/Simpson%27s_rule
    # Returns array a[i] = "integral on samples (x[0], y[0])..(x[i], y[i])".
    # Hence a[-1] (last value) equals to integral on full range.
    # First two values are approximated by Trapezoid Rule, as Simpson's Rule needs at least 3 points.
    # Needs: python -m pip install numpy
    
    import numpy as np
    
    dtype = np.dtype(dtype)
    
    x = np.array(x).ravel().astype(dtype)
    y = np.array(y).ravel().astype(dtype)
    assert x.size == y.size, (x.size, y.size)
    
    if x.size == 0:
        return np.array([], dtype = dtype)
    elif x.size == 1:
        return np.array([0], dtype = dtype)
        
    # Trapezoid approximation of first elements
    res = np.array([0, (y[0] + y[1]) * (x[1] - x[0]) / 2], dtype = dtype)
    
    N = x.size - 1
    h = np.diff(x)
    hc, hp = h[1:N:2], h[0:N-1:2]
    hs = hc + hp
        
    res = np.concatenate((res, (
        y[1:N  :2] * (hc ** 3 + hp ** 3 + 3 * hc * hp * hs) / (6 * hc * hp) +
        y[0:N-1:2] * (2 * hp ** 3 - hc ** 3 + 3 * hc * hp ** 2) / (6 * hp * hs) +
        y[2:N+1:2] * (2 * hc ** 3 - hp ** 3 + 3 * hp * hc ** 2) / (6 * hc * hs)
    ).cumsum().repeat(2)[:N-1]))
    
    hlp, hlpp = h[2:N:2], h[1:N-1:2]
    
    res[3::2] += (
        y[3:N+1:2] * (2 * hlp ** 2 + 3 * hlpp * hlp) / (6 * (hlpp + hlp)) +
        y[2:N  :2] * (hlp ** 2 + 3 * hlp * hlpp) / (6 * hlpp) -
        y[1:N-1:2] * hlp ** 3 / (6 * hlpp * (hlpp + hlp))
    )
    
    return res

# Usage Test
# Needs: python -m pip install numpy sympy
import numpy as np, sympy as sp
np.random.seed(0)
for l in range(1, 21):
    x = (np.random.random(l) / 10).cumsum()
    y = np.sin(x)
    spx = sp.Symbol('x')
    assert abs(IntegrateSimpson(x, y)[-1] - sp.N(sp.integrate(sp.sin(spx), (spx, x[0], x[-1])))) < 0.0001, l