// poly43.java  reads 35 coefficents (many may be zero) of polynomial
// in x, y, z, up to 4th power terms , all values and arithmetic is complex,
// the roots of the polynomial=0 are printed values of x, y, z
// real numbers are standard floating point, complex are (readl,imag)
// the possible coefficents c___ are:
// c + cx*x + cy*y + cz*z +                                         p1
// cxx*x^x  + cxy*x*y + cxz*x*z + cyy*y*y + cyz*y*z + czz*z*z +     p2
// cxxx*x*x*x + cxxy*x*x*y + cxxz*x*x*z + cxyy*x*y*y + cxyz*x*y*z + p3
// cxzz*x*z*z + cyyy*y*y*y + cyyz*y*y*z + cyzz*y*z*z + czzz*z*z*z + p3
// cxxxx*x*x*x*x + cxxxy*x*x*x*y + cxxxz*x*x*x*z + cxxyy*x*x*y*y +  p4
// cxxyz*x*x*y*z + cxxzz*x*x*z*z + cxyyy*x*y*y*y + cxyyz*x*y*y*z +  p4
// cxyzz*x*y*z*z + cxzzz*x*z*z*z + cyyyy*y*y*y*y + cyyyz*y*y*y*z +  p4
// cyyzz*y*y*z*z + cyzzz*y*z*z*z + czzzz*z*z*z*z                    p4
//
// these coefficents are stored in an array with subscrits P: term
// P[0] + P[1]*x + P[2]*y + P[3]*z + 
// P[4]*x^x  + P[5]*x*y + P[6]*x*z + P[7]*y*y + P[8]*y*z + P[9]*z*z +
// P[10]*x*x*x + P[11]*x*x*y + P[12]*x*x*z + P[13]*x*y*y + P[14]*x*y*z +
// P[15]*x*z*z + P[16]*y*y*y + P[17]*y*y*z + P[18]*y*z*z + P[19]*z*z*z +
// P[20]*x*x*x*x + P[21]*x*x*x*y + P[22]*x*x*x*z + P[23]*x*x*y*y +
// P[24]*x*x*y*z + P[25]*x*x*z*z + P[26]*x*y*y*y + P[27]*x*y*y*z +
// P[28]*x*y*z*z + P[29]*x*z*z*z + P[30]*y*y*y*y + P[31]*y*y*y*z +
// P[32]*y*y*z*z + P[33]*y*z*z*z + P[34]*z*z*z*z
//
// polynomial coefficients are stored in complex array  p  global
// precomputed values xp2=x*x xp3=x*x*x xp4=x*x*x*x  yp2,yp3,yp4  zp2,zp3,zp4
// err = eval(x, y, z) computes value of P at point x,y,z
// fpxv = fpx(x, y, z) computes derivative of P wrt x 
// fpyv = fpy(x, y, z) computes derivative of P wrt y 
// fpzv = fpz(x, y, z) computes derivative of P wrt z
//                     for newtons method
// sfpx(P, Der) computes 35 coefficients in Der that is derivative of P wrt x
// sfpy(P, Der) computes 35 coefficients in Der that is derivative of P wrt y
// sfpz(P, Der) computes 35 coefficients in Der that is derivative of P wrt z

import java.io.*;

