#include <alloc.h>
#include <math.h>

#define PI 3.1415926536
#define TNNearlyZero  1E-20
#define TNArraySize  30   /* maximum size of vectors */

typedef double TNvector;
typedef double *TNvectorPtr;
typedef int TNIntVector;
typedef struct _TNcomplex {
				   double Re, Im;
				 } TNcomplex;
typedef struct _TNquadratic {
				  TNcomplex A, B, C;
				} TNquadratic;

typedef TNcomplex TNCompVector;
typedef TNcomplex *TNCompVectorPtr;

double _abs(double x);

void Newt_Horn_Defl(int            InitDegree  ,
					TNvector	   *InitPoly   ,
					double	       Guess       ,
					double	       Tol         ,
					int	           MaxIter     ,
					int            *Degree     ,
					int            *NumRoots   ,
					TNvector       *Poly       ,
					TNvector       *Root        ,
					TNvector       *Imag        ,
					TNvector       *Value       ,
					TNvector       *Deriv       ,
					TNIntVector    *Iter        ,
					int  *Error);

/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-             Input: InitDegree, InitPoly, Guess, Tol, MaxIter             -*/
/*-            Output: Degree, NumRoots, Poly, Root, Imag, Value, Deriv      -*/
/*-                    Iter, Error                                           -*/
/*-                                                                          -*/
/*-           Purpose: This unit provides a procedure for finding several    -*/
/*-                    roots of a user specified polynomial given a user     -*/
/*-                    specified initial guess.  The procedure implements    -*/
/*-                    Newton-Horner's algorithm for finding a single        -*/
/*-                    root of a polynomial and deflation techniques for     -*/
/*-                    reducing a polynomial given one of its roots.         -*/
/*-                    Should the polynomial contain no more than two        -*/
/*-                    complex roots, they may also be determined.           -*/
/*-                    The user must specify the desired tolerance in the    -*/
/*-                    answer and the maximum number of iterations.          -*/
/*-                                                                          -*/
/*- Pre-Defined Types: TNvector    = array[0..TNArraySize] of real;          -*/
/*-                    TNIntVector = array[0..TNArraySize] of integer;       -*/
/*-                                                                          -*/
/*-  Global Variables: InitDegree : integer;  degree of user's polynomial    -*/
/*-                    InitPoly   : TNvector; coefficients of user's         -*/
/*-                                           polynomial where InitPoly[n]   -*/
/*-                                           is the coefficient of X^n      -*/
/*-                    Guess      : real;     user's estimate of root        -*/
/*-                    Tol        : real;     tolerance in answer            -*/
/*-                    Degree     : real;     degree of reduced polynomial   -*/
/*-                                           left when procedure is done    -*/
/*-                                           (>0 if all the roots were      -*/
/*-                                           not Found)                     -*/
/*-                    Poly       : TNvector; coefficients of reduced poly   -*/
/*-                    NumRoots   : integer;  number of roots calculated     -*/
/*-                    Root       : TNvector; real parts of calculated roots -*/
/*-                    Imag       : TNvector; imaginary parts of roots (non- -*/
/*-                                           zero for no more than 2)       -*/
/*-                    Value      : TNvector; value of the polynomial at     -*/
/*-                                           each root                      -*/
/*-                    Deriv      : TNvector; value of the derivative at     -*/
/*-                                           each root                      -*/
/*-                    Iter       : TNIntVector; number of iterations it     -*/
/*-                                              took to find each root      -*/
/*-                    Error      : byte;     flags if something went wrong  -*/
/*-                                                                          -*/
/*-            Errors: 0: No error                                           -*/
/*-                    1: Iter >= MaxIter                                    -*/
/*-                    2: The slope was zero at some point                   -*/
/*-                    3: Degree <= 0                                        -*/
/*-                    4: Tol <= 0                                           -*/
/*-                    5: MaxIter < 0                                        -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

void SynDiv(int Degree,
			TNvector	 *Poly,
			double	 X,
			TNvector  *NewPoly)

/*--------------------------------------------------------------------------*/
/*- Input:  Degree, Poly                                                   -*/
/*- Output: NewPoly                                                        -*/
/*-                                                                        -*/
/*- This void applies the technique of synthetic division             -*/
/*- to a polynomial, Poly, at the value X.  The kth element of NewPoly     -*/
/*- is the (k-1)th element of the new polynomial. The 0th element          -*/
/*- of the resulting polynomial, NewPoly, is the value of the polynomial   -*/
/*- at X                                                                   -*/
/*--------------------------------------------------------------------------*/

