/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-     Turbo C Numerical Methods Library                                    -*/
/*-     Copyright (c) 1986, 87 by Borland International, Inc. (kecy)               -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

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

#define PI 3.1415926536
#define TNNearlyZero  1e-16

typedef double TNvector;
typedef double *TNvectorPtr;

typedef int TNIntVector;
typedef int *TNIntVectorPtr;



typedef struct _TNcomplex {
				   double Re, Im;
				 } TNcomplex;

typedef TNcomplex TNCompVector;
typedef TNcomplex *TNCompVectorPtr;

typedef struct _TNquadratic {
				  TNcomplex A, B, C;
				} TNquadratic;

double UserFunction(double X,double (*FuncAddr)(double));
void  UserProcedure(TNcomplex X, TNcomplex *Y,
					TNcomplex (*ProcAddr)(TNcomplex));
double _abs(double x);

void Laguerre(int               *Degree    ,
			  TNCompVectorPtr  	 *Poly      ,
			  TNcomplex		     InitGuess  ,
			  double		     Tol        ,
			  int		         MaxIter    ,
			  int	             *NumRoots  ,
			  TNCompVectorPtr    *Roots     ,
			  TNCompVectorPtr    *yRoots    ,
			  TNIntVectorPtr	 *Iter      ,
			  int     *Error);



/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-            Input: Degree, Poly, InitGuess, Tol, MaxIter                  -*/
/*-           Output: Degree, Poly, NumRoots, Roots, yRoots, Iter, Error     -*/
/*-                                                                          -*/
/*-          Purpose: This unit provides a procedure for finding all the     -*/
/*-                   roots (real and complex) to a polynomial.              -*/
/*-                   Laguerre's method with deflation is used.              -*/
/*-                   The user must input the coefficients of the quadratic  -*/
/*-                   and the tolerance in the answers generated.            -*/
/*-                                                                          -*/
/*-  Pre-defined Types: TNcomplex    = record                                -*/
/*-                                      Re, Im : real;                      -*/
/*-                                    end;                                  -*/
/*-                                                                          -*/
/*-                     TNIntVector  = array[0..TNArraySize] of integer;     -*/
/*-                     TNCompVector = array[0..TNArraySize] of TNcomplex;   -*/
/*-                                                                          -*/
/*- Global Variables: Degree    : integer;      degree of deflated           -*/
/*-                                             polynomial                   -*/
/*-                   Poly      : TNCompVector; coefficients of deflated     -*/
/*-                                             polynomial where Poly[n] is  -*/
/*-                                             the coefficient of X^n       -*/
/*-                   InitGuess : TNcomplex;    initial guess to a root      -*/
/*-                                             (may be very crude)          -*/
/*-                   Tol       : real;         tolerance in the answer      -*/
/*-                   MaxIter   : integer;      number of iterations         -*/
/*-                   NumRoots  : integer;      number of roots calculated   -*/
/*-                   Roots     : TNCompVector; the roots calculated         -*/
/*-                   yRoots    : TNCompVector; the value of the function    -*/
/*-                                             at the calculated roots      -*/
/*-                   Iter      : TNIntVector;  number iteration it took to  -*/
/*-                                             find each root               -*/
/*-                   Error     : byte;         flags an error               -*/
/*-                                                                          -*/
/*-           Errors: 0: No error                                            -*/
/*-                   1: Iter > MaxIter                                      -*/
/*-                   2: Degree <= 0                                         -*/
/*-                   3: Tol <= 0                                            -*/
/*-                   4: MaxIter < 0                                         -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/



/*----------- Here are a few complex operations ------------*/

void _Conjugate(TNcomplex *C1, TNcomplex *C2)
{
  C2->Re = C1->Re;
  C2->Im = -C1->Im;
} /* void _Conjugate */

double _Modulus(TNcomplex *C1)
{
  return sqrt(pow(C1->Re,2) + pow(C1->Im,2));
} /* function _Modulus */

void _Add(TNcomplex *C1, TNcomplex *C2, TNcomplex *C3)
{
  C3->Re = C1->Re + C2->Re;
  C3->Im = C1->Im + C2->Im;
} /* void _Add */

void _Sub(TNcomplex *C1, TNcomplex *C2, TNcomplex *C3)
{
  C3->Re = C1->Re - C2->Re;
  C3->Im = C1->Im - C2->Im;
} /* void _Sub */

void _Mult(TNcomplex *C1, TNcomplex *C2, TNcomplex *C3)
{
  C3->Re = C1->Re * C2->Re - C1->Im * C2->Im;
  C3->Im = C1->Im * C2->Re + C1->Re * C2->Im;
} /* void _Mult */

