/****************************************************
*    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 <stdio.h>
#include <math.h>

#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<m+1; i++)
    C[i]=A[i]/A[m];
  // Start iteration
  // Set initial guess for the quadratic coefficients
  B[0]=0; B[1]=0;
e100: B[m-1]=C[m-1]-aa;
  B[m-2]=C[m-2]-aa*B[m-1]-bb;
  for (j=3; j<m+1; j++)
    B[m-j]=C[m-j]-aa*B[m+1-j]-bb*B[m+2-j];
  //Guard against divide by zero
  if (B[2]!=0) goto e200;
  aa=aa+1e-7; bb=bb-1e-7;
e200: a1=(C[1]-bb*B[3])/B[2];
  b1=C[0]/B[2];
  k++;
  // Test for the number of iterations
  if (k>=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<m+1; i++) {
    printf("     A(%d) = ",i); scanf("%lf",&A[i]);
  }
  printf("\n Convergence factor: "); scanf("%lf",&e);
  printf("\n Maximum number of iterations: "); scanf("%d",&n);

  aa=PI; bb=sqrt(2.0);

  LIN();       // Call LIN subroutione

  printf("\n\n The roots found are:\n\n");
  printf("      X1 = %f\n", x1);
  printf("      Y1 = %f\n\n", yy1);
  printf("      X2 = %f\n", x2);
  printf("      Y2 = %f\n\n", y2);
  printf(" The number of iterations was: %d\n\n", k);
}

// End of file Lin.cpp