{
 int  Term;

  *(NewPoly + Degree) = *(Poly + Degree);
  for (Term = Degree - 1;Term >= 0;Term--)
	*(NewPoly + Term) = *(NewPoly+Term + 1) * X + *(Poly + Term);
} /* void SynDiv */

void QuadraticFormula(double       A,
					  double	   B,
					  double	   C,
					  double *ReRoot1,
					  double *ImRoot1,
					  double *ReRoot2,
					  double *ImRoot2)

/*----------------------------------------------------------------*/
/*- Input:  A, B, C                                              -*/
/*- Output: ReRoot1, ImRoot1, ReRoot2, ImRoot2                   -*/
/*-                                                              -*/
/*- This void applies the quadratic formula to the equation -*/
/*- AX^2 + BX + C, where A,B,C are real.  It returns the real or -*/
/*- complex roots of the equation                                -*/
/*----------------------------------------------------------------*/

{
 double  Denominator, Discrim;

  Denominator = 2 * A;
  *ReRoot1 = -B / Denominator;
  *ReRoot2 = *ReRoot1;
  Discrim = B * B - 4 * A * C;
  if ( Discrim < 0 )
	{
	  *ImRoot1 = -sqrt(-Discrim) / Denominator;
	  *ImRoot2 =  sqrt(-Discrim) / Denominator;
	}
  else
	{
	  if ( B < 0 )  /* Choose ReRoot1 to have the greatest absolute value */
		*ReRoot1 = *ReRoot1 + sqrt(Discrim) / Denominator;
	  else
		*ReRoot1 = *ReRoot1 - sqrt(Discrim) / Denominator;
	  *ReRoot2 = C / (A * (*ReRoot1));   /* The product of the 2 roots is C/A */
	  *ImRoot1 = 0;
	  *ImRoot2 = 0;
	}
}    /* void QuadraticFormula */

void TestData(int          InitDegree,
			  TNvector     *InitPoly,
			  double	   Tol,
			  int	       MaxIter,
			  int          *Degree,
			  TNvectorPtr  *Poly,
			  int          *NumRoots,
			  TNvector     *Roots,
			  TNvector     *yRoots,
			  TNIntVector  *Iter,
			  int *Error)

/*----------------------------------------------------------*/
/*- Input:  InitDegree, InitPoly, Tol, MaxIter             -*/
/*- Output: Degree, Poly, NumRoots, Roots, yRoots,         -*/
/*-         Iter, Error                                    -*/
/*-                                                        -*/
/*- This void sets the initial value of the above     -*/
/*- variables.  This void also tests the tolerance    -*/
/*- (Tol), maximum number of iterations (MaxIter), and     -*/
/*- Degree for errors and returns the appropriate error    -*/
/*- code.  Finally, it examines the coefficients of Poly.  -*/
/*- If the constant term is zero, ) zero is one of the  -*/
/*- roots and the polynomial is deflated accordingly. Also -*/
/*- if ( the leading coefficient is zero, ) Degree is     -*/
/*- reduced until the leading coefficient is non-zero.     -*/
/*----------------------------------------------------------*/

{
  int Term;

  *Error = 0;
  *NumRoots = 0;
  *Degree = InitDegree;
  *Poly = ( TNvector *) calloc(TNArraySize,sizeof(TNvector));
  for (InitDegree = 0; InitDegree <= *Degree; InitDegree++)
   *(*Poly + InitDegree) =  *(InitPoly + InitDegree);
  if ( Tol <= 0 )
	*Error = 4;
  if ( MaxIter < 0 )
	*Error = 5;
  /* Reduce Degree until leading coefficient <> zero */
  while ( (Degree > 0) && (_abs(*(*Poly + *Degree) < TNNearlyZero)) )
	/* Reduce Degree until leading coefficient <> zero */
	(*Degree)--;
  if ( *Degree <= 0 )
	*Error = 3;
  /* Deflate polynomial until the constant term <> zero */
  while ( (_abs(**Poly) < TNNearlyZero) && (Degree > 0) )
  {
	(*NumRoots)++;
	*(Roots + (*NumRoots)) = 0;
	*(yRoots + (*NumRoots)) = 0;
	*(Iter + (*NumRoots)) = 0;
	(*Degree)--;
	for ( Term = 0; Term <= *Degree; Term++ )
	  *(Poly + Term) = *(Poly + Term + 1);
  }
} /* void TestData */

void FindValueAndDeriv(int Degree,
					  TNvector		*Poly,
					  double		X,
					  double *Value,
					  double *Deriv)

