/******************************************************* * Equations of degree 2 to 4 * * ---------------------------------------------------- * * This program calculates the real or complex roots of * * algebraic equations of degree 2, 3 and 4. * * ---------------------------------------------------- * * Reference: Mathématiques et statitiques - Programmes * * en BASIC. Editions du P.S.I., 1981 * * [BIBLI 13]. * * * * C++ version by J-P Moreau, Paris * * (www.jpmoreau.fr) * * ---------------------------------------------------- * * SAMPLE RUN: * * * * ROOTS OF EQUATIONS OF DEGREE 2, 3 OU 4 * * Example: X^4+4X^2-3X+3 = 0 * * * * INSTRUCTIONS: * * * * 1. Input degree of equation * * 2. Input coefficients beginning with highest exponent* * 3. Solutions are given under the form: * * (real part) + I (imaginary part) * * * * Degree of equation (2 to 4): 4 * * * * Coefficient of X^4 = 1 * * Coefficient of X^3 = 0 * * Coefficient of X^2 = 4 * * Coefficient of X^1 = -3 * * Coefficient of X^0 = 3 * * * * Root 1 = 0.449702657 + I 0.731070322 * * Root 2 = 0.449702657 - I 0.731070322 * * Root 1 = -0.449702657 + I 1.967231848 * * Root 2 = -0.449702657 - I 1.967231848 * * * *******************************************************/ #include #include double a[5][5], r[5]; int i,n; double aa,b,c,d,ii,im,k,l,m,q,rr,s,sw; double f(double x) { return x * x * x + a[3][2] * x * x + a[3][1] * x + a[3][0]; } /****************************************** * This Subroutine calculates the roots of * * X^2 + A(2,1)*X + A(2,0) = 0 * * W = Determinant of equation * ******************************************/ void Root_2() { double q1, q2, w; w = a[2][1] * a[2][1] - 4 * a[2][0]; q1 = -a[2][1] / 2.0; if (w < 0) { q2 = sqrt(-w) / 2.0; im = q2; r[1] = q1; r[2] = q1; } else if (w==0) { r[1] = q1; r[2] = q1; im = 0.0; } else if (w > 0) { q2 = sqrt(w) / 2.0; im = 0.0; r[1] = q1 + q2; r[2] = q1 - q2; } } /******************************************* * This subroutine calculates the roots of * * X^3 + A(3,2)*X^2 + A(3,1)*X + A(3,0) = 0 * *******************************************/ void Root_3() { //Labels e100,e200,e500 double am, er, te, tt, xa, xb, xc, y1, y2; // one root equals zero if (a[3][0] == 0) { xc = 0; goto e500; } // looking for( maximum coefficient in fabsolute magnitude am = fabs(a[3][0]); for (i = 1; i<4; i++) { tt = fabs(a[3][i]); if (am < tt) am = tt; } //Define interval where a real root exists // according to sign of A(3,0) if (a[3][0] > 0) { aa = -am - 1; b = 0; goto e100; } aa = 0; b = am + 1; e100: // Searching for( xc = real root in interval (a,b) // (by Bisection method) // Define segment (xa,xb) including the root xa = aa; xb = b; y1 = f(aa); y2 = f(b); // Desired précision er = 0.000001; if (y1 == 0) { xc = aa; goto e500; } if (y2 == 0) { xc = b; goto e500; } // xc : middle of segment e200: xc = (xa + xb) / 2.0; te = f(xc); if (te == 0) goto e500; if ((xb - xa) < er) goto e500; if (f(xa) * te > 0) { xa = xc; goto e200; } xb = xc; goto e200; // r[3] is a real root e500: r[3] = xc; if (sw = -1) return; // Calculates the roots of remaining quadratic equation // Define the equation coefficients a[2][1] = a[3][2] + xc; a[2][0] = a[3][1] + a[3][2] * xc + xc * xc; Root_2(); return; } // Root search main subroutine void Root_4() { // Labels: e100,e200,e300 // Normalization of coefficients for (i = 0; i= 0)) goto e100; // Particular case when this equation is of 2nd order a[2][1] = a[3][2]; a[2][0] = a[3][1]; Root_2(); // One takes the positive root r[3] = r[1]; goto e200; // Calling Root_3 with sw=-1 to calculate one root only e100: sw = -1; Root_3(); // real root of above equation e200: k = sqrt(r[3]); // Calculate L and M if k=0 if (k == 0) { rr = sqrt(q * q - 4 * s); goto e300; } q = q + (4 * r[3]); rr = rr / (2 * k); e300: l = (q - rr) / 2.0; m = (q + rr) / 2.0; // Solving two equations of degree 2 a[2][1] = 2 * k; a[2][0] = l; // 1st equation Root_2(); // Transferring solutions in r[3], r[4], ii r[3] = r[1] - (a[4][3] / 4); r[4] = r[2] - (a[4][3] / 4.0); ii = im; a[2][1] = (double) -2 * k; a[2][0] = m; // 2nd equation Root_2(); r[2] = r[2] - (a[4][3] / 4.0); r[1] = r[1] - (a[4][3] / 4.0); } } void main() { printf("\n ROOTS OF EQUATIONS OF DEGREE 2, 3 OR 4\n"); printf(" EXAMPLE : X^4+4X^2-3X+3 = 0\n\n"); printf(" INSTRUCTIONS:\n\n"); printf(" 1. Input degree of equation\n"); printf(" 2. Input coefficients beginning by highest exponent\n"); printf(" 3. Solutions are given under the form:\n"); printf(" (real part) + I (imaginary part)\n\n"); printf(" DEGREE OF EQUATION [2 to 4]: "); scanf("%d",&n); if (n < 2) n = 2; if (n > 4) n = 4; printf("\n"); for (i = n; i>-1; i--) { printf(" COEFFICIENT OF X^%d = ",i); scanf("%lf",&a[n][i]); } printf("\n"); // calling root search main subroutine Root_4(); // writing results if (n==2) { //case n=2 if (im == 0) { printf(" Root 1 = %12.9f\n",r[1]); printf(" Root 2 = %12.9f\n",r[2]); } else { printf(" Root 1 = %12.9f +I %12.9f\n",r[1],im); printf(" Root 2 = %12.9f -I %12.9f\n",r[2],im); } } else if (n==3) { //case n=3 printf(" Root 1 = %12.9f\n",r[3]); if (im == 0) { printf(" Root 2 = %12.9f\n",r[1]); printf(" Root 3 = %12.9f\n",r[2]); } else { printf(" Root 2 = %12.9f +I %12.9f\n",r[1],im); printf(" Root 3 = %12.9f -I %12.9f\n",r[2],im); } } else if (n==4) { //case n=4 if (im==0) { printf(" Root 1 = %12.9f\n",r[1]); printf(" Root 2 = %12.9f\n",r[2]); } else { printf(" Root 1 = %12.9f +I %12.9f\n",r[1],im); printf(" Root 2 = %12.9f -I %12.9f\n",r[2],im); } if (ii==0) { printf(" Root 3 = %12.9f\n",r[3]); printf(" Root 4 = %12.9f\n",r[4]); } else { printf(" Root 3 = %12.9f +I %12.9f\n",r[3],ii); printf(" Root 4 = %12.9f -I %12.9f\n",r[4],ii); } } printf("\n"); } // End of file root4.cpp