// deriv_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 (read,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 deriv_poly43 { Complex P[] = new Complex[35]; // coeficients of polynomial Complex pwr[] = new Complex[35]; // value of term of polynomail Complex Der[] = new Complex[35]; // coeficients of derivative Complex C0 = new Complex(0.0,0.0); // Complex constants 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 mone = new Complex(-1.0,0.0); Complex eye = new Complex(0.0,1.0); Complex meye = new Complex(0.0,-1.0); Complex xp2, xp3, xp4, yp2, yp3, yp4, zp2, zp3, zp4; // precomputed double tol = 1.0e-6; // can be changed public deriv_poly43() // constructor { Complex x, y, z, tx, ty, tz, err, val; Complex bestx, besty, bestz, bestval, tbestval; Complex bestdx, bestdy, bestdz; Complex dx, dy, dz, tdx, tdy, tdz, tbdx, tbdy, tbdz; 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 = new Complex(1.0,1.0); // may be changed y = new Complex(1.0,1.0); z = new Complex(1.0,1.0); xp2 = x.multiply(x); // precomputed for efficency xp3 = xp2.multiply(x); xp4 = xp3.multiply(x); yp2 = y.multiply(y); yp3 = yp2.multiply(y); yp4 = yp3.multiply(y); zp2 = z.multiply(z); zp3 = zp2.multiply(z); zp4 = zp3.multiply(z); val = eval(x, y, z); System.out.println("x="+x+", y="+y+", z="+z); System.out.println("val="+val); // need derivatives fpx, fpy, fpz to build jacobian fpxv = fpx(x, y, z); System.out.println("fpxv="+fpxv); fpyv = fpy(x, y, z); System.out.println("fpyv="+fpyv); fpzv = fpz(x, y, z); System.out.println("fpzv="+fpzv); sfpx(P, Der); prtp("Deriv wrt x", Der); sfpy(P, Der); prtp("Deriv wrt y", Der); sfpz(P, Der); prtp("Deriv wrt z", Der); System.out.println("deriv_poly43 finished"); } void prtp(String name, Complex p[]) { System.out.println(" "); System.out.println(name); System.out.println("[0]c = "+p[0]); System.out.println(" "); } // end prtp 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 read_data() // big ugly input method { 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); // System.out.println("1.0 pwr[0]="+pwr[0]); // System.out.println("x pwr[1]="+pwr[1]); // System.out.println("y pwr[2]="+pwr[2]); // System.out.println("z pwr[3]="+pwr[3]); // System.out.println("xx pwr[4]="+pwr[4]); // System.out.println("yy pwr[7]="+pwr[7]); // System.out.println("zz pwr[9]="+pwr[9]); // System.out.println("xxx pwr[10]="+pwr[10]); // System.out.println("yyy pwr[16]="+pwr[16]); // System.out.println("zzz pwr[19]="+pwr[19]); // System.out.println("xxxx pwr[20]="+pwr[20]); // System.out.println("yyyy pwr[30]="+pwr[30]); // System.out.println("zzzz pwr[34]="+pwr[34]); for(int i=0; i<35; i++) { sum = sum.add(P[i].multiply(pwr[i])); } return sum; } // end eval public static void main (String args[]) { new deriv_poly43(); // construct and execute } } // end deriv_poly34.java