/*--------------------------------------------------------------------*/
/*- Input:  Degree, Poly, X                                          -*/
/*- Output: Value, Deriv                                             -*/
/*-                                                                  -*/
/*- This function applies the technique of synthetic division to     -*/
/*- determine both the Value and derivative of the polynomial, Poly, -*/
/*- at X.  The 0th element of the first synthetic division is the    -*/
/*- value of the polynomial at X, and the 1st element of the second  -*/
/*- synthetic division is the derivative of the polynomial at X.     -*/
/*--------------------------------------------------------------------*/

{
  TNvector  *Poly1, *Poly2;

  Poly1 = ( TNvector *) calloc(TNArraySize,sizeof(TNvector));
  Poly2 = ( TNvector *) calloc(TNArraySize,sizeof(TNvector));

  SynDiv(Degree, Poly, X, Poly1);
  *Value = *(Poly1);
  SynDiv(Degree, Poly1, X, Poly2);
  *Deriv = *(Poly2 + 1);

  free(Poly1);free(Poly2);
} /* void FindValueAndDeriv */

void CheckSlope(double Slope,int *Error)

/*---------------------------------------------------*/
/*- Input:  Slope                                   -*/
/*- Output: Error                                   -*/
/*-                                                 -*/
/*- This void checks the slope to see if ( it is -*/
/*- zero.  The Newton Raphson algorithm may ! be  -*/
/*- applied at a point where the slope is zero.     -*/
/*---------------------------------------------------*/

{
  if ( _abs(Slope) <= TNNearlyZero )
	*Error = 2;
} /* void CheckSlope */

void Initial(int      Degree   ,
			 TNvector	  *Poly     ,
			 double	  Guess    ,
			 double   *OldX     ,
			 double   *OldY     ,
			 double   *OldDeriv ,
			 char     *Found    ,
			 int      *Iter     ,
			 int   *Error)

/*-------------------------------------------------------------*/
/*- Input:  Degree, Poly, Guess                               -*/
/*- Output: OldX, OldY, OldDeriv, Found, Iter, Error          -*/
/*-                                                           -*/
/*- This void sets the initial values of the above       -*/
/*- variables. If OldY is zero, ) a root has been          -*/
/*- found and Found = TRUE.                                   -*/
/*-------------------------------------------------------------*/

{
  *Found = 0;
  *Iter = 0;
  *Error = 0;
  *OldX = Guess;
  FindValueAndDeriv(Degree, Poly, *OldX, OldY, OldDeriv);

  if ( _abs(*OldY) <= TNNearlyZero )
	*Found = 1;
  else
	CheckSlope(*OldDeriv, Error);
} /* void Initial */

int TestForRoot(double X, double OldX, double Y, double Tol)

/*----------------------------------------------------------------*/
/* These are the stopping criteria.  Four different ones are     -*/
/* provided.  If you wish to change the active criteria, simply  -*/
/* comment off the current criteria (including the preceding OR) -*/
/* and remove the comment brackets from the criteria (including  -*/
/* the following OR) you wish to be active.                      -*/
/*----------------------------------------------------------------*/

{
  return                              /*---------------------------*/
	(_abs(Y) <= TNNearlyZero)          /*- Y=0                     -*/
									  /*-                         -*/
		   ||                         /*-                         -*/
									  /*-                         -*/
	(_abs(X - OldX) < _abs(OldX*Tol))   /*- Relative change in X    -*/
									  /*-                         -*/
 /*        ||                      */ /*-                         -*/
 /*                                */ /*-                         -*/
 /* (_abs(OldX - X) < Tol)          */ /*- Absolute change in X    -*/
 /*                                */ /*-                         -*/
 /*        ||                      */ /*-                         -*/
 /*                                */ /*-                         -*/
 /* (_abs(Y) <= Tol)                */ /*- Absolute change in Y    -*/
			;						  /*---------------------------*/

/*-----------------------------------------------------------------------*/
/*- The first criteria simply checks to see if ( the value of the         -*/
/*- function is zero.  You should probably always keep this criteria    -*/
/*- active.                                                             -*/
/*-                                                                     -*/
/*- The second criteria checks the relative error in X. This criteria   -*/
/*- evaluates the fractional change in X between interations. Note      -*/
/*- that X has been multiplied through the inequality to avoid divide   -*/
/*- by zero errors.                                                     -*/
/*-                                                                     -*/
/*- The third criteria checks the absolute difference in X between      -*/
/*- iterations.                                                         -*/
/*-                                                                     -*/
/*- The fourth criteria checks the absolute difference between          -*/
/*- the value of the function and zero.                                 -*/
/*-----------------------------------------------------------------------*/

} /* void TestForRoot */