void _Divide(TNcomplex *C1, TNcomplex *C2, TNcomplex *C3)
{
  TNcomplex Dum1, Dum2;
  double E;

  _Conjugate(C2, &Dum1);
  _Mult(C1, &Dum1, &Dum2);
  E = pow(_Modulus(C2),2);
  C3->Re = Dum2.Re / E;
  C3->Im = Dum2.Im / E;
}  /* void _Divide */

void Laguerre_SquareRoot(TNcomplex *C1, TNcomplex *C2)
{
  double NearlyZero = 1E-015, R, Theta;

  R = sqrt(pow(C1->Re,2) + pow(C1->Im,2));
  if ( _abs(C1->Re) < NearlyZero )
	{
	  if ( C1->Im < 0 )
		Theta = PI / 2;
	  else
		Theta = -PI / 2;
	}
  else
	if ( C1->Re < 0 )
	  Theta = atan(C1->Im / C1->Re) + PI;
	else
	  Theta = atan(C1->Im / C1->Re);
  C2->Re = sqrt(R) * cos(Theta / 2);
  C2->Im = sqrt(R) * sin(Theta / 2);
} /* void Laguerre_SquareRoot */

void InitAndTest(int              *Degree      ,
				 TNCompVector     *Poly       ,
				 double		      Tol          ,
				 int		      MaxIter      ,
				 TNcomplex		  InitGuess    ,
				 int	          *NumRoots   ,
				 TNCompVector	  *Roots      ,
				 TNCompVector	  *yRoots     ,
				 TNIntVector	  *Iter       ,
				 TNcomplex	      *GuessRoot  ,
				 int	          *InitDegree ,
				 TNCompVectorPtr  *InitPoly   ,
				 int	  *Error)

/*----------------------------------------------------------*/
/*- Input:  Degree, Poly, Tol, MaxIter, InitGuess          -*/
/*- Output: InitDegree, InitPoly, Degree, Poly, NumRoots,  -*/
/*-         Roots, yRoots, Iter, GuessRoot, Error          -*/
/*-                                                        -*/
/*- This void sets the initial value of the above     -*/
/*- variables.  This void also tests the tolerance    -*/
/*- (Tol), maximum number of iterations (MaxIter), and     -*/
/*- 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;
  if ((*Degree) <= 0 )
	*Error = 2;      /* degree is less than 2 */
  if ( Tol <= 0 )
	*Error = 3;
  if ( MaxIter < 0 )
	*Error = 4;

  if (!(*Error))
  {
	*NumRoots = 0;
	*GuessRoot = InitGuess;
	*InitDegree = *Degree;

	*InitPoly = ( TNCompVector * ) calloc((*InitDegree) + 1,sizeof(TNcomplex));
	for (Term = 0; Term <= *InitDegree; Term++)
	  *(*(InitPoly) + Term) = *(Poly + Term);

	/* Reduce degree until leading coefficient <> zero */
	while ( (*Degree > 0) && (_Modulus((Poly + (*Degree))) < TNNearlyZero) )
	  (*Degree)--;
	/* Deflate polynomial until the constant term <> zero */
	while ( (_Modulus(Poly) == 0) && (*Degree > 0) )
	{
	  /* Zero is a root */
	  (*NumRoots)++;
	  (Roots + (*NumRoots))->Re = 0;
	  (Roots + (*NumRoots))->Im = 0;
	  (yRoots + (*NumRoots))->Re = 0;
	  (yRoots + (*NumRoots))->Im = 0;
	  *(Iter + (*NumRoots)) = 0;
	  (*Degree)--;
	  for (Term = 0; Term <= *Degree; Term++)
		*(Poly + Term) = *(Poly + Term + 1);
	}
  }
} /* InitAndTest */

void EvaluatePoly(int              Degree        ,
				  TNCompVector	   *Poly          ,
				  TNcomplex	       X             ,
				  TNcomplex	       *yValue       ,
				  TNcomplex	       *yPrime       ,
				  TNcomplex	       *yDoublePrime)

/*--------------------------------------------------------------------*/
/*- Input:  Degree, Poly, X                                          -*/
/*- Output: yValue, yPrime, yDoublePrime                             -*/
/*-                                                                  -*/
/*- This void applies the technique of synthetic division to    -*/
/*- determine value (yValue), first derivative (yPrime) and second   -*/
/*- derivative (yDoublePrime) of the  polynomial, Poly, at X.        -*/
/*- The 0th element of the first synthetic division is the           -*/
/*- value of Poly at X, the 1st element of the second synthetic      -*/
/*- division is the first derivative of Poly at X, and twice the     -*/
/*- 2nd element of the third synthetic division is the second        -*/
/*- derivative of Poly at X.                                         -*/
/*--------------------------------------------------------------------*/