public class poly43
{
  Complex P[]   = new Complex[35];
  Complex P1[]  = new Complex[35];
  Complex P2[]  = new Complex[35];
  Complex P3[]  = new Complex[35];
  Complex P4[]  = new Complex[35];
  Complex pwr[] = new Complex[35];
  Complex Der[] = new Complex[35];
  Complex Int[] = new Complex[35];
  Complex R[]   = new Complex[35];
  Complex C0 = new Complex(0.0,0.0);
  Complex C1 = new Complex(1.0,0.0);
  Complex C2 = new Complex(2.0,0.0);
  Complex C3 = new Complex(3.0,0.0);
  Complex C4 = new Complex(4.0,0.0);
  Complex Chalf = new Complex(0.5,0.0);
  Complex Cthird = new Complex(0.333333333333333,0.0);
  Complex Cquarter = new Complex(0.25,0.0);
  Complex xp2, xp3, xp4, yp2, yp3, yp4, zp2, zp3, zp4; // precomputed powers
  String eqn[] =    {"  c   [0]  ",  
                     "  cx  [1]  ",   
                     "  cy  [2]  ",
                     "  cz  [3]  ",
                     "  cxx [4]  ",
                     "  cxy [5]  ",
                     "  cxz [6]  ",
                     "  cyy [7]  ",
                     "  cyz [8]  ",
                     "  czz [9]  ",
                     "  cxxx [10]  ", 
                     "  cxxy [11]  ",
                     "  cxxz [12]  ",
                     "  cxyy [13]  ",
                     "  cxyz [14]  ",
                     "  cxzz [15]  ",
                     "  cyyy [16]  ",
                     "  cyyz [17]  ", 
                     "  cyzz [18]  ",
                     "  cyzzz [19]  ",   
                     "  cxxxx [20]  ", 
                     "  cxxxy [21]  ",  
                     "  cxxxz [22]  ",
                     "  cxxyy [23]  ",
                     "  cxxyz [24]  ",
                     "  cxxzz [25]  ",
                     "  cxyyy [26]  ",
                     "  cxyyz [27]  ",
                     "  cxyzz [28]  ",
                     "  cxzzz [29]  ",
                     "  cxyyy [30]  ",
                     "  cxyyz [31]  ",
                     "  cyyzz [32]  ",
                     "  cxxxy [33]  ",
                     "  cxxxx [34]  "};
  String ends = "end.................  ";

    
  public poly43() // constructor
  {
    Complex x, y, z, tx, ty, tz, err, val, tval;
    Complex bestx, besty, bestz, bestval, tbestval;
    Complex bestdx, bestdy, bestdz;
    Complex dx, dy, dz, tdx, tdy, tdz, tbdx, tbdy, tbdz;
    Complex mone = new Complex(-1.0,0.0);
    Complex eye = new Complex(0.0,1.0);
    Complex meye = new Complex(0.0,-1.0);
    double tol = 1.0e-6;
    Complex fpxv, fpyv, fpzv; // values of derivatives
    
    System.out.println("poly43 running");
    System.out.println("read data that defines polynomial");
    read_data();

    // initial guess
    x = C0; // one finds 0,0,0
    y = C0;
    z = C0;
    
    xp2 = x.multiply(x);
    xp3 = x.multiply(xp2);
    xp4 = x.multiply(xp3);
    yp2 = y.multiply(y);
    yp3 = y.multiply(yp2);
    yp4 = y.multiply(yp3);
    zp2 = z.multiply(z);
    zp3 = z.multiply(zp2);
    zp4 = z.multiply(zp3);

    // brute force compute roots of P
    val = eval(x, y, z);
    System.out.println("val="+val);
    bestval = val;
    bestx = x;
    besty = y;
    bestz = z;
    dx = C1; // factor * zero, one, none, eye, -eye
    dy = C1;
    dz = C1;
    tdx = C0;
    tdy = C0;
    tdz = C0;
    bestdx = C0;
    bestdy = C0;
    bestdz = C0;
    tx = C0;
    ty = C0;
    tz = C0;
    System.out.println("x="+x+", y="+y+", z="+z);	  
    System.out.println("val="+val);

    tbestval = bestval;
    tbdx = tdx;
    tbdy = tdy;
    tbdz = tdz;

    outerloop:
    for(int n=0; n<2; n++) // iterations
    {
      if(bestval.abs() < tol) break outerloop;
      for(int i=0; i<5; i++)
      {
        if(i==0) tdx = C0;                // just x
        if(i==1) tdx = dx;                // real
        if(i==2) tdx = dx.multiply(mone); // -real
        if(i==3) tdx = dx.multiply(eye);  // imag
        if(i==4) tdx = dx.multiply(meye); // -imag	     
        tx = x.add(tdx);
        for(int j=0; j<5; j++)
        {
	  if(j==0) tdy = C0;                // just y
          if(j==1) tdy = dy;                // real
          if(j==2) tdy = dy.multiply(mone); // -real
          if(j==3) tdy = dy.multiply(eye);  // imag
          if(j==4) tdy = dy.multiply(meye); // -imag	     
          ty = y.add(tdy);
	  for(int k=0; k<5; k++)
	  {
            if(k==0) tdz = C0;                // just z
            if(k==1) tdz = dz;                // real
            if(k==2) tdz = dz.multiply(mone); // -real
            if(k==3) tdz = dz.multiply(eye);  // imag
            if(k==4) tdz = dz.multiply(meye); // -imag	     
            tz = z.add(tdz);
            if(i!=0 | j!=0 | k!=0) // not start
	    {
	      val = eval(tx, ty, tz);
              System.out.println("tx="+tx+", ty="+ty+", tz="+tz);	  
              System.out.println("val="+val);
	      // check if best
	      if(val.abs()<tbestval.abs())
	      {
	        System.out.println("better");
                tbdx = tdx;
	        tbdy = tdy;
	        tbdz = tdz;
	        tbestval = val;
	      }
	      if(val.abs()<tol)
	      {
	        x = x.add(tbdx);
	        y = y.add(tbdy);
	        z = z.add(tbdz);
	        val = eval(x, y, z);
                System.out.println("x="+x+", y="+y+", z="+z);	  
                System.out.println("val="+val);
	        System.out.println("have solution");
	        break outerloop;
	      }
	    } // end if
	  } // end k
        } // end j
      } // end i
      x = x.add(tbdx);
      y = y.add(tbdy);
      z = z.add(tbdz);
      System.out.println("x="+x+", y="+y+", z="+z);	  
      System.out.println("tbdx="+tbdx+", tbdy="+tbdy+", tbdz="+tbdz);	  
      System.out.println("tbestval="+tbestval);
    } // n iterations  
    
    System.out.println(" ");
    System.out.println("Try Newton's method using derivative to find roots ");
    // initial guess
    x = new Complex(1.0, 0.0);
    y = new Complex(1.0, 0.0);
    z = new Complex(1.0, 0.0);
    
    xp2 = x.multiply(x);
    xp3 = x.multiply(xp2);
    xp4 = x.multiply(xp3);
    yp2 = y.multiply(y);
    yp3 = y.multiply(yp2);
    yp4 = y.multiply(yp3);
    zp2 = z.multiply(z);
    zp3 = z.multiply(zp2);
    zp4 = z.multiply(zp3);

    poly_roots(x,y,z,R); // first root of P into R start at x,y,x 

    poly_read("r1234.inp", P);
    poly_print("r1234.inp", P);
    x = new Complex(1.0, 0.0);
    y = new Complex(1.0, 0.0);
    z = new Complex(1.0, 0.0);
    poly_roots(x, y, z, R);
	
    System.out.println(" ");
    System.out.println("coming soon, deriv as new poly");
    sfpx(P, Der);
    poly_print("Deriv wrt x", Der);
    sfpy(P, Der);
    poly_print("Deriv wrt y", Der);
    sfpz(P, Der);
    poly_print("Deriv wrt z", Der);
    
    
    System.out.println("poly43 finisked");

  } // end poly43

