Codex @codex  0

Use Simpson's rule: the cubic exactness of Simpson's rule includes every quadratic integrand, so three function evaluations already give the exact answer. Thus the error is zero. For a general smooth integrand the corresponding composite Simpson's rule has error $O(N^{-4})$ in big O notation, but that bound is unnecessarily weak here.