{
  int Loop;
  TNcomplex Dummy, yDPdummy;
  TNCompVector *Deriv = NULL, *Deriv2 = NULL;

  Deriv = ( TNCompVector * ) calloc(Degree + 1,sizeof(TNcomplex));
  Deriv2 = ( TNCompVector * ) calloc(Degree + 1,sizeof(TNcomplex));

  *(Deriv + Degree) = *(Poly + Degree);
  for ( Loop = Degree - 1; Loop >= 0; Loop-- )
  {
	_Mult(Deriv + Loop + 1, &X, &Dummy);
	_Add(&Dummy, Poly + Loop, Deriv + Loop);
  }
  *yValue = *Deriv;    /* Value of Poly at X */

  *(Deriv2 + Degree) = *(Deriv + Degree);
  for ( Loop = Degree - 1; Loop >= 1; Loop-- )
  {
	_Mult(Deriv2 + Loop + 1, &X, &Dummy);
	_Add(&Dummy, Deriv + Loop, Deriv2 + Loop);
  }
  *yPrime = *(Deriv2 + 1);   /* 1st deriv. of Poly at X */

  yDPdummy = *(Deriv2 + Degree);
  for ( Loop = Degree - 1; Loop >= 2; Loop-- )
  {
	_Mult(&yDPdummy, &X, &Dummy);
	_Add(&Dummy, Deriv2 + Loop, &yDPdummy);
  }
  yDoublePrime->Re = 2 * yDPdummy.Re;    /* 2nd derivative of Poly at X */
  yDoublePrime->Im = 2 * yDPdummy.Im;

  free(Deriv);free(Deriv2);
} /* void EvaluatePoly */

void ConstructDifference(int          Degree       ,
						 TNcomplex	  yValue       ,
						 TNcomplex	  yPrime       ,
						 TNcomplex	  yDoublePrime ,
						 TNcomplex    *Dif)

/*------------------------------------------------------------------*/
/*- Input:  Degree, yValue, yPrime, yDoublePrime                   -*/
/*- Output: Dif                                                    -*/
/*-                                                                -*/
/*- This void computes the difference between approximations; -*/
/*- given information about the function and its first two         -*/
/*- derivatives.                                                   -*/
/*-----------------------------------------------------------------*/

{
  TNcomplex yPrimeSQR, yTimesyDPrime, Sum, SRoot;
  TNcomplex Numer1, Numer2, Numer, Denom;


  _Mult(&yPrime, &yPrime, &yPrimeSQR);
  yPrimeSQR.Re = pow( ( double ) ( Degree - 1 ) ,2) * yPrimeSQR.Re;
  yPrimeSQR.Im = pow( ( double ) ( Degree - 1 ) ,2) * yPrimeSQR.Im;
  _Mult(&yValue, &yDoublePrime, &yTimesyDPrime);
  yTimesyDPrime.Re = (Degree - 1) * Degree * yTimesyDPrime.Re;
  yTimesyDPrime.Im = (Degree - 1) * Degree * yTimesyDPrime.Im;
  _Sub(&yPrimeSQR, &yTimesyDPrime, &Sum);
  Laguerre_SquareRoot(&Sum, &SRoot);
  _Add(&yPrime, &SRoot, &Numer1);
  _Sub(&yPrime, &SRoot, &Numer2);
  if ( _Modulus(&Numer1) > _Modulus(&Numer2) )
	Numer = Numer1;
  else
	Numer = Numer2;
  Denom.Re = Degree * yValue.Re;
  Denom.Im = Degree * yValue.Im;
  if ( _Modulus(&Numer) < TNNearlyZero )
	{
	  Dif->Re = 0;
	  Dif->Im = 0;
	}
  else
	_Divide(&Denom, &Numer, Dif);  /* The difference is the   */
								/* inverse of the fraction */
} /* void ConstructDifference */

int _TestForRoot(double X, double Dif, 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 appropriate OR) -*/
/*-  and remove the comment brackets from the criteria (including    -*/
/*-  the appropriate OR) you wish to be active.                      -*/
/*--------------------------------------------------------------------*/

