simpson's rule area and volume integrals using buildq and define

Thanks to the examples of Robert Dodier and
Barton Willis, I have written Simpson's one
third rule code for area and volume integrals
as a first step in improving the software  associated
with ch. 8, Maxima by Example: Numerical Integration.
The code handles variable limits of integration
for the inner integral(s).
------------------------------------------
Maxima 5.25.1 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)

(%i1) load(simpson);
(%o1) "c:/work2/simpson.mac"

/* area with numerical limits */

(%i2) integrate(integrate(x*y^2,y,0,1),x,0,1);

(%o2) 1/6
(%i3) float(%);

(%o3) 0.16666666666667

(%i4) area(ex,uL,vL,nL) :=
 block([%%u,u1,u2,%%v,v1,v2,nu,nv,u0,hh],
   local(hh),
   [%%u,u1,u2] : uL,
   [%%v,v1,v2] : vL,
   [nu,nv] : nL,
   define(hh(u0),
         buildq([%expr:subst(u0,%%u,ex),%v1:subst(u0,%%u,v1),
                     %v2:subst(u0,%%u,v2)],
                      simp1(%expr,%%v,%v1,%v2,nv))),
   simp1(hh(%uu),%uu,u1,u2,nu))$

(%i5) area(x*y^2,[x,0,1],[y,0,1],[4,4]);

(%o5) 0.16666666666667

/* area with variable limits on the inner integral */

(%i6) integrate(integrate(x*y^2,y,0,x),x,0,1);

(%o6) 1/15
(%i7) float(%);

(%o7) 0.066666666666667

(%i8) area(x*y^2,[x,0,1],[y,0,x],[16,16]);

(%o8) 0.066667344835069

/*  volume with numerical limits */

(%i9) integrate(integrate(integrate(x*y^2*z^3,z,0,1),y,0,1),x,0,1);

(%o9) 1/24
(%i10) float(%);

(%o10) 0.041666666666667

(%i11) volume(ex,uL,vL,wL,nL) :=
 block([%%u,u1,u2,%%v,v1,v2,%%w,w1,w2,
              nu,nv,nw,u0,us,vs,hh,gg],
   local(hh,gg),
   [%%u,u1,u2] : uL,
   [%%v,v1,v2] : vL,
   [%%w,w1,w2] : wL,
   [nu,nv,nw] : nL,
   define(gg(us,vs),
        buildq([%%expr:subst([u=us,v=vs],ex),
                 %w1:subst([u=us,v=vs],w1),
                 %w2:subst([u=us,v=vs],w2) ],
         simp1(%%expr,%%w,%w1,%w2,nw))),
   define(hh(u0),
         buildq([%gg:subst(u0,%%u,gg(u,v)),
                  %v1:subst(u0,%%u,v1),
                  %v2:subst(u0,%%u,v2)],
            simp1(%gg,%%v,%v1,%v2,nv))),
   simp1(hh(%uu),%uu,u1,u2,nu))$

(%i12) volume(x*y^2*z^3,[x,0,1],[y,0,1],[z,0,1],[4,4,4]);

(%o12) 0.041666666666667

/* volume with variable limits on the two inner
     integrals  */

(%i13) integrate(integrate(integrate(x*y^2*z^3,z,0,y),y,0,x),x,0,1);

(%o13) 1/252
(%i14) float(%);

(%o14) 0.003968253968254

(%i15) volume(x*y^2*z^3,[x,0,1],[y,0,x],[z,0,y],[16,16,16]);

(%o15) 0.003969543636397
--------------------------------------------------------
The one dimensional simpson's rule simp1
from simpson.mac is:
------------------------------------------
simp1(expr,var,a,b,n):=
    block ([hh,_j%],local(_ff%),
     if not evenp(n) then return("n must be an even integer"),
     define(_ff%(var),expr),     
     hh : abs (b-a)/n,
     (hh/3)*( _ff%(a) + _ff%(b) + 
             4*sum( _ff%(a + (2*_j%-1)*hh), _j%, 1, n/2 ) +
             2*sum( _ff%(a + 2*_j%*hh), _j%, 1, (n/2) - 1 ) ),
     float(%%) )$
 --------------------------------------

Ted Woollett