Simpson's rule area and volume integrals using buildq & lambda methods

Following the suggestions of Barton Willis, the following
code seems to work on simple examples, including variable
limits.

Since I am a complete novice using buildq with lambda methods,
there are surely more efficient ways to do this, which I would be
glad to hear about.

code plus examples:
------------------------------
/* nintmd.mac */
/* simpson's one third rule code
    based on buildq-lambda forms  */



simp(expr,var,a,b,n):=
    block ([h : abs (b-a)/n,jj,
                  sf : buildq([expr,var],lambda([var],expr)) ],
         if not evenp(n) then return("n must be an even integer"),
         float((h/3)*( sf(a) + sf(b) +
                 4*sum( sf(a + (2*jj-1)*h), jj, 1, n/2 ) +
                   2*sum( sf(a + 2*jj*h), jj, 1, (n/2) - 1 ) )))$


/*
(%i1) load(nintmd);

(%o1) "c:/work2/nintmd.mac"
(%i2) float(integrate(sin(x),x,0,1));

(%o2) 0.45969769413186
(%i3) simp(sin(x),x,0,1,16);

(%o3) 0.45969773311905
*/


area(ex,uL,vL,nL) :=
 block([%%u,u1,u2,%%v,v1,v2,nu,nv,%uu,
        exu : buildq([%%u,ex,v1,v2,%%v,nv ],lambda([u0],
                   block([%expr,%v1,%v2],
                     %expr:subst(u0,%%u,ex),
                     %v1:subst(u0,%%u,v1),
                     %v2:subst(u0,%%u,v2),
                     simp(%expr,%%v,%v1,%v2,nv))))],
        [%%u,u1,u2] : uL,
        [%%v,v1,v2] : vL,
        [nu,nv] : nL,
        simp(exu(%uu),%uu,u1,u2,nu))$


/*
(%i4) float(integrate(integrate(x*y^2,y,0,1),x,0,1));

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

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

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

(%o7) 0.066667344835069

*/

volume(ex,uL,vL,wL,nL) :=
 block([%%u,u1,u2,%%v,v1,v2,%%w,w1,w2,
              nu,nv,nw,%uu,

       exuv : buildq([%%u,%%v,ex,w1,w2,%%w,nw ],lambda([us,vs],
                     block([%%expr,%w1,%w2 ],
                         %%expr:subst([%%u=us,%%v=vs],ex),
                         %w1:subst([%%u=us,%%v=vs],w1),
                          %w2:subst([%%u=us,%%v=vs],w2),
                            simp(%%expr,%%w,%w1,%w2,nw)))),

       exu : buildq([%%u,ex,v1,v2,%%v,nv ],lambda([u0],
                    block([%exuv,%v1,%v2],
                         %exuv:subst(u0,%%u,exuv(%%u,%%v)),
                         %v1:subst(u0,%%u,v1),
                         %v2:subst(u0,%%u,v2),
                          simp(%exuv,%%v,%v1,%v2,nv))))],

        [%%u,u1,u2] : uL,
        [%%v,v1,v2] : vL,
        [%%w,w1,w2] : wL,
        [nu,nv,nw] : nL,
        simp(exu(%uu),%uu,u1,u2,nu))$


/*


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

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

(%o9) 0.041666666666667
(%i10) float(integrate(integrate(integrate(x*y^2*z^3,z,0,y),y,0,x),x,0,1));

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

(%o11) 0.003969543636397

*/

-------------------------------------
Ted Woollett