cube_int.c integrate to find 4D volume 
Numerical quadrature, x=1 to x=3, y=0.5 to y=2, z=0.25 to z=2.25 

Using trapezoidal method 
f1(x,y,z)=x+2*y+3*z+4  np=4,  hx=0.666667, hy=0.5, hz=0.666667 
64 function evaluations 
computed volume=152.444, exact volume=73.5, error=78.9444 
f1(x,y,z)=x+2*y+3*z+4  np=16,  hx=0.133333, hy=0.1, hz=0.133333 
4096 function evaluations 
computed volume=76.4836, exact volume=73.5, error=2.98356 
f1(x,y,z)=x+2*y+3*z+4  np=64,  hx=0.031746, hy=0.0238095, hz=0.031746 
262144 function evaluations 
computed volume=73.6673, exact volume=73.5, error=0.167255 

f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10  np=4,  hx=0.666667, hy=0.5, hz=0.666667 
64 function evaluations 
computed volume=1146.26, exact volume=543.875, error=602.384 
f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10  np=16,  hx=0.133333, hy=0.1, hz=0.133333 
4096 function evaluations 
computed volume=566.603, exact volume=543.875, error=22.7284 
f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10  np=64,  hx=0.031746, hy=0.0238095, hz=0.031746 
262144 function evaluations 
computed volume=545.149, exact volume=543.875, error=1.27395 
f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10  np=256,  hx=0.00784314, hy=0.00588235, hz=0.00784314 
16777216 function evaluations 
computed volume=543.953, exact volume=543.875, error=0.0775517 

f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20 
          np=4,  hx=0.666667, hy=0.5, hz=0.666667 
64 function evaluations 
computed volume=6363.36, exact volume=2887.97, error=3475.39 
f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20 
          np=16,  hx=0.133333, hy=0.1, hz=0.133333 
4096 function evaluations 
computed volume=3018.55, exact volume=2887.97, error=130.577 
f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20 
          np=64,  hx=0.031746, hy=0.0238095, hz=0.031746 
262144 function evaluations 
computed volume=2895.29, exact volume=2887.97, error=7.31629 
f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20 
          np=256,  hx=0.00784314, hy=0.00588235, hz=0.00784314 
16777216 function evaluations 
computed volume=2888.41, exact volume=2887.97, error=0.445352 

fe(x,y,z)=exp(x)*exp(y)*exp(z) 
          np=4,  hx=0.666667, hy=0.5, hz=0.666667 
64 function evaluations 
computed volume=2376.9, exact volume=817.86, error=1559.04 
fe(x,y,z)=exp(x)*exp(y)*exp(z) 
          np=16,  hx=0.133333, hy=0.1, hz=0.133333 
4096 function evaluations 
computed volume=874.635, exact volume=817.86, error=56.7758 
fe(x,y,z)=exp(x)*exp(y)*exp(z) 
          np=64,  hx=0.031746, hy=0.0238095, hz=0.031746 
262144 function evaluations 
computed volume=821.032, exact volume=817.86, error=3.17229 
fe(x,y,z)=exp(x)*exp(y)*exp(z) 
          np=256,  hx=0.00784314, hy=0.00588235, hz=0.00784314 
16777216 function evaluations 
computed volume=818.053, exact volume=817.86, error=0.193006 

Using Gauss Legendre method 
f1(x,y,z)=x+2y+3z+4 np=4 
64 function evaluations 
computed volume=73.5, exact volume=73.5, error=-1.33582e-12 

f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ 
          7*x+8*y+9*z+10   np=4 
64 function evaluations 
computed volume=543.875, exact volume=543.875, error=-9.89075e-12 

f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + 
          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + 
          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=4 
64 function evaluations 
computed volume=2887.97, exact volume=2887.97, error=-5.50244e-11 

fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=4 
64 function evaluations 
computed volume=817.859, exact volume=817.86, error=-0.00021625 
fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=8 
512 function evaluations 
computed volume=817.86, exact volume=817.86, error=4.30371e-08