/* test_nuderiv3dg.c   test  nuderiv3dg  with max 37 points            */
/*                     order 1  minimum independent points 4           */
/*                     order 2  minimum independent points 10          */
/*                     order 3  minimum independent points 20          */
/*                     order 4  minimum independent points 35          */
/*                     going beyond minimum points does not help much  */
/*                     except points must be linearly independent      */
/*                     (automatically selected from extra points)      */
#include <stdio.h>
#include <stdlib.h>
#include "nuderiv3dg.h"
#include "inverse.h"

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef  max
#define max(x,y) ((x)>(y)?(x):(y))

/* 1st order test case */

static double f1(double x, double y, double z)
{
  return x + 2.0*y + 3.0*z + 4.0;
}
static double f1x(double x, double y, double z) /* df1/dx */
{
  return 1.0;
}
static double f1y(double x, double y, double z) /* df1/dy */
{
  return 2.0;
}
static double f1z(double x, double y, double z) /* df1/dz */
{
  return 3.0;
}

/* 2nd order test case */

static double f2(double x, double y, double z)
{
  return x*x + 2.0*x*y + 3.0*x*z + 4.0*y*y +
     5.0*y*z + 6.0*z*z + 7.0*x   + 8.0*y   +
    9.0*z   + 10.0;
}
static double f2x(double x, double y, double z) /* df2/dx */
{
  return 2.0*x + 2.0*y + 3.0*z + 7.0;
}
static double f2y(double x, double y, double z)
{
  return 2.0*x + 8.0*y + 5.0*z + 8.0;
}
static double f2z(double x, double y, double z)
{
  return 3.0*x + 5.0*y + 12.0*z + 9.0;
}
static double f2xx(double x, double y, double z) /* df2^2/dx^2 */
{
  return 2.0;
}
static double f2xy(double x, double y, double z) /* df2^2/dxdy */
{
  return 2.0;
}
static double f2xz(double x, double y, double z)
{
  return 3.0;
}
static double f2yy(double x, double y, double z)
{
  return 8.0;
}
static double f2yz(double x, double y, double z)
{
  return 5.0;
}
static double f2zz(double x, double y, double z)
{
  return 12.0;
}


/* 3rd order test case */
static double f3(double x, double y, double z)
{
  return x*x*x +  2.0*x*x*y +  3.0*x*x*z + 4.0*x*y*y +
     5.0*x*y*z +  6.0*x*z*z +  7.0*y*y*y + 8.0*y*y*z +
     9.0*y*z*z + 10.0*z*z*z + 11.0*x*x   + 12.0*x*y +
    13.0*x*z   + 14.0*y*y   + 15.0*y*z   + 16.0*z*z  +
    17.0*x     + 18.0*y     + 19.0*z     + 20.0;
}
static double f3x(double x, double y, double z)
{
  return 3.0*x*x +  4.0*x*y +  6.0*x*z +  4.0*y*y +
         5.0*y*z +  6.0*z*z +
                              22.0*x   + 12.0*y   +
        13.0*z  +
    17.0;
}
static double f3y(double x, double y, double z)
{
  return             2.0*x*x  +             8.0*x*y +
         5.0*x*z +              21.0*y*y + 16.0*y*z +
         9.0*z*z +                         12.0*x   +
                    28.0*y    + 15.0*z   +
    18.0;
}
static double f3z(double x, double y, double z)
{
  return                       3.0*x*x +
         5.0*x*y + 12.0*x*z +             8.0*y*y +
        18.0*y*z + 30*z*z   +
        13.0*x   +            15.0*y   + 32.0*z   +
    19.0;
}
static double f3xx(double x, double y, double z)
{
  return 6.0*x + 4.0*y +  6.0*z + 22.0;
}
static double f3xy(double x, double y, double z)
{
  return 4.0*x + 8.0*y + 5.0*z + 12.0  ;
}
static double f3xz(double x, double y, double z)
{
  return 6.0*x + 5.0*y + 12.0*z + 13.0;
}
static double f3yy(double x, double y, double z)
{
  return 8.0*x + 42.0*y + 16.0*z + 28.0;
}
static double f3yz(double x, double y, double z)
{
  return 5.0*x + 16.0*y + 18.0*z + 15.0;
}
static double f3zz(double x, double y, double z)
{
  return 12.0*x + 18.0*y + 60.0*z + 32.0;
}
static double f3xxx(double x, double y, double z)
{
  return 6.0;
}
static double f3xxy(double x, double y, double z)
{
  return 4.0;
}
static double f3xxz(double x, double y, double z)
{
  return 6.0;
}
static double f3xyy(double x, double y, double z)
{
  return 8.0;
}
static double f3xyz(double x, double y, double z)
{
  return 5.0;
}
static double f3xzz(double x, double y, double z)
{
  return 12.0;
}
static double f3yyy(double x, double y, double z)
{
  return 42.0;
}
static double f3yyz(double x, double y, double z)
{
  return 16.0;
}
static double f3yzz(double x, double y, double z)
{
  return 18.0;
}
static double f3zzz(double x, double y, double z)
{
  return 60.0;
}