  void poly_read(String fname, Complex Pa[]) // reads file into Pa
  {
    int n_input = 0;
    Complex a;
    
    try
    {
      FileReader fin = new FileReader(fname);
      BufferedReader in = new BufferedReader(fin);
      String input_line;
      int len;
      int index;
      int last;

      System.out.println("read_data running");

      for(int i=0; i<35; i++)
      {
        input_line = in.readLine(); // 35 p values
        System.out.println("input: "+input_line);
        index = 0;
        last = input_line.indexOf(' ',index+1);
	// System.out.println("index="+index+", last="+last);
        Pa[i] = Complex.parseComplex(input_line.substring(index,last));
        System.out.println("Pa["+i+"]="+Pa[i]);
      }
    }
    catch(IOException exception)
    {
      System.out.println(exception);
    }
  } // end read_poly
    
  void read_data() // read from command line into P
  {
    int n_input = 0;
    Complex a;
    
    try
    {
      BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
      String input_line;
      int len;
      int index;
      int last;

      System.out.println("read_data running");

      for(int i=0; i<35; i++)
      {
        input_line = in.readLine(); // 35 p values
        System.out.println("input: "+input_line);
        index = 0;
        last = input_line.indexOf(' ',index+1);
	// System.out.println("index="+index+", last="+last);
        P[i] = Complex.parseComplex(input_line.substring(index,last));
        System.out.println("P["+i+"]="+P[i]);
      }
    }
    catch(IOException exception)
    {
      System.out.println(exception);
    }
  } // end read data

