/**************************************************** * Program to demonstrate the BAIRSTOW subroutine * * ------------------------------------------------- * * Reference: BASIC Scientific Subroutines, Vol. II * * By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * * * * C++ version by J-P Moreau, Paris. * * (www.jpmoreau.fr) * * ------------------------------------------------- * * Example: Find two complex roots of polynomial: * * f(x) = x^5-10x^4+35x^3-50x^2+24x * * * * SAMPLE RUN: * * * * Input order of polynomial: 5 * * * * Input the polynomial coefficients: * * * * A(0) = 0 * * A(1) = 24 * * A(2) = -50 * * A(3) = 35 * * A(4) = -10 * * A(5) = 1 * * * * Convergence factor: 0.000001 * * Maximum number of iterations: 50 * * * * The roots found are: * * * * X1 = 1.000000 * * Y1 = 0.000000 * * * * X2 = 0.000000 * * Y2 = 0.000000 * * * * The number of iterations was: 11 * * * ****************************************************/ #include #include double A[10],B[10],C[10],D[10]; int i,k,m,n; double aa,bb,e,PI,x1,x2,yy1,y2; /************************************************* * Bairstow complex root subroutine * * ---------------------------------------------- * * This routine finds the complex conjugate roots * * of a polynomial having real coefficients. * * ---------------------------------------------- * * Reference: Computer Methods for Science and * * Engineering by R.L. Lafara. * * ---------------------------------------------- * * INPUTS: * * Polynomial coefficients : A(0) to A(m) * * Order of polynomial (>=4) : m * * Initial guess : a and b * * Convergence factor : e * * Maximum number of iterations : n * * OUTPUTS: * * Two conjugate complex roots : x1,y1 x2,y2 * * Number of iterations : k * *************************************************/ void Bairstow() { // Labels: e100,e200,e300,e400 int i,j; double a1,b1,cc,dd,d2; //Normalize the A[i] series for (i=0; i=n) goto e200; //Test for convergence if (fabs(a1)+fabs(b1)>e*e) goto e100; //Extract roots from quadratic equation e200: cc=aa*aa-4.0*bb; //Test to see if a complex root if (cc>0) goto e300; x1=-aa; x2=x1; yy1=sqrt(-cc); y2=-yy1; goto e400; e300: x1=-aa+sqrt(cc); x2=-aa-sqrt(cc); yy1=0; y2=yy1; e400: x1=x1/2.0; x2=x2/2.0; yy1=yy1/2.0; y2=y2/2.0; } void main() { PI = 4*atan(1); printf("\n Order of polynomial: "); scanf("%d",&m); printf("\n Input the polynomial coefficients:\n\n"); for (i=0; i