/* fourd_int.c  integrate a function over a four dimensional cube     */
/*             this is numerical quadrature                           */
/*             basically compute a 5D volume under 4D function        */
/*             compare methods for accuracy                           */

/* first trapezoid rule, varying step size, int f(x,y,z,t) dx dx dz dt*/
/* then using Gauss-Legendre coordinates and weights                  */

/* remember np points or nodes means np*np*np function evaluations    */
/* or (np-1)^4  four D elements                                       */

#include <stdio.h>
#include <math.h>
#undef  abs
#define abs(x)  ((x)<0.0?(-(x)):(x))

static void gaulegf(double x1, double x2, double x[], double w[], int n);

static double f1(double x, double y, double z, double t)
{
  return x+2.0*y+3.0*z+4.0*t+5.0; /* linear, first order */
}

static double f2(double x, double y, double z, double t)
{
  return x*x+2.0*y*y+3.0*z*z+4.0*t*t+5.0*x*y+6.0*x*z+7.0*y*z+
         8.0*x*t+9.0*y*t+10.0*z*t+
         11.0*x+12.0*y+23.0*z+14.0*t+15.0; /* quadratic, second order */
}

static double f3(double x, double y, double z, double t)
{
  return x*x*x+2.0*y*y*y+3.0*z*z*z+4.0*t*t*t+5.0*x*x*y+
         6.0*x*x*z+7.0*x*y*y+8.0*x*y*z+9.0*x*z*z+
         10.0*y*y*z+11.0*y*z*z+12.0*x*x*t+13.0*x*t*t+14.0*x*y*t+
         15.0*y*y*t+16.0*y*t*t+17.0*y*z*t+18.0*z*z*t+19.0*z*t*t+
         20.0*x*z*t+21.0*x*x+22.0*y*y+23.0*z*z+24.0*x*y+25.0*x*z+
         26.0*y*z+27.0*x*t+28.0*y*t+29.0*z*t+
         30.0*x+31.0*y+32.0*z+33.0*t+34.0; /* cubic, third order */
}

static double fe(double x, double y, double z, double t)
{
  return exp(x)*exp(y)*exp(z)*exp(t);
}