{
  return                              /*---------------------------*/
	(_abs(Y) <= TNNearlyZero)         /*- Y=0                     -*/
									  /*-                         -*/
		   ||                         /*-                         -*/
									  /*-                         -*/
	(_abs(Dif) < _abs(X * Tol))       /*- Relative change in X    -*/
									  /*-                         -*/
									  /*-                         -*/
 /*       ||                      */  /*-                         -*/
 /*                               */  /*-                         -*/
 /*	(_abs(Dif) < 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 throught 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 */



void ReducePoly(int               *Degree ,
				TNCompVectorPtr   *Poly   ,
				TNcomplex	      Root)

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

{
  int Term;
  TNCompVector *NewPoly = NULL;
  TNcomplex Dummy;

  NewPoly = ( TNCompVector * ) calloc(*Degree,sizeof(TNCompVector));
  *(NewPoly + (*Degree) - 1) = (*((*Poly) + (*Degree)));
  for (Term = (*Degree) - 1; Term >= 1; Term--)
  {
	_Mult(NewPoly + Term, &Root, &Dummy);
	_Add(&Dummy, (*Poly) + Term, NewPoly + Term - 1);
  }
  (*Degree)--;

  free(*Poly);   /***/
  *Poly = NewPoly;
} /* void ReducePoly */



void FindOneRoot(int              Degree     ,
				 TNCompVector	  *Poly       ,
				 TNcomplex	      GuessRoot  ,
				 double	          Tol        ,
				 int	          MaxIter    ,
				 TNcomplex        *Root      ,
				 TNcomplex        *yValue    ,
				 int              *Iter      ,
				 int    *Error)

/*-------------------------------------------------------------------*/
/*- Input:  Degree, Poly, GuessRoot, Tol, MaxIter                   -*/
/*- Output: Root, yValue, Iter, Error                               -*/
/*-                                                                 -*/
/*- This void approximates 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 (yValue), and the number of iterations   -*/
/*- (Iter) are returned. If no root is found, the appropriate error -*/
/*- code (Error) is returned.                                       -*/
/*-------------------------------------------------------------------*/

{
  char Found;
  TNcomplex Dif;
  TNcomplex yPrime, yDoublePrime;

  /* void FindOneRoot */
  *Root = GuessRoot;
  Found = 0;
  *Iter = 0;
  EvaluatePoly(Degree, Poly, *Root, yValue, &yPrime, &yDoublePrime);

  while ( ((*Iter) < MaxIter) && !(Found) )
  {
	(*Iter)++;
	ConstructDifference(Degree, *yValue, yPrime, yDoublePrime, &Dif);

	_Sub(Root, &Dif, Root);
	EvaluatePoly(Degree, Poly, *Root, yValue, &yPrime, &yDoublePrime);
	Found = _TestForRoot(_Modulus(Root), _Modulus(&Dif), _Modulus(yValue), Tol);

  }
  if ( !Found ) *Error = 1;   /* Iterations execeeded MaxIter */
} /* void FindOneRoot */

void Laguerre(int               *Degree    ,
			  TNCompVectorPtr  	*Poly      ,
			  TNcomplex		    InitGuess  ,
			  double		    Tol        ,
			  int		        MaxIter    ,
			  int	            *NumRoots  ,
			  TNCompVectorPtr	*Roots     ,
			  TNCompVectorPtr	*yRoots    ,
			  TNIntVectorPtr	*Iter      ,
			  int     *Error)
{
  int             AddIter;
  int             InitDegree;
  TNCompVector    *InitPoly = NULL;
  TNcomplex       GuessRoot;

  *Roots = ( TNCompVector * ) calloc((*Degree) + 1,sizeof(TNcomplex));
  *yRoots = ( TNCompVector * ) calloc((*Degree) + 1,sizeof(TNcomplex));
  *Iter = ( int * ) calloc((*Degree) + 1,sizeof(TNIntVector));

  InitAndTest(Degree, *Poly, Tol, MaxIter, InitGuess, NumRoots, *Roots,
			  *yRoots, *Iter, &GuessRoot, &InitDegree, &InitPoly, Error);

  while ( ((*Degree) > 0) && !(*Error) )
  {
	FindOneRoot(*Degree, *Poly, GuessRoot, Tol, MaxIter,
				(*Roots) + (*NumRoots) + 1, (*yRoots) + (*NumRoots) + 1,
				(*Iter) + (*NumRoots) + 1, Error);

	if (!(*Error))
	{
	  /*------------------------------------------------------*/
	  /*- 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                -*/
	  /*------------------------------------------------------*/
	  FindOneRoot(InitDegree, InitPoly, *((*Roots) + (*NumRoots) + 1),
				  Tol, MaxIter, (*Roots) + (*NumRoots) + 1,
				  (*yRoots) + (*NumRoots) + 1, &AddIter, Error);
	  *((*Iter) + (*NumRoots) + 1) = *((*Iter) + (*NumRoots) + 1) + AddIter;
	  (*NumRoots)++;
	  ReducePoly(Degree, Poly, *((*Roots) + (*NumRoots))); /* Reduce polynomial */

	}
	GuessRoot = *((*Roots) + (*NumRoots));
  }
  free(InitPoly);
  if (!(*Error)) free(*Poly);
} /* Laguerre */

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

double UserFunction(double X, double (*ProcAddr)(double))
{
 return ProcAddr(X);
}

void  UserProcedure(TNcomplex X, TNcomplex *Y,
					TNcomplex (*ProcAddr)(TNcomplex))
{
 *Y = ProcAddr(X);
}
