// test_nuderiv3dg.java   test  nuderiv3dg  with 36 points 
//          order 1  minimum independent points 4   Ux ... on first degree
//          order 2  minimum independent points 10  Uxx ... on second degree
//          order 3  minimum independent points 20  Uxxx ... on third degree
//          order 4  minimum independent points 35  Uxxxx ... on fourth degree
//          going beyond minimum points does not improve accuracy
//          May need more than 35 if not linearly independent
 
class test_nuderiv3dg
{
  test_nuderiv3dg()
  {
    int debug = 0;
    int m, order;
    int n = 36; // max 4th order this test 
    int ordpt[] = {0, 4, 10, 20, 35};
    double x[] = {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.92, 0.98};
    double y[] = {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.93, 0.99};
    double z[] = {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.65, 0.82};
    double U[] = new double[36];
    double est, actual, dx, dy;
    double cx[] = new double[36];
    double cy[] = new double[36];
    double cz[] = new double[36];
    double cxx[] = new double[36];
    double cxy[] = new double[36];
    double cxz[] = new double[36];
    double cyy[] = new double[36];
    double cyz[] = new double[36];
    double czz[] = new double[36];
    double cxxx[] = new double[36];
    double cxxy[] = new double[36];
    double cxxz[] = new double[36];
    double cxyy[] = new double[36];
    double cxyz[] = new double[36];
    double cxzz[] = new double[36];
    double cyyy[] = new double[36];
    double cyyz[] = new double[36];
    double cyzz[] = new double[36];
    double czzz[] = new double[36];
    double cxxxx[] = new double[36];
    double cxxxy[] = new double[36];
    double cxxxz[] = new double[36];
    double cxxyy[] = new double[36];
    double cxxyz[] = new double[36];
    double cxxzz[] = new double[36];
    double cxyyy[] = new double[36];
    double cxyyz[] = new double[36];
    double cxyzz[] = new double[36];
    double cxzzz[] = new double[36];
    double cyyyy[] = new double[36];
    double cyyyz[] = new double[36];
    double cyyzz[] = new double[36];
    double cyzzz[] = new double[36];
    double czzzz[] = new double[36];
    int ip[] = new int[36];
    double error, maxerror;

    nuderiv3dg N = new nuderiv3dg();

    System.out.println("test_nuderiv3dg.java running");

    // f3(x,y,z) is 3th order
    order = 3;
    n = ordpt[order]+1; // give extra point to prevent singular
    for(int i=0; i<n; i++) U[i] = f3(x[i], y[i], z[i]);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Ux 
    {
      m = N.nu3dx(order, n, point, x, y, z, cx, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cx[i]*U[ip[i]];
      actual = f3x(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Ux 
    System.out.println("Ux order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uy 
    {
      m = N.nu3dy(order, n, point, x, y, z, cy, ip);
      // System.out.println("order="+order+", n="+n+", point="+point);
      est = 0.0;
      for(int i=0; i<m; i++) est += cy[i]*U[ip[i]];
      actual = f3y(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
      // System.out.println("cy  est="+est+", actual="+actual+
      //                   ", error="+ error);
    } // end Uy 
    System.out.println("Uy order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uz 
    {
      m = N.nu3dz(order, n, point, x, y, z, cz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cz[i]*U[ip[i]];
      actual = f3z(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uz
    System.out.println("Uz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxx 
    {
      m = N.nu3dxx(order, n, point, x, y, z, cxx, ip);
       est = 0.0;
      for(int i=0; i<m; i++) est += cxx[i]*U[ip[i]];
      actual = f3xx(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxx 
    System.out.println("Uxx order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxy 
    {
      m = N.nu3dxy(order, n, point, x, y, z, cxy, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxy[i]*U[ip[i]];
      actual = f3xy(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxy 
    System.out.println("Uxy order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxz 
    {
      m = N.nu3dxz(order, n, point, x, y, z, cxz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxz[i]*U[ip[i]];
      actual = f3xz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxz 
    System.out.println("Uxz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uyy 
    {
      m = N.nu3dyy(order, n, point, x, y, z, cyy, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cyy[i]*U[ip[i]];
      actual = f3yy(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uyy 
    System.out.println("Uyy order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uyz 
    {
      m = N.nu3dyz(order, n, point, x, y, z, cyz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cyz[i]*U[ip[i]];
      actual = f3yz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uyz
    System.out.println("Uyz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uzz 
    {
      m = N.nu3dzz(order, n, point, x, y, z, czz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += czz[i]*U[ip[i]];
      actual = f3zz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uzz
    System.out.println("Uzz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxxx
    {
      m = N.nu3dxxx(order, n, point, x, y, z, cxxx, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxxx[i]*U[ip[i]];
      actual = f3xxx(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxxx 
    System.out.println("Uxxx order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxxy 
    {
      m = N.nu3dxxy(order, n, point, x, y, z, cxxy, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxxy[i]*U[ip[i]];
      actual = f3xxy(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxxy 
    System.out.println("Uxxy order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxxz 
    {
      m = N.nu3dxxz(order, n, point, x, y, z, cxxz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxxz[i]*U[ip[i]];
      actual = f3xxz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxxz 
    System.out.println("Uxxz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxyy 
    {
      m = N.nu3dxyy(order, n, point, x, y, z, cxyy, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxyy[i]*U[ip[i]];
      actual = f3xyy(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxyy 
    System.out.println("Uxyy order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxyz 
    {
      m = N.nu3dxyz(order, n, point, x, y, z, cxyz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxyz[i]*U[ip[i]];
      actual = f3xyz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxyz 
    System.out.println("Uxyz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uyyy 
    {
      m = N.nu3dyyy(order, n, point, x, y, z, cyyy, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cyyy[i]*U[ip[i]];
      actual = f3yyy(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uyyy 
    System.out.println("Uyyy order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uyyz 
    {
      m = N.nu3dyyz(order, n, point, x, y, z, cyyz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cyyz[i]*U[ip[i]];
      actual = f3yyz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uyyz 
    System.out.println("Uyyz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uyzz 
    {
      m = N.nu3dyzz(order, n, point, x, y, z, cyzz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cyzz[i]*U[ip[i]];
      actual = f3yzz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uyzz 
    System.out.println("Uyzz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uzzz 
    {
      m = N.nu3dzzz(order, n, point, x, y, z, czzz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += czzz[i]*U[ip[i]];
      actual = f3zzz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uzzz 
    System.out.println("Uzzz order="+order+", n="+n+", maxerr="+maxerror);

    // fourth order derivatives
    n = 36;
    order = 4;
    for(int i=0; i<n; i++) U[i] = f4(x[i], y[i], z[i]);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxxxx 
    {
      m = N.nu3dxxxx(order, n, point, x, y, z, cxxxx, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxxxx[i]*U[ip[i]];
      actual = f4xxxx(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxxxx 
    System.out.println("Uxxxx order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxxxy 
    {
      m = N.nu3dxxxy(order, n, point, x, y, z, cxxxy, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxxxy[i]*U[ip[i]];
      actual = f4xxxy(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxxxy 
    System.out.println("Uxxxy order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxxxz 
    {
      m = N.nu3dxxxz(order, n, point, x, y, z, cxxxz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxxxz[i]*U[ip[i]];
      actual = f4xxxz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxxxy 
    System.out.println("Uxxxz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxxyy 
    {
      m = N.nu3dxxyy(order, n, point, x, y, z, cxxyy, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxxyy[i]*U[ip[i]];
      actual = f4xxyy(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxxyy 
    System.out.println("Uxxyy order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxxyz 
    {
      m = N.nu3dxxyz(order, n, point, x, y, z, cxxyz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxxyz[i]*U[ip[i]];
      actual = f4xxyz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxxyz
    System.out.println("Uxxyz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxxzz 
    {
      m = N.nu3dxxzz(order, n, point, x, y, z, cxxzz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxxzz[i]*U[ip[i]];
      actual = f4xxzz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxxzz
    System.out.println("Uxxzz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxyyy 
    {
      m = N.nu3dxyyy(order, n, point, x, y, z, cxyyy, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxyyy[i]*U[ip[i]];
      actual = f4xyyy(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxyyy 
    System.out.println("Uxyyy order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxyyz 
    {
      m = N.nu3dxyyz(order, n, point, x, y, z, cxyyz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxyyz[i]*U[ip[i]];
      actual = f4xyyz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxyyz
    System.out.println("Uxyyz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxyzz 
    {
      m = N.nu3dxyzz(order, n, point, x, y, z, cxyzz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxyzz[i]*U[ip[i]];
      actual = f4xyzz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxyzz
    System.out.println("Uxyzz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uxzzz 
    {
      m = N.nu3dxzzz(order, n, point, x, y, z, cxzzz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cxzzz[i]*U[ip[i]];
      actual = f4xzzz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uxzzz
    System.out.println("Uxzzz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uyyyy 
    {
      m = N.nu3dyyyy(order, n, point, x, y, z, cyyyy, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cyyyy[i]*U[ip[i]];
      actual = f4yyyy(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uyyyy
    System.out.println("Uyyyy order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uyyyz
    {
      m = N.nu3dyyyz(order, n, point, x, y, z, cyyyz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cyyyz[i]*U[ip[i]];
      actual = f4yyyz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uyyyz
    System.out.println("Uyyyz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uyyzz
    {
      m = N.nu3dyyzz(order, n, point, x, y, z, cyyzz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cyyzz[i]*U[ip[i]];
      actual = f4yyzz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uyyzz
    System.out.println("Uyyzz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uyzzz
    {
      m = N.nu3dyzzz(order, n, point, x, y, z, cyzzz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += cyzzz[i]*U[ip[i]];
      actual = f4yzzz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uyzzz
    System.out.println("Uyzzz order="+order+", n="+n+", maxerr="+maxerror);

    maxerror = 0.0;
    for(int point=0; point<n; point++) // Uzzzz
    {
      m = N.nu3dzzzz(order, n, point, x, y, z, czzzz, ip);
      est = 0.0;
      for(int i=0; i<m; i++) est += czzzz[i]*U[ip[i]];
      actual = f4zzzz(x[point],y[point],z[point]);
      error = actual-est;
      maxerror = Math.max(maxerror,Math.abs(error));
    } // end Uzzzz
    System.out.println("Uzzzz order="+order+", n="+n+", maxerr="+maxerror);

    System.out.println("test_nuderiv3dg.java finished");
  }

  // test functions
  // 1st order test case 
  double f1(double x, double y, double z)
  {
    return x + 2.0*y + 3.0*z + 4.0;
  }
  double f1x(double x, double y, double z)
  {
    return 1.0;
  }
  double f1y(double x, double y, double z)
  {
    return 2.0;
  }
  double f1z(double x, double y, double z)
  {
    return 3.0;
  }

  // 2nd order test case 
  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;
  }
  double f2x(double x, double y, double z)
  {
    return 2.0*x + 2.0*y + 3.0*z + 7.0;
  }
  double f2y(double x, double y, double z)
  {
    return 2.0*x + 8.0*y + 5.0*z + 8.0;
  }
  double f2z(double x, double y, double z)
  {
   return 3.0*x + 5.0*y + 12.0*z + 9.0;
  }
  double f2xx(double x, double y, double z)
  {
    return 2.0;
  }
  double f2xy(double x, double y, double z)
  {
    return 2.0;
  }
  double f2xz(double x, double y, double z)
  {
    return 3.0;
  }
  double f2yy(double x, double y, double z)
  {
    return 8.0;
  }
  double f2yz(double x, double y, double z)
  {
    return 5.0;
  }
  double f2zz(double x, double y, double z)
  {
    return 12.0;
  }

  // 3rd order test case 
  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;
  }
  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;
  }
  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;
  }
  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;
  }
  double f3xx(double x, double y, double z)
  {
    return 6.0*x + 4.0*y +  6.0*z + 22.0;
  }
  double f3xy(double x, double y, double z)
  {
    return 4.0*x + 8.0*y + 5.0*z + 12.0  ;
  }
  double f3xz(double x, double y, double z)
  {
    return 6.0*x + 5.0*y + 12.0*z + 13.0;
  }
  double f3yy(double x, double y, double z)
  {
    return 8.0*x + 42.0*y + 16.0*z + 28.0;
  }
  double f3yz(double x, double y, double z)
  {
    return 5.0*x + 16.0*y + 18.0*z + 15.0;
  }
  double f3zz(double x, double y, double z)
  {
    return 12.0*x + 18.0*y + 60.0*z + 32.0;
  }
  double f3xxx(double x, double y, double z)
  {
    return 6.0;
  }
  double f3xxy(double x, double y, double z)
  {
    return 4.0;
  }
  double f3xxz(double x, double y, double z)
  {
    return 6.0;
  }
  double f3xyy(double x, double y, double z)
  {
    return 8.0;
  }
  double f3xyz(double x, double y, double z)
  {
    return 5.0;
  }
  double f3xzz(double x, double y, double z)
  {
    return 12.0;
  }
  double f3yyy(double x, double y, double z)
  {
    return 42.0;
  }
  double f3yyz(double x, double y, double z)
  {
    return 16.0;
  }
  double f3yzz(double x, double y, double z)
  {
    return 18.0;
  }
  double f3zzz(double x, double y, double z)
  {
    return 60.0;
  }

  // 4th order test case  
  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;
  }
  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;
  }
  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;
  }
  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;
  }
  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;
  }
  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;
  }
  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;
  }
  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;
  }
  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;
  }
  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;
  }
  double f4xxx(double x, double y, double z) // df4/dxxx 
  {
    return 24.0*x + 12.0*y + 18.0*z + 96.0;
  }
    double f4xxy(double x, double y, double z) // df4/dxxy 
  {
    return 12.0*x + 16.0*y + 10.0*z + 34.0;
  }
  double f4xxz(double x, double y, double z) // df4/dxxz 
  {
    return 18.0*x + 10.0*y + 24.0*z + 36.0;
  }
  double f4xyy(double x, double y, double z) // df4/dxyy 
  {
    return  16.0*x + 42.0*y + 16.0*z + 38.0;
  }
  double f4xyz(double x, double y, double z) // df4/dxyz 
  {
    return 10.0*x + 16.0*y +  18.0*z + 20.0;
  }
  double f4xzz(double x, double y, double z) // df4/dxzz 
  {
    return  24.0*x + 18.0*y +  60.0*z + 42.0;
  }
  double f4yyy(double x, double y, double z) // df4/dyyy 
  {
    return 42.0*x + 264.0*y + 72.0*z + 132.0;
  }
  double f4yyz(double x, double y, double z) // df4/dyyz 
  {
    return  16.0*x + 72.0*y + 52.0*z + 46.0;
  }
  double f4zzz(double x, double y, double z) // df4/dzzz 
  {
    return   60.0*x + 84.0*y + 360.0*z + 150.0;
  }
  double f4xxxx(double x, double y, double z) // df4/dxxxx 
  {
    return 24.0;
  }
  double f4xxxy(double x, double y, double z) // df4/dxxxy 
  {
    return 12.0;
  }
  double f4xxxz(double x, double y, double z) // df4/dxxxz 
  {
    return 18.0;
  }
  double f4xxyy(double x, double y, double z) // df4/dxxyy 
  {
    return 16.0;
  }
  double f4xxyz(double x, double y, double z) // df4/dxxyz 
  {
    return 10.0;
  }
  double f4xxzz(double x, double y, double z) // df4/dxxzz 
  {
    return  24.0;
  }
  double f4xyyy(double x, double y, double z) // df4/dxyyy 
  {
    return  42.0;
  }
  double f4xyyz(double x, double y, double z) // df4/dxyyz 
  {
    return  16.0;
  }
  double f4xyzz(double x, double y, double z) // df4/dxyzz 
  {
    return  18.0;
  }
  double f4xzzz(double x, double y, double z) // df4/dxzzz 
  {
    return  60.0;
  }
  double f4yyyy(double x, double y, double z)  // df4/dyyyy 
  {
    return 264.0;
  }
  double f4yyyz(double x, double y, double z)  // df4/dyyyz 
  {
    return  72.0;
  }
  double f4yyzz(double x, double y, double z)  // df4/dyyzz 
  {
    return  52.0;
  }
  double f4yzzz(double x, double y, double z)  // df4/dyzzz 
  {
    return  84.0;
  }
  double f4zzzz(double x, double y, double z) // df4/dzzzz 
  {
    return 360.0;
  }

  public static void main (String[] args)
  {
    new test_nuderiv3dg();
  }
} // end  class test_nuderiv3dg