/* 4th order test case  */
static double f4(double x, double y, double z)
{
  return x*x*x*x +  2.0*x*x*x*y +  3.0*x*x*x*z +
     4.0*x*x*y*y +  5.0*x*x*y*z +  6.0*x*x*z*z +
     7.0*x*y*y*y +  8.0*x*y*y*z +  9.0*x*y*z*z +
    10.0*x*z*z*z + 11.0*y*y*y*y + 12.0*y*y*y*z +
    13.0*y*y*z*z + 14.0*y*z*z*z + 15.0*z*z*z*z +
    16.0*x*x*x   + 17.0*x*x*y   + 18.0*x*x*z   + 19.0*x*y*y +
    20.0*x*y*z   + 21.0*x*z*z   + 22.0*y*y*y   + 23.0*y*y*z +
    24.0*y*z*z   + 25.0*z*z*z   + 26.0*x*x     + 27.0*x*y +
    28.0*x*z     + 29.0*y*y     + 30.0*y*z     + 31.0*z*z  +
    32.0*x       + 33.0*y       + 34.0*z       + 35.0;
}
static double f4x(double x, double y, double z) /* df4/dx */
{
  return 4.0*x*x*x +  6.0*x*x*y +  9.0*x*x*z +
         8.0*x*y*y + 10.0*x*y*z + 12.0*x*z*z +
         7.0*y*y*y +  8.0*y*y*z +  9.0*y*z*z +
        10.0*z*z*z + 48.0*x*x   + 34.0*x*y   +
        36.0*x*z   + 19.0*y*y   + 20.0*y*z   +
        21.0*z*z   + 52.0*x     + 27.0*y     +
        28.0*z     + 32.0;
}
static double f4y(double x, double y, double z) /* df4/dx */
{
  return  2.0*x*x*x +  8.0*x*x*y +  5.0*x*x*z + 
         21.0*x*y*y + 16.0*x*y*z +  9.0*x*z*z +
         44.0*y*y*y + 36.0*y*y*z + 26.0*y*z*z +
         14.0*z*z*z + 17.0*x*x   + 38.0*x*y   +
         20.0*x*z   + 66.0*y*y   + 46.0*y*z   +
         24.0*z*z   + 27.0*x     + 58.0*y     +
         30.0*z     + 33.0;
}
static double f4z(double x, double y, double z)
{
  return 3.0*x*x*x +  5.0*x*x*y + 12.0*x*x*z +
         8.0*x*y*y + 18.0*x*y*z + 30.0*x*z*z +
        12.0*y*y*y + 26.0*y*y*z + 42.0*y*z*z +
        60.0*z*z*z + 18.0*x*x   + 20.0*x*y   +
        42.0*x*z   + 23.0*y*y   + 48.0*y*z   +
        75.0*z*z   + 28.0*x     +   30.0*y   +
        62.0*z     + 34.0;
}
static double f4xx(double x, double y, double z) /* df4/dxx */
{
  return 12.0*x*x + 12.0*x*y + 18.0*x*z +
          8.0*y*y + 10.0*y*z + 12.0*z*z  +
         96.0*x   + 34.0*y   + 36.0*z + 52.0;
}
static double f4xy(double x, double y, double z) /* df4/dxy */
{
  return  6.0*x*x +  16.0*x*y + 10.0*x*z + 
         21.0*y*y +  16.0*y*z +  9.0*z*z +
         34.0*x   +  38.0*y   + 20.0*z   + 27.0;
}
static double f4xz(double x, double y, double z) /* df4/dxz */
{
  return 9.0*x*x + 10.0*x*y + 24.0*x*z +
         8.0*y*y + 18.0*y*z + 30.0*z*z +
        36.0*x   + 20.0*y   + 42.0*z   + 28.0;
}
static double f4yy(double x, double y, double z) /* df4/dyy */
{
  return  8.0*x*x +  42.0*x*y + 16.0*x*z +
        132.0*y*y +  72.0*y*z + 26.0*z*z + 
         38.0*x   + 132.0*y   + 46.0*z   + 58.0;
}
static double f4yz(double x, double y, double z) /* df4/dyz */
{
  return  5.0*x*x + 16.0*x*y + 18.0*x*z +
         36.0*y*y + 52.0*y*z + 42.0*z*z +
         20.0*x   + 46.0*y   + 48.0*z   + 30.0;
}
static double f4zz(double x, double y, double z)
{
  return  12.0*x*x + 18.0*x*y +  60.0*x*z +
          26.0*y*y + 84.0*y*z + 180.0*z*z +
          42.0*x   + 48.0*y   + 150.0*z   + 62.0;
}
static double f4xxx(double x, double y, double z) /* df4/dxxx */
{
  return 24.0*x + 12.0*y + 18.0*z + 96.0;
}
static double f4xxy(double x, double y, double z) /* df4/dxxy */
{
  return 12.0*x + 16.0*y + 10.0*z + 34.0;
}
static double f4xxz(double x, double y, double z) /* df4/dxxz */
{
  return 18.0*x + 10.0*y + 24.0*z + 36.0;
}
static double f4xyy(double x, double y, double z) /* df4/dxyy */
{
  return  16.0*x + 42.0*y + 16.0*z + 38.0;
}
static double f4xyz(double x, double y, double z) /* df4/dxyz */
{
  return 10.0*x + 16.0*y +  18.0*z + 20.0;
}
static double f4xzz(double x, double y, double z) /* df4/dxzz */
{
  return  24.0*x + 18.0*y +  60.0*z + 42.0;
}
static double f4yyy(double x, double y, double z) /* df4/dyyy */
{
  return 42.0*x + 264.0*y + 72.0*z + 132.0;
}
static double f4yyz(double x, double y, double z) /* df4/dyyz */
{
  return  16.0*x + 72.0*y + 52.0*z + 46.0;
}
static double f4yzz(double x, double y, double z) /* df4/dyzz */
{
  return  18.0*x + 52.0*y + 84.0*z + 48.0;
}
static double f4zzz(double x, double y, double z) /* df4/dzzz */
{
  return   60.0*x + 84.0*y + 360.0*z + 150.0;
}
static double f4xxxx(double x, double y, double z) /* df4/dxxxx */
{
  return 24.0;
}
static double f4xxxy(double x, double y, double z) /* df4/dxxxy */
{
  return 12.0;
}
static double f4xxxz(double x, double y, double z) /* df4/dxxxz */
{
  return 18.0;
}
static double f4xxyy(double x, double y, double z) /* df4/dxxyy */
{
  return 16.0;
}
static double f4xxyz(double x, double y, double z) /* df4/dxxyz */
{
  return 10.0;
}
static double f4xxzz(double x, double y, double z) /* df4/dxxzz */
{
  return  24.0;
}
static double f4xyyy(double x, double y, double z) /* df4/dxyyy */
{
  return  42.0;
}
static double f4xyyz(double x, double y, double z) /* df4/dxyyz */
{
  return  16.0;
}
static double f4xyzz(double x, double y, double z) /* df4/dxyzz */
{
  return  18.0;
}
static double f4xzzz(double x, double y, double z) /* df4/dxzzz */
{
  return  60.0;
}
static double f4yyyy(double x, double y, double z)  /* df4/dyyyy */
{
  return 264.0;
}
static double f4yyyz(double x, double y, double z)  /* df4/dyyyz */
{
  return  72.0;
}
static double f4yyzz(double x, double y, double z)  /* df4/dyyzz */
{
  return  52.0;
}
static double f4yzzz(double x, double y, double z)  /* df4/dyzzz */
{
  return  84.0;
}
static double f4zzzz(double x, double y, double z)/* df4/dzzzz */
{
  return 360.0;
}