int main(int argc, char *argv[])
{
  /* simple trapezoidal integration from x=1 to x=3*/
  double volume;
  double x, y, z, t, hx, hy, hz, ht;
  double xmin = 1.0;
  double xmax = 3.0;
  double ymin = 0.5;
  double ymax = 2.0;
  double zmin = 0.25;
  double zmax = 2.25;
  double tmin = 0.5;
  double tmax = 2.5;
  double f1exact = 231.0;
  double f2exact = 2684.0;
  double f3exact = 17976.25;
  double feexact = 8615.146615676;
  double xx[200], wx[200]; /* for Gauss-Legendre */
  double yy[200], wy[200]; /* for Gauss-Legendre */
  double zz[200], wz[200]; /* for Gauss-Legendre */
  double tt[200], wt[200]; /* for Gauss-Legendre */
  int i, j, k, m, np;

  printf("fourd_int.c integrate to find 5D volume under 4D surface \n");
  printf("Numerical quadrature, x=%g to x=%g, y=%g to y=%g, \n z=%g to z=%g, t=%g to t=%g \n",
          xmin, xmax, ymin, ymax, zmin, zmax, tmin, tmax);
  printf("\n");

  printf("Using trapezoidal method \n");
  np = 4; /* f1 */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("f1(x,y,z,t)=x+2*y+3*z+4*t+5 \n");
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (f1(x,y,z,tmin) + f1(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (f1(x,y,zmin,t) + f1(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f1(x,ymin,z,t) + f1(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f1(xmin,y,z,t) + f1(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + f1(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);

  np = 16; /* f1 */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("f1(x,y,z,t)=x+2*y+3*z+4*t+5 \n");
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (f1(x,y,z,tmin) + f1(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (f1(x,y,zmin,t) + f1(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f1(x,ymin,z,t) + f1(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f1(xmin,y,z,t) + f1(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + f1(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);

  np = 64; /* f1 */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("f1(x,y,z,t)=x+2*y+3*z+4*t+5 \n");
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (f1(x,y,z,tmin) + f1(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (f1(x,y,zmin,t) + f1(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f1(x,ymin,z,t) + f1(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f1(xmin,y,z,t) + f1(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + f1(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);

  printf("\n");


  np = 4; /* f2 */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("f2(x,y,z,t)=x*x+2*y*y+3*z*z+4*t*t+5*x*y+6*x*z+7*y*z+ \n");
  printf("            8*x*t+9*y*t+10*z*t+ \n");
  printf("            11*x+12*y+23*z+14*t+15 \n");
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (f2(x,y,z,tmin) + f2(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (f2(x,y,zmin,t) + f2(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f2(x,ymin,z,t) + f2(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f2(xmin,y,z,t) + f2(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + f2(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);

  np = 16; /* f2 */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("f2(x,y,z,t)=x*x+2*y*y+3*z*z+4*t*t+5*x*y+6*x*z+7*y*z+ \n");
  printf("            8*x*t+9*y*t+10*z*t+ \n");
  printf("            11*x+12*y+23*z+14*t+15 \n");
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (f2(x,y,z,tmin) + f2(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (f2(x,y,zmin,t) + f2(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f2(x,ymin,z,t) + f2(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f2(xmin,y,z,t) + f2(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + f2(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);

  np = 64; /* f2 */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("f2(x,y,z,t)=x*x+2*y*y+3*z*z+4*t*t+5*x*y+6*x*z+7*y*z+ \n");
  printf("            8*x*t+9*y*t+10*z*t+ \n");
  printf("            11*x+12*y+23*z+14*t+15 \n");
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (f2(x,y,z,tmin) + f2(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (f2(x,y,zmin,t) + f2(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f2(x,ymin,z,t) + f2(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f2(xmin,y,z,t) + f2(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + f2(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);

  printf("\n");


  np = 4; /* f3 */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("f3(x,y,z,t)=x*x*x+2*y*y*y+3*z*z*z+4*t*t*t+5*x*x*y+ \n");
  printf("            6*x*x*z+7*x*y*y+8*x*y*z+9*x*z*z+ \n");
  printf("            10*y*y*z+11*y*z*z+12*x*x*t+13*x*t*t+14*x*y*t+ \n");
  printf("            15*y*y*t+16*y*t*t+17*y*z*t+18*z*z*t+19*z*t*t+ \n");
  printf("            20*x*z*t+21*x*x+22*y*y+23*z*z+24*x*y+25*x*z+ \n");
  printf("            26*y*z+27*x*t+28*y*t+29*z*t+ \n");
  printf("            30*x+31*y+32*z+33*t+34 \n");
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (f3(x,y,z,tmin) + f3(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (f3(x,y,zmin,t) + f3(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f3(x,ymin,z,t) + f3(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f3(xmin,y,z,t) + f3(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + f3(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);

  np = 16; /* f3 */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("f3(x,y,z,t)=x*x*x+2*y*y*y+3*z*z*z+4*t*t*t+5*x*x*y+ \n");
  printf("            6*x*x*z+7*x*y*y+8*x*y*z+9*x*z*z+ \n");
  printf("            10*y*y*z+11*y*z*z+12*x*x*t+13*x*t*t+14*x*y*t+ \n");
  printf("            15*y*y*t+16*y*t*t+17*y*z*t+18*z*z*t+19*z*t*t+ \n");
  printf("            20*x*z*t+21*x*x+22*y*y+23*z*z+24*x*y+25*x*z+ \n");
  printf("            26*y*z+27*x*t+28*y*t+29*z*t+ \n");
  printf("            30*x+31*y+32*z+33*t+34 \n");
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (f3(x,y,z,tmin) + f3(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (f3(x,y,zmin,t) + f3(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f3(x,ymin,z,t) + f3(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f3(xmin,y,z,t) + f3(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + f3(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);

  np = 64; /* f3 */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("f3(x,y,z,t)=x*x*x+2*y*y*y+3*z*z*z+4*t*t*t+5*x*x*y+ \n");
  printf("            6*x*x*z+7*x*y*y+8*x*y*z+9*x*z*z+ \n");
  printf("            10*y*y*z+11*y*z*z+12*x*x*t+13*x*t*t+14*x*y*t+ \n");
  printf("            15*y*y*t+16*y*t*t+17*y*z*t+18*z*z*t+19*z*t*t+ \n");
  printf("            20*x*z*t+21*x*x+22*y*y+23*z*z+24*x*y+25*x*z+ \n");
  printf("            26*y*z+27*x*t+28*y*t+29*z*t+ \n");
  printf("            30*x+31*y+32*z+33*t+34 \n");
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (f3(x,y,z,tmin) + f3(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (f3(x,y,zmin,t) + f3(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f3(x,ymin,z,t) + f3(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (f3(xmin,y,z,t) + f3(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + f3(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);

  printf("\n");

  np = 4; /* fe */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)*exp(t)    np=%d \n", np);
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (fe(x,y,z,tmin) + fe(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (fe(x,y,zmin,t) + fe(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (fe(x,ymin,z,t) + fe(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (fe(xmin,y,z,t) + fe(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + fe(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);

  np = 16; /* fe */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)*exp(t)    np=%d \n", np);
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (fe(x,y,z,tmin) + fe(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (fe(x,y,zmin,t) + fe(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (fe(x,ymin,z,t) + fe(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (fe(xmin,y,z,t) + fe(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + fe(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);

  np = 64; /* fe */
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  ht = (tmax-tmin)/(double)(np-1);
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)*exp(t)    np=%d \n", np);
  printf("np=%d,  hx=%g, hy=%g, hz=%g, ht=%g \n", np, hx, hy, hz, ht);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(k=0; k<np; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + (fe(x,y,z,tmin) + fe(x,y,z,tmax))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;  
        volume = volume + (fe(x,y,zmin,t) + fe(x,y,zmax,t))/2.0;
      }
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (fe(x,ymin,z,t) + fe(x,ymax,z,t))/2.0;
      }
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      for(m=0; m<np; m++)
      {
        t = tmin+(double)m*ht;
        volume = volume + (fe(xmin,y,z,t) + fe(xmax,y,z,t))/2.0;
      }
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        for(m=1; m<np-1; m++)
	{
          t = tmin+(double)m*ht;
          volume = volume + fe(x,y,z,t);
	}
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*(tmax-tmin)*volume/
           (double)((np-1)*(np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);

  printf("\n");


  printf("Using Gauss Legendre method \n");

  np = 4;
  printf("Gauss Legendre  np=%d \n", np);
  printf("f1(x,y,z,t)=x+2*y+3*z+4*t+5; \n");
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  gaulegf(tmin, tmax, tt, wt, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        for(m=1; m<=np; m++)
	{
          volume = volume + wx[i]*wy[j]*wz[k]*wt[m]*
                        f1(xx[i], yy[j], zz[k], tt[m]);
	}
      }
    }
  }
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);
  printf("\n");

  np = 4;
  printf("Gauss Legendre  np=%d \n", np);
  printf("f2(x,y,z,t)=x*x+2*y*y+3*z*z+4*t*t+5*x*y+6*x*z+7*y*z+ \n");
  printf("            8*x*t+9*y*t+10*z*t+ \n");
  printf("            11*x+12*y+23*z+14*t+15 \n");
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  gaulegf(tmin, tmax, tt, wt, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        for(m=1; m<=np; m++)
	{
          volume = volume + wx[i]*wy[j]*wz[k]*wt[m]*
                        f2(xx[i], yy[j], zz[k], tt[m]);
	}
      }
    }
  }
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);
  printf("\n");

  np = 4;
  printf("Gauss Legendre  np=%d \n", np);
  printf("f3(x,y,z,t)=x*x*x+2*y*y*y+3*z*z*z+4*t*t*t+5*x*x*y+ \n");
  printf("            6*x*x*z+7*x*y*y+8*x*y*z+9*x*z*z+ \n");
  printf("            10*y*y*z+11*y*z*z+12*x*x*t+13*x*t*t+14*x*y*t+ \n");
  printf("            15*y*y*t+16*y*t*t+17*y*z*t+18*z*z*t+19*z*t*t+ \n");
  printf("            20*x*z*t+21*x*x+22*y*y+23*z*z+24*x*y+25*x*z+ \n");
  printf("            26*y*z+27*x*t+28*y*t+29*z*t+ \n");
  printf("            30*x+31*y+32*z+33*t+34 \n");
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  gaulegf(tmin, tmax, tt, wt, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        for(m=1; m<=np; m++)
	{
          volume = volume + wx[i]*wy[j]*wz[k]*wt[m]*
                        f3(xx[i], yy[j], zz[k], tt[m]);
	}
      }
    }
  }
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);

  np = 8;
  printf("Gauss Legendre  np=%d \n", np);
  printf("f3(x,y,z,t)=x*x*x+2*y*y*y+3*z*z*z+4*t*t*t+5*x*x*y+ \n");
  printf("            6*x*x*z+7*x*y*y+8*x*y*z+9*x*z*z+ \n");
  printf("            10*y*y*z+11*y*z*z+12*x*x*t+13*x*t*t+14*x*y*t+ \n");
  printf("            15*y*y*t+16*y*t*t+17*y*z*t+18*z*z*t+19*z*t*t+ \n");
  printf("            20*x*z*t+21*x*x+22*y*y+23*z*z+24*x*y+25*x*z+ \n");
  printf("            26*y*z+27*x*t+28*y*t+29*z*t+ \n");
  printf("            30*x+31*y+32*z+33*t+34 \n");
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  gaulegf(tmin, tmax, tt, wt, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        for(m=1; m<=np; m++)
	{
          volume = volume + wx[i]*wy[j]*wz[k]*wt[m]*
                        f3(xx[i], yy[j], zz[k], tt[m]);
	}
      }
    }
  }
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);
  printf("\n");

  printf("\n");

  np = 4;
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)*exp(t)    np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  gaulegf(tmin, tmax, tt, wt, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        for(m=1; m<=np; m++)
	{
          volume = volume + wx[i]*wy[j]*wz[k]*wt[m]*
                        fe(xx[i], yy[j], zz[k], tt[m]);
	}
      }
    }
  }
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);

  np = 8;
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)*exp(t)    np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  gaulegf(tmin, tmax, tt, wt, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        for(m=1; m<=np; m++)
	{
          volume = volume + wx[i]*wy[j]*wz[k]*wt[m]*
                        fe(xx[i], yy[j], zz[k], tt[m]);
	}
      }
    }
  }
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);

  np = 16;
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)*exp(t)    np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  gaulegf(tmin, tmax, tt, wt, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        for(m=1; m<=np; m++)
	{
          volume = volume + wx[i]*wy[j]*wz[k]*wt[m]*
                        fe(xx[i], yy[j], zz[k], tt[m]);
	}
      }
    }
  }
  printf("%d function evaluations \n", np*np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);
  printf("\n");

  return 0;
  } /* end main of fourd_int.c */

static void gaulegf(double x1, double x2, double x[], double w[], int n)
{
  int i, j, m;
  double eps = 3.0E-14;
  double p1, p2, p3, pp, xl, xm, z, z1;
  
  m = (n+1)/2;
  xm = 0.5*(x2+x1);
  xl = 0.5*(x2-x1);
  for(i=1; i<=m; i++)
  {
    z = cos(3.141592654*((double)i-0.25)/((double)n+0.5));
    while(1)
    {
      p1 = 1.0;
      p2 = 0.0;
      for(j=1; j<=n; j++)
      {
        p3 = p2;
        p2 = p1;
        p1 = ((2.0*(double)j-1.0)*z*p2-((double)j-1.0)*p3)/
              (double)j;
      }
      pp = (double)n*(z*p1-p2)/(z*z-1.0);
      z1 = z;
      z = z1 - p1/pp;
      if(abs(z-z1) <= eps) break;
    }
    x[i] = xm - xl*z;
    x[n+1-i] = xm + xl*z;
    w[i] = 2.0*xl/((1.0-z*z)*pp*pp);
    w[n+1-i] = w[i];
  }
} /* end gaulegf */
/* end fourd_int.c */