int FindOneRoot_1_TestForRoot(double X, double OldX, double Y, double Tol)

/*----------------------------------------------------------------*/
/* These are the stopping criteria.  Four different ones are     -*/
/* provided.  If you wish to change the active criteria, simply  -*/
/* comment off the current criteria (including the preceding OR) -*/
/* and remove the comment brackets from the criteria (including  -*/
/* the following OR) you wish to be active.                      -*/
/*----------------------------------------------------------------*/

{
  return                              /*---------------------------*/
	(_abs(Y) <= TNNearlyZero)          /*- Y=0                     -*/
									  /*-                         -*/
		   ||                         /*-                         -*/
									  /*-                         -*/
	(_abs(X - OldX) < _abs(OldX*Tol))   /*- Relative change in X    -*/
									  /*-                         -*/
 /*        ||                      */ /*-                         -*/
 /*                                */ /*-                         -*/
 /* (_abs(OldX - X) < Tol)          */ /*- Absolute change in X    -*/
 /*                                */ /*-                         -*/
 /*        ||                      */ /*-                         -*/
 /*                                */ /*-                         -*/
 /* (_abs(Y) <= Tol)                */ /*- Absolute change in Y    -*/
		   ;							  /*---------------------------*/

/*-----------------------------------------------------------------------*/
/*- The first criteria simply checks to see if ( the value of the       -*/
/*- function is zero.  You should probably always keep this criteria    -*/
/*- active.                                                             -*/
/*-                                                                     -*/
/*- The second criteria checks the relative error in X. This criteria   -*/
/*- evaluates the fractional change in X between interations. Note      -*/
/*- that X has been multiplied through the inequality to avoid divide   -*/
/*- by zero errors.                                                     -*/
/*-                                                                     -*/
/*- The third criteria checks the absolute difference in X between      -*/
/*- iterations.                                                         -*/
/*-                                                                     -*/
/*- The fourth criteria checks the absolute difference between          -*/
/*- the value of the function and zero.                                 -*/
/*-----------------------------------------------------------------------*/

} /* void FindOneRoot_1_TestForRoot */



void FindOneRoot_2(int Degree  ,
				 TNvector	  *Poly    ,
				 double	  Guess   ,
				 double	  Tol     ,
				 int	  MaxIter ,
				 double   *Root    ,
				 double   *Value   ,
				 double   *Deriv   ,
				 int      *Iter    ,
				 int  *Error)

/*-------------------------------------------------------------------*/
/*- Input:  Degree, Poly, Guess, Tol, MaxIter                       -*/
/*- Output: Root, Value, Deriv, Iter, Error                         -*/
/*-                                                                 -*/
/*- A single root of the polynomial Poly.  The root must be         -*/
/*- approximated within MaxIter iterations to a tolerance of Tol.   -*/
/*- The root, value of the polynomial at the root (Value), and the  -*/
/*- value of the derivative of the polynomial at the root (Deriv),  -*/
/*- and the number of iterations (Iter) are returned. If no root    -*/
/*- is found, the appropriate error code (Error) is returned.       -*/
/*-------------------------------------------------------------------*/

{
  char Found;
  double OldX, OldY, OldDeriv, NewX, NewY, NewDeriv;

  /* void FindOneRoot_2 */
  Initial(Degree, Poly, Guess,
					  &OldX, &OldY, &OldDeriv, &Found,
					  Iter, Error);

  while ( !Found && !(*Error) && (*Iter<MaxIter) )
  {
	(*Iter)++;
	NewX = OldX - OldY / OldDeriv;
	FindValueAndDeriv(Degree, Poly, NewX, &NewY, &NewDeriv);

	Found = TestForRoot(NewX, OldX, NewY, Tol);
	OldX = NewX;
	OldY = NewY;
	OldDeriv = NewDeriv;
	if ( !(Found) )
	  CheckSlope(OldDeriv, Error);
  }
  *Root = OldX;
  *Value = OldY;
  *Deriv = OldDeriv;
  if ( !Found && !(*Error) && (*Iter >= MaxIter) )
	*Error = 1;
  if ( Found )
	*Error = 0;
} /* void FindOneRoot_2 */

void ReducePolynomial(int      *Degree ,
					  TNvector *Poly   ,
					  double   Root)

/*------------------------------------------------------*/
/*- Input:  Degree, Poly, Root                         -*/
/*- Output: Degree, Poly                               -*/
/*-                                                    -*/
/*- This function deflates the polynomial Poly by      -*/
/*- factoring out the Root.  Degree is reduced by one. -*/
/*------------------------------------------------------*/