int main(int argc, char *argv[])
{
  int i, j;
  int debug = 0;
  int n, npoint;
  int order; // minimum independent npoints  1, 4, 10, 20, 35
  double x[38] = {0.05, 0.14, 0.27, 0.42, 0.57, 0.71, 0.86,
                  0.04, 0.18, 0.31, 0.43, 0.63, 0.82, 0.95,
                  0.07, 0.16, 0.33, 0.41, 0.59, 0.82, 0.99,
                  0.11, 0.21, 0.35, 0.45, 0.65, 0.75, 0.97,
                  0.06, 0.17, 0.32, 0.46, 0.61, 0.79, 0.83, 0.87, 0.93, 0.98};
  double y[38] = {0.06, 0.05, 0.07, 0.10, 0.13, 0.15, 0.17,
                  0.30, 0.33, 0.36, 0.39, 0.31, 0.32, 0.37,
                  0.48, 0.49, 0.52, 0.51, 0.47, 0.50, 0.53,
                  0.65, 0.77, 0.71, 0.67, 0.67, 0.75, 0.76,
                  0.92, 0.97, 0.93, 0.98, 0.87, 0.91, 0.94, 0.88, 0.77, 0.65};
  double z[38] = {0.07, 0.28, 0.47, 0.69, 0.91,
                  0.15, 0.45, 0.68, 0.93, 0.05,
                  0.26, 0.50, 0.97, 0.08, 0.14,
                  0.39, 0.65, 0.76, 0.93, 0.27,
                  0.56, 0.82, 0.99, 0.11, 0.33,
                  0.72, 0.94, 0.09, 0.22, 0.65,
                  0.95, 0.10, 0.24, 0.42, 0.51, 0.82, 0.93, 0.99};
  double U[38];
  int ip[38];
  double est, actual, dx, dy, dz;
  double cx[38], cy[38], cz[38], cxx[38], cxy[38], cxz[38];
  double cyy[38], cyz[38], czz[38], cxxx[38], cxxy[38];
  double cxxz[38], cxyy[38], cxyz[38], cxzz[38], cyyy[38], cyyz[38];
  double cyzz[38], czzz[38], cxxxx[38], cxxxy[38], cxxxz[38];
  double cxxyy[38], cxxyz[38], cxxzz[38], cxyyy[38], cxyyz[38];
  double cxyzz[38], cxzzz[38], cyyyy[38], cyyyz[38], cyyzz[38];
  double cyzzz[38], czzzz[38];
  double error, maxerr;
  int point;
  int m = 0;

  printf("test_nuderiv3dg.c running \n");
  self_test(); /* nuderiv3dg.c */
  n = 36; // minimum for order=4, can use 36 or 37
  npoint = 6; // minimum 4 independent for order 1
  for(i=0; i<npoint; i++)
  {
    U[i] = f1(x[i],y[i],z[i]);
    dx = f1x(x[i],y[i],z[i]);
    dy = f1y(x[i],y[i],z[i]);
    dz = f1z(x[i],y[i],z[i]);
    if(debug)
      printf("point %2d  x=%5.3f, y=%5.3f, z=%5.3f, U=%7.3f, dx=%f, dy=%f, dz=%f \n",
	     i, x[i], y[i], z[i], U[i], dx, dy, dz);
  } 

  /* test f1 dx, dy, dz */
  order = 1;
  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f1 dx */
  {
    m = nu3dx(order, npoint, point, x, y, z, cx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cx[i]*U[ip[i]];
    actual = f1x(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f1 dx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f1 dx */
  printf("f1 dx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f1 dy */
  {
    m = nu3dy(order, npoint, point, x, y, z, cy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cy[i]*U[ip[i]];
    actual = f1y(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug)printf("f1 dy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f1 dy */
  printf("f1 dy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f1 dz */
  {
    m = nu3dz(order, npoint, point, x, y, z, cz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cz[i]*U[ip[i]];
    actual = f1z(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f1 dz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f1 dz */
  printf("f1 dz max error=%g \n", maxerr);

  /* test f2 dx, dy, dz, dxx, dxy, dxz, dyy, dyz, dzz */
  npoint = 11; // minimum 10 independent
  for(i=0; i<npoint; i++)  U[i] = f2(x[i],y[i],z[i]);
  order = 2;
  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dx */
  {
    m = nu3dx(order, npoint, point, x, y, z, cx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cx[i]*U[ip[i]];
    actual = f2x(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f2 dx */
  printf("f2 dx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dy */
  {
    m = nu3dy(order, npoint, point, x, y, z, cy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cy[i]*U[ip[i]];
    actual = f2y(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f2 dy */
  printf("f2 dy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dz */
  {
    m = nu3dz(order, npoint, point, x, y, z, cz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cz[i]*U[ip[i]];
    actual = f2z(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f2 dz */
  printf("f2 dz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dxx */
  {
    m = nu3dxx(order, npoint, point, x, y, z, cxx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxx[i]*U[ip[i]];
    actual = f2xx(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dxx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f2 dxx */
  printf("f2 dxx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dxy */
  {
    m = nu3dxy(order, npoint, point, x, y, z, cxy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxy[i]*U[ip[i]];
    actual = f2xy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dxy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f2 dxy */
  printf("f2 dxy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dxz */
  {
    m = nu3dxz(order, npoint, point, x, y, z, cxz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxz[i]*U[ip[i]];
    actual = f2xz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dxz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f2 dxz */
  printf("f2 dxz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dyy */
  {
    m = nu3dyy(order, npoint, point, x, y, z, cyy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyy[i]*U[ip[i]];
    actual = f2yy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dyy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f2 dyy */
  printf("f2 dyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dyz */
  {
    m = nu3dyz(order, npoint, point, x, y, z, cyz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyz[i]*U[ip[i]];
    actual = f2yz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dyz %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f2 dyz */
  printf("f2 dyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f2 dzz */
  {
    m = nu3dzz(order, npoint, point, x, y, z, czz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += czz[i]*U[ip[i]];
    actual = f2zz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f2 dzz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f2 dzz */
  printf("f2 dzz max error=%g \n", maxerr);


  /* test f3 dx, dy, dz, dxx, dxy, dxz, dyy, dyz, dzz */
  /*         dxxx, dxxy, dxxz, dxyy, dxyz, dxzz, dyyy, dyyz, dyzz, dzzz */
  npoint = 21; // minimum 20 independent
  for(i=0; i<npoint; i++)  U[i] = f3(x[i],y[i],z[i]);
  order = 3;
  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dx */
  {
    m = nu3dx(order, npoint, point, x, y, z, cx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cx[i]*U[ip[i]];
    actual = f3x(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dx */
  printf("f3 dx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dy */
  {
    m = nu3dy(order, npoint, point, x, y, z, cy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cy[i]*U[ip[i]];
    actual = f3y(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f3 dy */
  printf("f3 dy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dz */
  {
    m = nu3dz(order, npoint, point, x, y, z, cz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cz[i]*U[ip[i]];
    actual = f3z(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f3 dz */
  printf("f3 dz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxx */
  {
    m = nu3dxx(order, npoint, point, x, y, z, cxx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxx[i]*U[ip[i]];
    actual = f3xx(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxx */
  printf("f3 dxx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxy */
  {
    m = nu3dxy(order, npoint, point, x, y, z, cxy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxy[i]*U[ip[i]];
    actual = f3xy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxy */
  printf("f3 dxy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxz */
  {
    m = nu3dxz(order, npoint, point, x, y, z, cxz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxz[i]*U[ip[i]];
    actual = f3xz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxz */
  printf("f3 dxz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dyy */
  {
    m = nu3dyy(order, npoint, point, x, y, z, cyy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyy[i]*U[ip[i]];
    actual = f3yy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dyy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f3 dyy */
  printf("f3 dyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dyz */
  {
    m = nu3dyz(order, npoint, point, x, y, z, cyz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyz[i]*U[ip[i]];
    actual = f3yz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dyz %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f3 dyz */
  printf("f3 dyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dzz */
  {
    m = nu3dzz(order, npoint, point, x, y, z, czz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += czz[i]*U[ip[i]];
    actual = f3zz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dzz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f3 dzz */
  printf("f3 dzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxxx */
  {
    m = nu3dxxx(order, npoint, point, x, y, z, cxxx, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxx[i]*U[ip[i]];
    actual = f3xxx(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxxx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxxx */
  printf("f3 dxxx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxxy */
  {
    m = nu3dxxy(order, npoint, point, x, y, z, cxxy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxy[i]*U[ip[i]];
    actual = f3xxy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxxy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxxy */
  printf("f3 dxxy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxxz */
  {
    m = nu3dxxz(order, npoint, point, x, y, z, cxxz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxz[i]*U[ip[i]];
    actual = f3xxz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxxz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxxz */
  printf("f3 dxxz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxyy */
  {
    m = nu3dxyy(order, npoint, point, x, y, z, cxyy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxyy[i]*U[ip[i]];
    actual = f3xyy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxyy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxyy */
  printf("f3 dxyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxyz */
  {
    m = nu3dxyz(order, npoint, point, x, y, z, cxyz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxyz[i]*U[ip[i]];
    actual = f3xyz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxyz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f3 dxyz */
  printf("f3 dxyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dxzz */
  {
    m = nu3dxzz(order, npoint, point, x, y, z, cxzz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cxzz[i]*U[ip[i]];
    actual = f3xzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dxzz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f3 dxzz */
  printf("f3 dxzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dyyy */
  {
    m = nu3dyyy(order, npoint, point, x, y, z, cyyy, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyyy[i]*U[ip[i]];
    actual = f3yyy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dyy %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f3 dyyy */
  printf("f3 dyyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dyyz */
  {
    m = nu3dyyz(order, npoint, point, x, y, z, cyyz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyyz[i]*U[ip[i]];
    actual = f3yyz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dyyz %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f3 dyyz */
  printf("f3 dyyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dyzz */
  {
    m = nu3dyzz(order, npoint, point, x, y, z, cyzz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += cyzz[i]*U[ip[i]];
    actual = f3yzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dyzz %2d, est=%g, actual=%g, error=%g \n",
           point, est, actual, error);
  } /* end f3 dyzz */
  printf("f3 dyzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f3 dzzz */
  {
    m = nu3dzzz(order, npoint, point, x, y, z, czzz, ip);
    if(debug) printf("order=%1d, n=%2d, m=%2d, point=%2d \n",
                     order, npoint, m, point);
    est = 0.0;
    for(i=0; i<m; i++) est += czzz[i]*U[ip[i]];
    actual = f3zzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f3 dzzz %2d, est=%g, actual=%g, error=%g \n",
              point, est, actual, error);
  } /* end f3 dzzz */
  printf("f3 dzzz max error=%g \n", maxerr);

  /* test f4 dx, dy, dz, dxx, dxy, dxz, dyy, dyz, dzz */
  /*         dxxx, dxxy, dxxz, dxyy, dxyz, dxzz, dyyy, dyyz, dyzz, dzzz */
  /*         dxxxx, dxxxy, dxxxz, dxxyy, dxxyz, dxxzz, dxyyy, dxyyz     */
  /*         dxyzz, dxzzz, dyyyy, dyyyz, dyyzz, dyzzz, dzzzz            */
  npoint = 36; // minimum 35 independent
  for(i=0; i<npoint; i++)  U[i] = f4(x[i],y[i],z[i]);
  order = 4; // order of fit, not derivative
  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dx */
  {
    m = nu3dx(order, npoint, point, x, y, z, cx, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cx[i]*U[ip[i]];
    actual = f4x(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dx */
  printf("f4 dx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dy */
  {
    m = nu3dy(order, npoint, point, x, y, z, cy, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cy[i]*U[ip[i]];
    actual = f4y(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dy */
  printf("f4 dy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dz */
  {
    m = nu3dz(order, npoint, point, x, y, z, cz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cz[i]*U[ip[i]];
    actual = f4z(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dz */
  printf("f4 dz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxx */
  {
    m = nu3dxx(order, npoint, point, x, y, z, cxx, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxx[i]*U[ip[i]];
    actual = f4xx(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxx */
  printf("f4 dxx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxy */
  {
    m = nu3dxy(order, npoint, point, x, y, z, cxy, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxy[i]*U[ip[i]];
    actual = f4xy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxy */
  printf("f4 dxy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxz */
  {
    m = nu3dxz(order, npoint, point, x, y, z, cxz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxz[i]*U[ip[i]];
    actual = f4xz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxz */
  printf("f4 dxz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dyy */
  {
    m = nu3dyy(order, npoint, point, x, y, z, cyy, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cyy[i]*U[ip[i]];
    actual = f4yy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dyy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dyy */
  printf("f4 dyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dyz */
  {
    m = nu3dyz(order, npoint, point, x, y, z, cyz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cyz[i]*U[ip[i]];
    actual = f4yz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dyz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dyz */
  printf("f4 dyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dzz */
  {
    m = nu3dzz(order, npoint, point, x, y, z, czz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += czz[i]*U[ip[i]];
    actual = f4zz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dzz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dzz */
  printf("f4 dzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxxx */
  {
    m = nu3dxxx(order, npoint, point, x, y, z, cxxx, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxx[i]*U[ip[i]];
    actual = f4xxx(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxxx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxxx */
  printf("f4 dxxx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxxy */
  {
    m = nu3dxxy(order, npoint, point, x, y, z, cxxy, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxy[i]*U[ip[i]];
    actual = f4xxy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxxy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxxy */
  printf("f4 dxxy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxxz */
  {
    m = nu3dxxz(order, npoint, point, x, y, z, cxxz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxz[i]*U[ip[i]];
    actual = f4xxz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxxz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxxz */
  printf("f4 dxxz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxyy */
  {
    m = nu3dxyy(order, npoint, point, x, y, z, cxyy, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxyy[i]*U[ip[i]];
    actual = f4xyy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxyy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxyy */
  printf("f4 dxyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxyz */
  {
    m = nu3dxyz(order, npoint, point, x, y, z, cxyz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxyz[i]*U[ip[i]];
    actual = f4xyz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxyz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxyz */
  printf("f4 dxyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxzz */
  {
    m = nu3dxzz(order, npoint, point, x, y, z, cxzz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxyz[i]*U[ip[i]];
    actual = f4xzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxzz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxzz */
  printf("f4 dxzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dyyy */
  {
    m = nu3dyyy(order, npoint, point, x, y, z, cyyy, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cyyy[i]*U[ip[i]];
    actual = f4yyy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dyyy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dyyy */
  printf("f4 dyyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dyyz */
  {
    m = nu3dyyz(order, npoint, point, x, y, z, cyyz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cyyz[i]*U[ip[i]];
    actual = f4yyz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dyyz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dyyz */
  printf("f4 dyyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dyzz */
  {
    m = nu3dyzz(order, npoint, point, x, y, z, cyzz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cyzz[i]*U[ip[i]];
    actual = f4yzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dyzz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dyzz */
  printf("f4 dyzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dzzz */
  {
    m = nu3dzzz(order, npoint, point, x, y, z, czzz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += czzz[i]*U[ip[i]];
    actual = f4zzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dzzz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dzzz */
  printf("f4 dzzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxxxx */
  {
    m = nu3dxxxx(order, npoint, point, x, y, z, cxxxx, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxxx[i]*U[ip[i]];
    actual = f4xxxx(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxxxx %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxxxx */
  printf("f4 dxxxx max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxxxy */
  {
    m = nu3dxxxy(order, npoint, point, x, y, z, cxxxy, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxxy[i]*U[ip[i]];
    actual = f4xxxy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxxxy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxxxy */
  printf("f4 dxxxy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxxxz */
  {
    m = nu3dxxxz(order, npoint, point, x, y, z, cxxxz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxxz[i]*U[ip[i]];
    actual = f4xxxz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxxxz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxxxz */
  printf("f4 dxxxz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxxyy */
  {
    m = nu3dxxyy(order, npoint, point, x, y, z, cxxyy, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxyy[i]*U[ip[i]];
    actual = f4xxyy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxxyy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxxyy */
  printf("f4 dxxyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxxyy */
  {
    m = nu3dxxyz(order, npoint, point, x, y, z, cxxyz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxyz[i]*U[ip[i]];
    actual = f4xxyz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxxyz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxxyz */
  printf("f4 dxxyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxxzz */
  {
    m = nu3dxxzz(order, npoint, point, x, y, z, cxxzz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxxzz[i]*U[ip[i]];
    actual = f4xxzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxxzz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxxzz */
  printf("f4 dxxzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxyyy */
  {
    m = nu3dxyyy(order, npoint, point, x, y, z, cxyyy, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxyyy[i]*U[ip[i]];
    actual = f4xyyy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxyyy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxyyy */
  printf("f4 dxyyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxyyz */
  {
    m = nu3dxyyz(order, npoint, point, x, y, z, cxyyz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxyyz[i]*U[ip[i]];
    actual = f4xyyz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxyyz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxyyz */
  printf("f4 dxyyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxyzz */
  {
    m = nu3dxyzz(order, npoint, point, x, y, z, cxyzz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxyzz[i]*U[ip[i]];
    actual = f4xyzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxyzz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxyzz */
  printf("f4 dxyzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dxzzz */
  {
    m = nu3dxzzz(order, npoint, point, x, y, z, cxzzz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cxzzz[i]*U[ip[i]];
    actual = f4xzzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dxzzz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dxzzz */
  printf("f4 dxzzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dyyyy */
  {
    m = nu3dyyyy(order, npoint, point, x, y, z, cyyyy, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cyyyy[i]*U[ip[i]];
    actual = f4yyyy(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dyyyy %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dyyyy */
  printf("f4 dyyyy max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dyyyz */
  {
    m = nu3dyyyz(order, npoint, point, x, y, z, cyyyz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cyyyz[i]*U[ip[i]];
    actual = f4yyyz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dyyyz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dyyyz */
  printf("f4 dyyyz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dyyzz */
  {
    m = nu3dyyzz(order, npoint, point, x, y, z, cyyzz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cyyzz[i]*U[ip[i]];
    actual = f4yyzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dyyzz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dyyzz */
  printf("f4 dyyzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dyzzz */
  {
    m = nu3dyzzz(order, npoint, point, x, y, z, cyzzz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += cyzzz[i]*U[ip[i]];
    actual = f4yzzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dyzzz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dyzzz */
  printf("f4 dyzzz max error=%g \n", maxerr);

  maxerr = 0.0;
  for(point=0; point<npoint; point++) /* f4 dzzzz */
  {
    m = nu3dzzzz(order, npoint, point, x, y, z, czzzz, ip);
    est = 0.0;
    for(i=0; i<m; i++) est += czzzz[i]*U[ip[i]];
    actual = f4zzzz(x[point],y[point],z[point]);
    error = actual-est;
    maxerr = max(maxerr,abs(error));
    if(debug) printf("f4 dzzzz %2d, est=%g, actual=%g, error=%g \n",
                     point, est, actual, error);
  } /* end f4 dzzzz */
  printf("f4 dzzzz max error=%g \n", maxerr);

  printf("end test_nuderiv3dg.c\n");
  return 0;
} /* end  test_nuderiv3dg.c */