/**************************************************** * Program to demonstrate the LIN subroutine * * ------------------------------------------------- * * Reference: BASIC Scientific Subroutines, Vol. II * * by F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981. * * * * C++ version by J-P Moreau, Paris. * * ------------------------------------------------- * * Example: Find two complex roots of polynomial: * * f(x) = x^5+5x^4+10x^3+10x^2+5x+1 * * * * SAMPLE RUN: * * * * Input order of polynomial: 5 * * * * Input the polynomial coefficients: * * * * A(0) = ? 1 * * A(1) = ? 5 * * A(2) = ? 10 * * A(3) = ? 10 * * A(4) = ? 5 * * A(5) = ? 1 * * * * Convergence factor: 0 * * Maximum number of iterations: 1000 * * * * The roots found are: * * * * X1 = -0.937617 * * Y1 = 0.025812 * * * * X2 = -0.937617 * * Y2 = -0.025812 * * * * The number of iterations was: 1000 * * * ****************************************************/ #include #include #define PI 4*atan(1) double A[10],B[10],C[10]; int i,k,m,n; double aa,bb,e,x1,x2,yy1,y2; /************************************************ * Polynomial complex roots subroutine * * --------------------------------------------- * * This routine uses the Lin's method * * --------------------------------------------- * * Reference: A Practical Guide to Computer * * Method for Engineers by T.E. Shoup * * --------------------------------------------- * * INPUTS: * * Polynomial coefficients : A(0) to A(m) * * Order of polynomial : m * * Initial guess : a and b * * Convergence factor : e * * Maximum number of iterations : n * * OUTPUTS: * * Two complex roots : x,y1 x2,y2 * * Number of iterations : k * ************************************************/ void LIN() { //Labels: e100,e200,e300,e400,e500 int i,j; double a1,b1,cc; // Normalize the A(i) series for (i=0; i=n) goto e300; // Test for convergence if (fabs(aa-a1)+fabs(bb-b1) < e*e) goto e300; aa=a1; bb=b1; // Next iteration goto e100; e300: aa=a1; bb=b1; cc=aa*aa-4.0*bb; // Is there an imaginary part? if (cc>0) goto e400; yy1=sqrt(-cc); y2=-yy1; x1=-aa; x2=x1; goto e500; e400: yy1=0; y2=yy1; x1=-aa+sqrt(cc); x2=-aa-sqrt(cc); e500: x1=x1/2.0; x2=x2/2.0; yy1=yy1/2.0; y2=y2/2.0; } //LIN() void main() { printf("\n Order of polynomial: "); scanf("%d",&m); printf("\n Input the polynomial coefficients:\n\n"); for (i=0; i