#include
#include

#define u_0 2
#define v_0 2
#define tol 0.0000001
/* tolerance in required Procedure Newton */

/* master 040195 */


typedef double coeffType[4] ;  /* 0..4 */

typedef coeffType rootsType ;

enum enumTypeRoot { realr , imagr } ;
typedef enum enumTypeRoot typeRoot ;

static struct TypeRec {
    /* represents a root, complex or real */
    double a , b;
    typeRoot root ;
} ;

typedef struct TypeRec typeRec ;

/* repair 050195 */
/* typedef typeRec rootArr[2][2] ; */ /*keeps all the roots of the polynomial*/


static void Newton ( double *u, double *v , double r, double s,
	      double c_0,double c_1,double c_2)
{
/* Finds next iterates of u and v */

double denom ;

 denom=c_1*c_1 - c_0*c_2;
 if ( denom==0 )
 {
  denom=0.00001 ;
 }

 *u=*u-(r*c_1 - s*c_2)/denom;
 *v=*v-(s*c_1 - r*c_0)/denom;

} /* NEWTON */



static void FindFact(double a[] , double b[] , double *u, double *v,double *r,
	      double *s , int *count, int *n)
{
/* Finds a quadratic factor of the polynomial using Bairstow's algorithm */

 double oldv, oldu ;
 double c[5] ;
 int stop ; /* Boolean */
 double r_u, r_v, s_u, s_v ;
 int index ;
 double tmp;

 stop = -1 ; /* false */
 *count=0;

 *u=u_0;
 *v=v_0;

 do
 {
   /* get b's, r and s */

     b[*n-2] = a[*n];
     b[*n-3] = a[*n-1] + (*u)*b[*n-2];
     b[*n-3] = a[*n-1] + (*u)*b[*n-2];

     for (index=*n-4 ; index >= 0 ; index = index - 1 )
     {
     b[index]=a[index+2] + *u * b[index+1] + * v * b[index+2] ;
     } /* ??? */

     *r=a[1] + *u * b[0] + *v * b[1] ;
     *s=a[0] + *u * *r + *v * b[0] ;

   /* end get b's, r and s */
   /* get c's */

   c[*n-1]= b[*n-2];
   c[*n-2]= b[*n-3] + *u * c[*n-1];

   for ( index=*n-3 ; index >= 1 ; index = index - 1 )
   {
     c[index]=b[index-1] + *u * c[index+1] + *v * c[index+2];
   } ; /* ??? */

     c[0]=*r + *u * c[1] + *v * c[2] ;

    /* end get c's */

   /* get partials */

 /*  r_u=c[1];
     s_u=c[0];
     r_v=c[2];
     s_v=c[1]; */

   /* end get partials */

   /* find next iterates of u and v */

    oldu=*u ;
    oldv=*v ;

    Newton(u,v,*r,*s,c[0],c[1],c[2]);

   /* end find */

   if ( (fabs(*u-oldu) < tol ) && (fabs(*v-oldv)=-0.000000001) && (delta <= 0.0) )
 {
  delta=0.0;
 }

 if (delta >= 0.0)
 { /* real roots */
  roots[*num][1].a=(u + sqrt(delta))/2;
  roots[*num][1].b=0;
  roots[*num][2].a=(u - sqrt(delta))/2;
  roots[*num][2].b=0;
  roots[*num][1].root=realr;
  roots[*num][2].root=realr ;
 } else {
   /* delta<0 */
   roots[*num][1].root=imagr;
   roots[*num][2].root=imagr;
   roots[*num][1].a=u/2;
   roots[*num][2].a=roots[*num][1].a;
   roots[*num][1].b=-sqrt(fabs(delta))/2;
   roots[*num][2].b=-roots[*num][1].b ;
  } ;

} /* FindRoot */ ;



void QuadRoots ( double *p , int *numReal , double realRoots[] )
{
 double u, v, oldv, oldu, r, s ;
 double a[5], b[5], c[5] ;
 int index, count, num, n ;
 int  one, already ; /* Boolean */
 double  r_u, r_v, s_u, s_v ;
 typeRec roots[3][3];
 double tmpRoots[5];
 int i , j , duplicate , rootsPnt , numTmp;


 n=4 ;

 for ( index=0 ; index<= n ; index++ )
 {
    a[index] = *(p+index) ;
 } ;

  num=0;
 
 /* Find all roots */

 while (n>2)
 {
  num++;   /* keeps track of the number of pairs of roots */

  FindFact(&(a[0]),&(b[0]),&u,&v,&r,&s,&count,&n);

  FindRoot(&num,u,v,roots);

  n=n-2;
  for (index=0 ; index<=n ; index++ )
  {
  a[index]=b[index];
  } ;

 } ; /* END WHILE */


 /* Find roots of remaining quadratic or linear term */

 num=num+1;

 one=-1 ; /* false , specifies whether there is a single root or a pair */

 if ( (n==2) && (a[2] != 0) )
 {
  FindRoot(&num,-a[1]/a[2],-a[0]/a[2],roots) ;
 } else
  /* repair */
  if ( (n==1)||(a[2]==0) )
  {

   if (a[1]==0.0)
   {
   num = num-1 ;
   } else
   {
    roots[num][1].a=-a[0]/a[1];
    roots[num][1].b=0;
    roots[num][1].root=realr;
    one=+1 ; /* true */
   }

  }

 *numReal=0 ;

 for ( index=1 ; index <= num ; index++ )
 {
  if (roots[index][1].root==realr)
  {
     if ( (one==1) && (index==num) )
     {
      (*numReal)++ ;
      realRoots[*numReal] = roots[index][1].a ;
     } else
     {
      (*numReal)++ ;
      realRoots[*numReal]= roots[index][1].a ;
      (*numReal)++ ;
      realRoots[*numReal] = roots[index][2].a ;
     } ; /* END ELSE */

  } ; /* ENDIF */

 } ; /* END FOR */


 /* Eliminate repeated roots. */

 numTmp=0;
 for ( i=1 ; i<=*numReal-1 ; i++ )
 {
  duplicate=-1; /*false*/
  for ( j=i+1 ; j<=*numReal ; j++ )
  {
    if ( fabs(realRoots[i]-realRoots[j]) < 0.00001 )
     {
      duplicate=1; /*true*/
     }
  }
  if (duplicate==-1)
  {
   numTmp++;
   tmpRoots[numTmp]=realRoots[i];
  }

 }

 if(*numReal>0){
  numTmp++;
  tmpRoots[numTmp]=realRoots[*numReal];
 }

 if (*numReal != numTmp)
 {
  *numReal=numTmp;
  for ( i=1 ; i<=*numReal ; i++ )
  realRoots[i] = tmpRoots[i] ;
 }


} ; /* END  QuadRoots */