To recap, the 3point Simpson rule divides the interval of integration (x_{0}, x_{1}, x_{2}) into 2 equal panels of width h, and the 3 ordinates at the panel boundaries are combined using simple arithmetic
S = (h/3) (y_{0} + 4y_{1} + y_{2}) 
to approximate the actual integral I. Next, the rule is applied to the two subintervals (x_{0}, x_{1}) and (x_{1}, x_{2}), and the sum of these results, S_{ next}, i.e., the new best value, is checked to see if it is “close enough” to S, the most recent old best value. If so, S_{ next} is returned. If not, further subdivision is called for. This is the socalled “adaptive quadrature”, as the frenzy of subdividing adapts to those sections of the interval that need it most. (Note that when subdividing the interval (x_{0}, x_{1}, x_{2}) into (x_{0}, x_{0.5}, x_{1}) and (x_{1}, x_{1.5}, x_{2}), two new ordinates are needed at x_{0.5} and x_{1.5}, but the old ordinates at x_{0}, x_{1}, and x_{2} can be reused.)
NewtonCotes formulas with more points give better results. Simpson’s rule (3point) has an error term that is proportial to the 5th power of the interval times the 4th derivative of the integrand at some point ξ within the interval. What this means is that if the integrand does not have a 4th derivative (i.e., it is a polynomial of degree ≤ 3), Simpson’s rule will give accurate quadratures at the first try! If one chooses the 5point Boole’s rule, the error term is proportional to the 7th power of the interval times the 6th derivative of the integrand. So polynomials of degree ≤ 5 would get accurate quadratures at first try. In general, a NewtonCotes formula has an error term proportional to a power p of the interval, with p being 5 for Simpson, and 7 for Boole.
The ordinate summation for Boole’s rule is
(2/45) h [7y_{0} + 32y_{1} + 12y_{2} + 32y_{3} + 7y_{4}] 
which is no tougher than Simpson’s rule to program. (See Abramowitz & Stegun, Handbook of Mathematical Functions, sec. 25.4.14, Bode’s [sic] rule.)
Furthermore, Boole’s rule starts out with 4 panels, and thus is convenient for intervalhalving and ordinate sharing, just like Simpson’s rule (which starts with 2 panels).
A crucial part of any adaptive NewtonCotes quadrature is determining when the difference between S and S_{ next}, S_{ next}  S, is small enough. It stands to reason that the rule at every intervalhalving should also halve the error spec ε, so that the subintervals equitably share the burden. However, what does that mean for the S_{ next}  S check? Ward & Cheney (Numerical Mathematics and Computing, 2/e, sec. 5.4, An Adaptive Simpson’s Scheme) show how S_{ next}  S must be ≤ 15ε for Simpson’s rule. Let’s try the same steps without specifying which NewtonCotes formula is in use.
Let’s say at a certain step, for an interval of width h, the actual integral is I; the calculated integral is S; and the error spec (for that interval) is ε. As we’ve noted, the error term in S, error_term(S) = I  S is proportial to the p’th power of h.
error_term(S) = k * h^{p} 
Now, on halving the interval, we get S_{ next} = the sum of two quadratures, each of interval h/2. So, we have

In other words,
error_term(S) = (2^{p1})* error_term(S_{ next}) 
Obviously,  error_term(S_{ next}) is quite a bit smaller than  error_term(S), which is good, but what we really want is for  error_term(S_{ next}) to shrink to at least ε. Let’s now investigate the difference between S and S_{ next}.

So, if  error_term(S_{ next}) is to be ≤ ε, then simple multiplication gives
(2^{p1}  1)* error_term(S_{ next}) must be ≤ (2^{p1}  1)ε 
In other words,
S_{ next}  S must be ≤ (2^{p1}  1)ε 
For Simpson’s rule, p = 5, so we confirm
S_{ next}  S must indeed be ≤ 15ε 
For Boole’s rule, p = 7, so
S_{ next}  S must be ≤ 63ε 
One final wrinkle is due to the fact that in the course of our subdividing, we may go past the smallest representable nonzero float in our Lisp. Thus, if we find that either h or ε has decayed to machine 0, we will return whatever result we got at that point.
The entire code is in quadboole.lisp.
Compare quadsimpson.lisp.
*quadratureerror*
(i.e. the initial ε) of 1e12, i.e., 11 correct places
after the decimal.
? (quadrature (lambda (x) (/ (+ 1 (* x x)))) 0 1) ; = π/4 0.7853981633974341 ? (/ pi 4) 0.7853981633974483 ? (defun normpdf (x) (/ (exp (* 1/2 x x)) #.(sqrt (* 2L0 pi)))) NORMPDF ? (quadrature #'normpdf 0 1) 0.3413447460685443 ;octave q eval 'printf("%.20f", normcdf(1)  1/2)' gives ;0.34134474606854292578 ? (quadrature (lambda (x) (exp ( x))) 0 nil) 1.0 ? (quadrature #'exp nil 0) 1.0 ? (defun gamma (n) (cond ((< n 1) (/ (gamma (+ n 1)) n)) ((= n 1) 1) ((< 1 n 2) (let ((n1 ( n 1))) (quadrature (lambda (x) (* (expt x n1) (exp ( x)))) 0 nil))) (t (let ((n1 ( n 1))) (* n1 (gamma n1)))))) GAMMA ? (gamma 3/2) ; = (1/2)√π 0.8862269254527402 ? (* 1/2 (sqrt pi)) 0.8862269254527579