{
  TNvector *NewPoly;
  int Term;

  NewPoly = ( TNvector *) calloc(TNArraySize,sizeof(TNvector));

  SynDiv(*Degree, Poly, Root, NewPoly);
  (*Degree)--;
  for ( Term = 0; Term <= *Degree; Term++ )
	*(Poly + Term) = *(NewPoly + Term+1);

  free(NewPoly);
} /* void ReducePolynomial */


void Newt_Horn_Defl(int            InitDegree  ,
					TNvector	   *InitPoly    ,
					double	       Guess       ,
					double	       Tol         ,
					int	           MaxIter     ,
					int            *Degree     ,
					int            *NumRoots   ,
					TNvector       *Poly        ,
					TNvector       *Root        ,
					TNvector       *Imag        ,
					TNvector       *Value       ,
					TNvector       *Deriv       ,
					TNIntVector    *Iter        ,
					int  *Error)
{
 int  Iter1;


  /* void Newt_Horn_Defl */
  TestData(InitDegree, InitPoly, Tol, MaxIter,
		   Degree, &Poly, NumRoots, Root, Value, Iter, Error);

  while ( !(*Error) && (*Degree>2) )
  {
	FindOneRoot_2(*Degree, Poly, Guess, Tol, MaxIter, Root + (*NumRoots+1),
				Value + (*NumRoots)+1, Deriv + (*NumRoots)+1,
				Iter + (*NumRoots)+1, Error);

	if ( !(*Error) )
	{
	  (*NumRoots)++;
	  /*------------------------------------------------------*/
	  /*- The next statement refines the approximate root by -*/
	  /*- plugging it into the original polynomial.  This    -*/
	  /*- eliminates a lot of the round-off error            -*/
	  /*- accumulated through many iterations                -*/
	  /*------------------------------------------------------*/
	  if ( *NumRoots > 1 )
	  {
		Iter1 = 0;
		FindOneRoot_2(InitDegree, InitPoly, *(Root + (*NumRoots)),
					Tol, MaxIter, Root + (*NumRoots), Value + (*NumRoots),
					Deriv + (*NumRoots), &Iter1, Error);
		*(Iter + (*NumRoots)) = *(Iter + (*NumRoots)) + Iter1;
	  }
	  ReducePolynomial(Degree, Poly, *(Root + (*NumRoots)));
	  Guess = *(Root + (*NumRoots));
	}
  }
  switch (*Degree) {
	case 1 :            /* Solve this linear */
			*Degree = 0;
			(*NumRoots)++;
			Root[*NumRoots] = -Poly[0] / Poly[1];
			FindOneRoot_2(InitDegree, InitPoly, Root[*NumRoots], Tol,
						MaxIter, Root + (*NumRoots), Value + (*NumRoots),
						Deriv + (*NumRoots), Iter + (*NumRoots), Error);
			break;

	case 2 :           /* Solve this quadratic */
			*Degree = 0;
			*(NumRoots) += 2;
			QuadraticFormula(Poly[2], Poly[1], Poly[0],
							 Root + (*NumRoots) - 1, Imag + (*NumRoots) - 1,
							 Root + (*NumRoots), Imag + (*NumRoots));
			if ( _abs(*(Imag + (*NumRoots))) < TNNearlyZero )
			  /* if ( the roots are real, they can be     */
			  /* made more accurate using Newton-Horner */
			  {
				FindOneRoot_2(InitDegree, InitPoly, *(Root + (*NumRoots)-1), Tol,
							MaxIter, Root + (*NumRoots-1), Value + (*NumRoots)-1,
							Deriv + (*NumRoots)-1, Iter + (*NumRoots)-1, Error);

				FindOneRoot_2(InitDegree, InitPoly, *(Root + (*NumRoots)), Tol,
							MaxIter, Root + (*NumRoots), Value + (*NumRoots),
							Deriv + (*NumRoots), Iter + (*NumRoots), Error);
			  }
			else
			  /* If the roots are complex, ) assign    */
			  /* the value to be zero (which is true,     */
			  /* except for some roundoff error) and the  */
			  /* derivative to be zero (which is usually  */
			  /* FALSE; the derivative is usually complex */
			  {
				*(Value + (*NumRoots)-1) = 0;    *(Value + (*NumRoots)) = 0;
				*(Deriv + (*NumRoots)-1) = 0;    *(Deriv + (*NumRoots)) = 0;
				*(Iter + (*NumRoots)-1) = 0;     *(Iter + (*NumRoots)) = 0;
			  }

  } /* switch */
} /* void Newt_Horn_Defl */


double _abs(double x)
{
 return ( x < 0 ? -x : x);
}