  void make_pwr(Complex x, Complex y, Complex z)
  {
    pwr[0] = new Complex(1.0,0.0); // c
    pwr[1] = x;                    // x 
    pwr[2] = y;                    // y 
    pwr[3] = z;                    // z
    pwr[4] = x.multiply(x);        // xx
    pwr[5] = x.multiply(y);        // xy
    pwr[6] = x.multiply(z);        // xz
    pwr[7] = y.multiply(y);        // yy
    pwr[8] = y.multiply(z);        // yz
    pwr[9] = z.multiply(z);        // zz
    pwr[10] = pwr[4].multiply(x);  // xxx
    pwr[11] = pwr[4].multiply(y);  // xxy
    pwr[12] = pwr[4].multiply(z);  // xxz
    pwr[13] = pwr[5].multiply(y);  // xyy
    pwr[14] = pwr[5].multiply(z);  // xyz
    pwr[15] = pwr[6].multiply(z);  // xzz
    pwr[16] = pwr[7].multiply(y);  // yyy
    pwr[17] = pwr[7].multiply(z);  // yyz
    pwr[18] = pwr[8].multiply(z);  // yzz
    pwr[19] = pwr[9].multiply(z);  // zzz
    pwr[20] = pwr[10].multiply(x); // xxxx
    pwr[21] = pwr[10].multiply(y); // xxxy
    pwr[22] = pwr[10].multiply(z); // xxxz
    pwr[23] = pwr[4].multiply(pwr[7]); // xxyy
    pwr[24] = pwr[4].multiply(pwr[8]); // xxyz
    pwr[25] = pwr[5].multiply(pwr[9]); // xxzz
    pwr[26] = x.multiply(pwr[16]);     // xyyy
    pwr[27] = pwr[5].multiply(pwr[8]); // xyyz
    pwr[28] = pwr[5].multiply(pwr[9]); // xyzz
    pwr[29] = pwr[6].multiply(pwr[9]); // xzzz
    pwr[30] = pwr[16].multiply(y);     // yyyy
    pwr[31] = pwr[16].multiply(z);     // yyyz
    pwr[32] = pwr[7].multiply(pwr[9]); // yyzz
    pwr[33] = y.multiply(pwr[9]);      // yzzz
    pwr[34] = pwr[19].multiply(z);     // zzzz
  } // end make pwr

  Complex eval(Complex x, Complex y, Complex z) // uses P and pwr
  {
    Complex sum = new Complex(0.0,0.0);
    // build pwr
    make_pwr(x, y, z);
    for(int i=0; i<35; i++)
    {
      sum = sum.add(P[i].multiply(pwr[i]));
    }
    return sum;
  } // end eval

  void poly_zero(Complex Pa[])
  {
    for(int i=0; i<35; i++) Pa[i] = C0;
  } // end poly_zero
  
  void poly_add(Complex Pa[], Complex Pb[], Complex Pc[]) // Pc = Pa + Pb
  {
    for(int i=0; i<35; i++) Pc[i] = Pa[i].add(Pb[i]);
  } // emd poly_add
  
  void poly_sub(Complex Pa[], Complex Pb[], Complex Pc[]) // Pc = Pa - Pb
  {
    for(int i=0; i<35; i++) Pc[i] = Pa[i].subtract(Pb[i]);
  } // emd poly_sub
  
  void poly_mul(Complex Pa[], Complex Pb[], Complex Pp[]) // Pp = Pa * Pb
  {
    int k;
    poly_zero(Pp);
    for(int i=0; i<20; i++) // [0] *c  [1] *x  [2] *y  [3] * z [4] *xx
    {
      for(int j=0; j<35; j++) // [0] *c  [1] *x  [2] *y  [3] * z [4] *xx
      {
	k = term(i, j); // get i*j term = k
	if(k<0)
	{
	  // check for invalid overflow  ???
	  break; // not a valid term 0..34
	} // end if
	// System.out.println("Pp["+k+"]+=Pa["+i+"]*Pb["+j+"]");
	// System.out.println(Pp[k]+"+="+Pa[i]+"*"+Pb[j]);
        Pp[k] = Pp[k].add(Pa[i].multiply(Pb[j]));
      } // end j
    } // end i
  } // emd poly_mul
  
  void poly_div(Complex Pa[], Complex Pb[], Complex Pq[], Complex Pr[])
  {	                                    // Pq = Pa / Pb  quotienrt
                                            // pr = Pa % Pb  remainder

  } // emd poly_div



  void poly_print(String name, Complex Pa[])
  {
    System.out.println(" ");
    System.out.println("poly "+name+" is: ");
    int n = 35;
    for(int i=34; i>0; i--)
    {
      if(Pa[i].abs()==0)
      {
        n = n-1;
      }
      else
      {
        break;
      }
    } // end i
    for(int i=0; i<n; i++)
    {
      System.out.println(Pa[i]+eqn[i]);
    } // end i
    System.out.println("above ["+n+"] are zero ");
    System.out.println(" ");
  } // end poly_print

  int term(int i, int j)
  {
    int tt[][] =
    {{ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,16,17,18,19,20,
      21,22,23,24,25,26,27,28,29,30,31,23,33,34}, // *C  [0]
     { 1, 4, 5, 6,10,11,12,13,14,15,20,21,22,23,24,25,26,27,28,29,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *x  [1]
     { 2, 5, 7, 8,11,13,14,16,17,18,21,23,24,26,27,28,30,31,32,33,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *y  [2]
     { 3, 6, 8, 9,12,14,15,17,18,19,22,24,25,27,28,29,31,32,33,34,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *z  [3]
     { 4,10,11,12,20,21,22,23,24,25,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *xx [4]
     { 5,11,13,14,21,23,24,26,27,28,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *xy [5]
     { 6,12,14,15,22,24,25,27,28,29,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *xz [6]
     { 7,13,16,17,23,26,27,30,31,32,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *yy [7]
     { 8,14,17,18,24,27,28,31,32,33,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *yz [8]
     { 9,15,18,19,25,28,29,32,33,34,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *zz [9]
     {10,20,21,22,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *xxx [10]
     {11,21,23,24,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *xxy [11]
     {12,22,24,25,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *xxz [12]
     {13,23,26,27,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *xyy [13]
     {14,24,27,28,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *xyz [14]
     {15,25,28,29,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *xzz [15]
     {16,26,30,31,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *yyy [16]
     {17,27,31,32,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *yyz [17]
     {18,28,32,33,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}, // *yzz [18]
     {19,29,33,34,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}};// *zzz [19]
     // no i in 20..34
    if(i>19) return -1;
    return tt[i][j];
  }

  Complex pr2(Complex c1, Complex c2) // product 2
  {
      return c1.multiply(c2);
  }

  Complex pr3(Complex c1, Complex c2, Complex c3) // product 3
  {
      return c1.multiply(c2.multiply(c3));
  }

  Complex pr4(Complex c1, Complex c2, Complex c3, Complex c4) // product 4
  {
      return c1.multiply(c2.multiply(c3.multiply(c4)));
  }

  Complex pr5(Complex c1, Complex c2, Complex c3, Complex c4, Complex c5)
                                                              // product 5
  {
      return c1.multiply(c2.multiply(c3.multiply(c4.multiply(c5))));
  }

  //  fpx is derivative of polynomial P with respect to  x
  Complex fpx(Complex x, Complex y, Complex z)
  {                                         // term  of derivative
                                            // 0     f' P[0]c  wrt x
    Complex deriv = P[1];                   // P[1]  f' P[1]x
                                            // 0     f' P[2]y
                                            // 0     f' P[3]z
    deriv=deriv.add(pr3(P[4],C2,x));        // 2x    f' P[4]xx
    deriv=deriv.add(pr2(P[5],y));           // y     f' P[5]xy
    deriv=deriv.add(pr2(P[6],z));           // z     f' P[6]xz
                                            // 0     f' P[7]yy
                                            // 0     f' P[8]yz
                                            // 0     f' P[9]zz
    deriv=deriv.add(pr3(P[10],C3,xp2));     // 3xx   f' P[10]xxx
    deriv=deriv.add(pr4(P[11],C2,x,y));     // 2xy   f' P[11]xxy
    deriv=deriv.add(pr4(P[12],C2,x,z));     // 2xz   f' P[12]xxz
    deriv=deriv.add(pr2(P[13],yp2));        // yy    f' P[13]xyy
    deriv=deriv.add(pr3(P[14],y,z));        // yz    f' P[14]xyz
    deriv=deriv.add(pr2(P[15],zp2));        // zz    f' P[15]xzz
                                            // 0     f' P[16]yyy
                                            // 0     f' P[17]yyz
                                            // 0     f' P[18]yzz
                                            // 0     f' P[19]zzz
    deriv=deriv.add(pr3(P[20],C4,xp3));     // 4xxx  f' P[20]xxxx
    deriv=deriv.add(pr4(P[21],C3,xp2,y));   // 3xxy  f' P[21]xxxy
    deriv=deriv.add(pr4(P[22],C3,xp2,z));   // 3xxz  f' P[22]xxxz
    deriv=deriv.add(pr4(P[23],C2,x,yp2));   // 2xyy  f' P[23]xxyy
    deriv=deriv.add(pr5(P[24],C2,x,y,z));   // 2xyz  f' P[24]xxyz
    deriv=deriv.add(pr4(P[25],C2,x,zp2));   // 2xzz  f' P[25]xxzz
    deriv=deriv.add(pr2(P[26],yp3));        // yyy   f' P[26]xyyy
    deriv=deriv.add(pr3(P[27],yp2,z));      // yyz   f' P[27]xyyz
    deriv=deriv.add(pr3(P[28],y,zp2));      // yzz   f' P[28]xyzz
    deriv=deriv.add(pr2(P[29],zp3));        // zzz   f' P[29]xzzz
                                            // 0     f' P[30]yyyy
                                            // 0     f' P[31]yyyz
                                            // 0     f' P[32]yyzz
                                            // 0     f' P[33]yzzz
                                            // 0     f' P[34]zzzz
    return deriv;
  } // end fpx

  //  fpy is derivative of polynomial P with respect to  y
  Complex fpy(Complex x, Complex y, Complex z)
  {
                                            // 0     f' P[0]c  wrt y
                                            // 0     f' P[1]x
    Complex deriv = P[2];                   // c     f' P[2]y
                                            // 0     f' P[3]z
                                            // 0     f' P[4]xx
    deriv=deriv.add(pr2(P[5],x));           // x     f' P[5]xy
                                            // 0     f' P[6]xz
    deriv=deriv.add(pr3(P[7],C2,y));        // 2y    f' P[7]yy
    deriv=deriv.add(pr2(P[8],z));           // z     f' P[8]yz
                                            // 0     f' P[9]zz
                                            // 0     f' P[10]xxx
    deriv=deriv.add(pr2(P[11],xp2));        // xx    f' P[11]xxy
                                            // 0     f' P[12]xxz
    deriv=deriv.add(pr4(P[13],x,C2,y));     // x2y   f' P[13]xyy
    deriv=deriv.add(pr3(P[14],x,z));        // xz    f' P[14]xyz
                                            // 0     f' P[15]xzz
    deriv=deriv.add(pr3(P[16],C3,yp2));     // 3yy   f' P[16]yyy
    deriv=deriv.add(pr4(P[17],C2,y,z));     // 2yz   f' P[17]yyz
    deriv=deriv.add(pr2(P[18],zp2));        // zz    f' P[18]yzz
                                            // 0     f' P[19]zzz
                                            // 0     f' P[20]xxxx
    deriv=deriv.add(pr2(P[21],xp3));        // xxx   f' P[21]xxxy
                                            // 3xxz  f' P[22]xxxz
    deriv=deriv.add(pr4(P[23],xp2,C2,y));   // xx2y  f' P[23]xxyy
    deriv=deriv.add(pr3(P[24],xp2,z));      // xxz   f' P[24]xxyz
                                            // 0     f' P[25]xxzz
    deriv=deriv.add(pr4(P[26],x,C3,yp2));   // x3yy  f' P[26]xyyy
    deriv=deriv.add(pr5(P[27],x,C2,y,z));   // x2yz  f' P[27]xyyz
    deriv=deriv.add(pr3(P[28],x,zp2));      // xzz   f' P[28]xyzz
                                            // 0     f' P[29]xzzz
    deriv=deriv.add(pr3(P[30],C4,yp3));     // 4yyy  f' P[30]yyyy   
    deriv=deriv.add(pr4(P[31],C3,yp2,z));   // 3yyz  f' P[31]yyyz
    deriv=deriv.add(pr4(P[32],C2,y,zp2));   // 2yzz  f' P[32]yyzz
    deriv=deriv.add(pr2(P[33],zp3));        // zzz   f' P[33]yzzz
                                            // 0     f' P[34]zzzz
    return deriv;
  } // end fpy
    
  //  fpz is derivative of polynomial P with respect to  z
  Complex fpz(Complex x, Complex y, Complex z)
  {                                          // term
                                             // 0     f' P[0]c  wrt z
                                             // 0     f' P[1]x
                                             // 0     f' P[2]y
    Complex deriv = P[3];                    // c     f' P[3]z
                                             // 0     f' P[4]xx
                                             // 0     f' P[5]xy
    deriv=deriv.add(pr2(P[6],x));            // x     f' P[6]xz
                                             // 0     f' P[7]yy
    deriv=deriv.add(pr2(P[8],z));            // y     f' P[8]yz
    deriv=deriv.add(pr3(P[9],C2,z));         // 2z    f' P[9]zz
                                             // 0     f' P[10]xxx
                                             // 0     f' P[11]xxy
    deriv=deriv.add(pr2(P[12],xp2));         // xx    f' P[12]xxz
                                             // 0     f' P[13]xyy
    deriv=deriv.add(pr3(P[14],x,y));         // xyz   f' P[14]xyz
    deriv=deriv.add(pr4(P[15],x,C2,z));      // x2z   f' P[15]xzz
                                             // 0     f' P[16]yyy
    deriv=deriv.add(pr2(P[17],yp2));         // yy    f' P[17]yyz
    deriv=deriv.add(pr4(P[18],y,C2,z));      // y2z   f' P[18]yzz
    deriv=deriv.add(pr3(P[19],C3,zp2));      // 3zz   f' P[19]zzz
                                             // 0     f' P[20]xxxx
                                             // 0     f' P[21]xxxy
    deriv=deriv.add(pr2(P[22],xp3));         // xxx   f' P[22]xxxz
                                             // 0     f' P[23]xxyy
    deriv=deriv.add(pr3(P[24],xp2,y));       // xxy   f' P[24]xxyz
    deriv=deriv.add(pr4(P[25],xp2,C2,z));    // xx2z  f' P[25]xxzz
                                             // 0     f' P[26]xyyy
    deriv=deriv.add(pr3(P[27],x,yp2));       // xyy   f' P[27]xyyz
    deriv=deriv.add(pr5(P[28],x,y,C2,z));    // xy2z  f' P[28]xyzz
    deriv=deriv.add(pr4(P[29],x,C3,z));      // x3zz  f' P[29]xzzz
                                             // 0     f' P[30]yyyy   
    deriv=deriv.add(pr2(P[31],yp3));         // yyy   f' P[31]yyyz
    deriv=deriv.add(pr4(P[32],yp2,C2,z));    // yy2z  f' P[32]yyzz
    deriv=deriv.add(pr4(P[33],y,C3,z));      // y3zz  f' P[33]yzzz
    deriv=deriv.add(pr3(P[34],C4,z));        // 4zzz  f' P[34]zzzz
    
    return deriv;
  } // end fpz

  // sfpx(P, Der) computes 35 coefficients in der that is derivative of P wrt x
  void sfpx(Complex P[], Complex Der[]) // ?? at y,z  needed or P[2],P[3]
  {
    for(int i=0; i<35; i++) Der[i] = C0;    // for none
    Der[0] = P[1];                          // [0]c = f'x
      	                                    // none f'y
      	                                    // none f'y
      	                                    // none f'z
    Der[1] = C2;                            // [1]x = f'xx ??
    Der[1] = Der[1].add(P[2]); // add(y)    // [1]x = f'xy ??  
  } // end sfpx
    
  // sfpy(P, Der) computes 35 coefficients in der that is derivative of P wrt y
  void sfpy(Complex P[], Complex Der[]) // ?? at y,z  needed or P[2],P[3]
  {
    for(int i=0; i<35; i++) Der[i] = C0;    // for none
    Der[0] = P[1];                          // [0]c = f'x
      	                                    // none f'y
      	                                    // none f'y
      	                                    // none f'z
    Der[2] = C2;                            // [1]x = f'yy ??
    Der[2] = Der[1].add(P[2]); // add(y)    // [1]x = f'xy ??  
  } // end sfpy

  // sfpz(P, Der) computes 35 coefficients in der that is derivative of P wrt z
  void sfpz(Complex P[], Complex Der[]) // ?? at x,y  needed or P[2],P[3]
  {
    for(int i=0; i<35; i++) Der[i] = C0;    // for none
    Der[0] = P[1];                          // [0]c = f'x
      	                                    // none f'y
      	                                    // none f'y
      	                                    // none f'z
    Der[3] = C2;                            // [1]x = f'xx ??
    Der[3] = Der[1].add(P[2]); // add(y)    // [1]x = f'xy ??  
  } // end sfpz


  void poly_roots(Complex x, Complex y, Complex z, Complex R[]) // of P
  {    
    Complex tx, ty, tz, err, val, tval;
    Complex bestx, besty, bestz, bestval, tbestval;
    Complex fpxv, fpyv, fpzv;
    double tol = 0.1;
    
    System.out.println("poly_roots start x="+x+", y="+y+", z="+z);	  
    val = eval(x, y, z);
    System.out.println("at initial x,y,z val="+val);
    bestval = val;
    bestx = x;
    besty = y;
    bestz = z;
    tx = x;
    ty = y;
    tz = z;

    fpxv = fpx(bestx, besty, bestz);
    System.out.println("fpxv="+fpxv);
    fpyv = fpy(bestx, besty, bestz);
    System.out.println("fpyv="+fpyv);
    fpzv = fpz(bestx, besty, bestz);
    System.out.println("fpzv="+fpzv);

    if(fpxv.abs()>tol) tx = x.subtract(val.divide(fpxv));
    if(fpyv.abs()>tol) ty = y.subtract(val.divide(fpyv));
    if(fpzv.abs()>tol) tz = z.subtract(val.divide(fpzv));
    tval = eval(tx, ty, tz);
    System.out.println("tx="+tx);
    System.out.println("ty="+ty);
    System.out.println("tz="+tz);	  
    System.out.println("tval="+tval);
    if(tval.abs()< bestval.abs())
    {
      bestval=tval;
      bestx = tx;
      besty = ty;
      bestz = tz;
    }
    else
    {
      System.out.println("did not converge");
      R[0] = bestx;
      R[1] = besty;
      R[2] = bestz;
      return;
    }

    for(int k=0; k<10; k++)
    {	    
      fpxv = fpx(bestx, besty, bestz);
      System.out.println("fpxv="+fpxv);
      fpyv = fpy(bestx, besty, bestz);
      System.out.println("fpyv="+fpyv);
      fpzv = fpz(bestx, besty, bestz);
      System.out.println("fpzv="+fpzv);

      if(fpxv.abs()>tol) tx = bestx.subtract(tval.divide(fpxv));
      if(fpyv.abs()>tol) ty = besty.subtract(tval.divide(fpyv));
      if(fpzv.abs()>tol) tz = bestz.subtract(tval.divide(fpzv));
      tval = eval(tx, ty, tz);
      System.out.println("tx="+tx);
      System.out.println("ty="+ty);
      System.out.println("tz="+tz);	  
      System.out.println("tval="+tval);
      if(tval.abs()< bestval.abs())
      {
        bestval=tval;
        bestx = tx;
        besty = ty;
        bestz = tz;
      }
      else
      {
        System.out.println("did not converge");
        R[0] = bestx;
        R[1] = besty;
        R[2] = bestz;
        return;
      }
    } // end k
    R[0] = bestx;
    R[1] = besty;
    R[2] = bestz;
  } // end poly_root


  public static void main (String args[])
  {
    new poly43(); // construct and execute
  }
} // end poly34.java