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

#define PI 3.1415926536
#define TNNearlyZero  1E-20

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 UserFunction(double X,double (*FuncAddr)(double));
void  UserProcedure(TNcomplex X, TNcomplex *Y,
					TNcomplex (*ProcAddr)(TNcomplex));
double _abs(double x);

void Muller(TNcomplex      Guess   ,
			double	       Tol      ,
			int	           MaxIter  ,
			TNcomplex      *Answer  ,
			TNcomplex      *yAnswer ,
			int            *Iter    ,
			int  *Error   ,
			TNcomplex (*ProcPtr)(TNcomplex));

/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-               Input: Guess, Tol, MaxIter                                 -*/
/*-              Output: Answer, yAnswer, Iter, Error                        -*/
/*-                                                                          -*/
/*-             Purpose: This program uses Muller's method to find a root    -*/
/*-                      of a user defined function Y=TNTargetF given an     -*/
/*-                      initial approximation.  The root may be complex.    -*/
/*-                                                                          -*/
/*-                                                                          -*/
/*-   User-Defined                                                           -*/
/*-          Procedures: TNTargetF(X : TNcomplex; VAR Y : TNcomplex);        -*/
/*-                                                                          -*/
/*-  User-Defined Types: TNcomplex = record                                  -*/
/*-                                    Re, Im : real;                        -*/
/*-                                  end;                                    -*/
/*-                                                                          -*/
/*-    Global Variables: Guess   : real;            initial guess            -*/
/*-                      Tol     : real;            tolerance in the         -*/
/*-                                                 answer                   -*/
/*-                      MaxIter : integer;         maximum number of        -*/
/*-                                                 iterations               -*/
/*-                      Answer  : TNcomplex;       a root of the            -*/
/*-                                                 polynomial               -*/
/*-                      yAnswer : TNcomplex;       value of the             -*/
/*-                                                 polynomial at the        -*/
/*-                                                 root (close to zero)     -*/
/*-                      Iter    : integer;         number of iterations     -*/
/*-                                                 it took to find root     -*/
/*-                      Error   : byte;            flags an error           -*/
/*-                                                                          -*/
/*-              Errors: 0: No errors                                        -*/
/*-                      1: Iter > MaxIter                                   -*/
/*-                      2: parabola could not                               -*/
/*-                         be formed                                        -*/
/*-                      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 SquareRoot(TNcomplex C1, TNcomplex *C2)
{
  double R, Theta;

  R = sqrt(pow(C1.Re,2) + pow(C1.Im,2));
  if ( _abs(C1.Re) < TNNearlyZero )
	{
	  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 SquareRoot */


void MakeParabola(TNcomplex    X0     ,
				  TNcomplex	   X1     ,
				  TNcomplex	   Approx ,
				  TNquadratic  *Factor ,
				  int  *Error,
				  TNcomplex (*ProcPtr)(TNcomplex))

/*--------------------------------------------------------*/
/*- Input:  X0, X1, Approx                               -*/
/*- Output: Factor, Error                                -*/
/*-                                                      -*/
/*- This void constructs a parabola to fit the      -*/
/*- three points X0, X1, Approx.  The intersection of    -*/
/*- this parabola with the x-axis will yield the next    -*/
/*- approximation.  If the parabola is a horizontal line -*/
/*- ) Error = 2 since a horizontal line will !      -*/
/*- intersect the x-axis.                                -*/
/*--------------------------------------------------------*/

{
  TNcomplex  H1, H2, H3, H;
  TNcomplex  Delta1, Delta2;
  TNcomplex  Dum1, Dum2, Dum3, Dum4;     /* Dummy variables */

  Sub(X0, Approx, &H1);
  Sub(X1, Approx, &H2);
  Sub(X0, X1, &H3);
  Mult(H2, H3, &Dum1);
  Mult(H1, Dum1, &H);
  if ( Modulus(&H) < TNNearlyZero )
	*Error = 2;                     /* Can't fit a quadratic to these points */
  UserProcedure(X1, &Dum1, ProcPtr);
  Sub(Dum1, Factor->C, &Delta1);             /* C was passed in */
  UserProcedure(X0, &Dum1, ProcPtr);
  Sub(Dum1, Factor->C, &Delta2);

  if (!(*Error))                 /* Calculate coefficients of quadratic */
  {
	Mult(H1, H1, &Dum1);             /* Calculate B */
	Mult(Dum1, Delta1, &Dum2);
	Mult(H2, H2, &Dum1);
	Mult(Dum1, Delta2, &Dum3);
	Sub(Dum2, Dum3, &Dum4);
	Divide(Dum4, H, &Factor->B);

	Mult(H2, Delta2, &Dum1);         /* Calculate A */
	Mult(H1, Delta1, &Dum2);
	Sub(Dum1, Dum2, &Dum3);
	Divide(Dum3, H, &Factor->A);
  }

  if ( (Modulus(&Factor->A) <= TNNearlyZero) &&
	   (Modulus(&Factor->B) <= TNNearlyZero) )
	*Error = 2;            /* Test if ( parabola is actually a constant */
} /* void MakeParabola */

void Initial(TNcomplex  Guess      ,
			 double	    Tol        ,
			 int	    MaxIter    ,
			 int *Error   ,
			 int        *Iter       ,
			 char       *Found      ,
			 TNcomplex  *NewApprox  ,
			 TNcomplex  *yNewApprox ,
			 TNcomplex  *X0         ,
			 TNcomplex  *X1         ,
			 TNcomplex  *OldApprox  ,
			 TNquadratic  *Factor,
			 TNcomplex (*ProcPtr)(TNcomplex))

/*------------------------------------------------------------------*/
/*- Input:  Guess, Tol, MaxIter                                    -*/
/*- Output: Error, Iter, Found, NewApprox, yNewApprox,             -*/
/*-         X0, X1, OldApprox, Factor                              -*/
/*-                                                                -*/
/*- This function initializes all the above variables. It          -*/
/*- sets OldApprox equal to Guess, X0 and X1 are set close to      -*/
/*- Guess.  The void also checks the tolerance (Tol) and      -*/
/*- maximum number of iterations (MaxIter) for errors.             -*/
/*------------------------------------------------------------------*/

{
  TNcomplex  cZero, yOldApprox;


  cZero.Re = 0;    /* Complex zero */
  cZero.Im = 0;    /* Complex zero */
  *Error = 0;
  *Found = 0;
  *Iter = 0;
  *NewApprox = cZero;
  *yNewApprox = cZero;
		  /* X0 and X1 are points which are close to Guess */
  X0->Re = -0.75 * Guess.Re - 1;     X0->Im = Guess.Im;
  X1->Re = 0.75 * Guess.Re;          X1->Im = 1.2 * Guess.Im + 1;
  *OldApprox = Guess;
  UserProcedure(*OldApprox, &yOldApprox, ProcPtr);  /* Evaluate the function at OldApprox */
  Factor->A = cZero;
  Factor->B = cZero;
  Factor->C = yOldApprox;
  MakeParabola(*X0, *X1, *OldApprox, Factor, Error, ProcPtr);
  if ( Tol <= 0 )
	*Error = 3;
  if ( MaxIter < 0 )
	*Error = 4;
} /* void Initial */

void Muller_QuadraticFormula(TNquadratic  Factor,
					  TNcomplex	   OldApprox,
					  TNcomplex    *NewApprox)

/*-------------------------------------------------------------*/
/*- Input:  Factor, OldApprox                                 -*/
/*- Output: NewApprox                                         -*/
/*-                                                           -*/
/*- This void applies the complex quadratic formula      -*/
/*- to the quadratic Factor to determine where the parabola   -*/
/*- represented by Factor intersects the x-axis. The solution -*/
/*- of the quadratic formula is subtracted from OldApprox to  -*/
/*- yield NewApprox.                                          -*/
/*-------------------------------------------------------------*/

{
  TNcomplex  Discrim, Difference, Dum1, Dum2, Dum3;

  Mult(Factor.B, Factor.B, &Dum1);                     /* B^2 */
  Mult(Factor.A, Factor.C, &Dum2);
  Discrim.Re = Dum1.Re - 4 * Dum2.Re;  /* B^2 - 4AC */
  Discrim.Im = Dum1.Im - 4 * Dum2.Im;
  SquareRoot(Discrim, &Discrim);
  Sub(Factor.B, Discrim, &Dum1);                /* B +/- sqrt(B^2 - 4AC) */
  Add(Factor.B, Discrim, &Dum2);
  /* Choose the root with B +/- Discrim greatest */
  if ( Modulus(&Dum1) < Modulus(&Dum2) )
	Dum1 = Dum2;
  Add(Factor.C, Factor.C, &Dum3);
  if ( Modulus(&Dum1) < TNNearlyZero )  /* if ( B +/- sqrt(B^2 - 4AC) = 0 */
	{
	  NewApprox->Re = 0;
	  NewApprox->Im = 0;
	}
  else
	{                               /* 2C/[B +/- sqrt(B^2 - 4AC)] */
	  Divide(Dum3, Dum1, &Difference);
	  /* Calculate NewApprox */
	  Sub(OldApprox, Difference, NewApprox);
	}
} /* void Muller_QuadraticFormula */

int Muller_TestForRoot(TNcomplex  X    ,
				TNcomplex  OldX ,
				TNcomplex  Y    ,
				double	   Tol)

{
  TNcomplex Dif, FracDif;

/*-----------------------------------------------------------------*/
/*- 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.                      -*/
/*-----------------------------------------------------------------*/


  Sub(X, OldX, &Dif);
  FracDif.Re = X.Re * Tol;
  FracDif.Im = X.Im * Tol;

  return                              /*---------------------------*/
	((Modulus(&Y) <= TNNearlyZero))      /*- Y = 0                   -*/
									  /*-                         -*/
		   ||                         /*-                         -*/
									  /*-                         -*/
	(Modulus(&Dif) < Modulus(&FracDif)) /*- Relative change in X    -*/
									  /*-                         -*/
									  /*-                         -*/
/*         ||                      */ /*-                         -*/
/*                                 */ /*-                         -*/
/*  (Modulus(&Dif) < Tol)           */ /*- Absolute change in X    -*/
/*                                 */ /*-                         -*/
/*         ||                      */ /*-                         -*/
/*                                 */ /*-                         -*/
/*  (Modulus(&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.                                 -*/
/*-----------------------------------------------------------------------*/

} /* function Muller_TestForRoot */


void Muller(TNcomplex       Guess   ,
			double	       Tol      ,
			int	           MaxIter  ,
			TNcomplex      *Answer  ,
			TNcomplex      *yAnswer ,
			int            *Iter    ,
			int  *Error   ,
			TNcomplex (*ProcPtr)(TNcomplex))
{



  TNcomplex X0, X1, OldApprox,
  NewApprox, yNewApprox;                /* Iteration variables */
  TNquadratic Factor  ;                 /* Factor of polynomial */
  char Found ;                          /* Flags that a factor */
										/* has been found */


  /* function  Muller */
  Initial(Guess, Tol, MaxIter, Error, Iter, &Found, &NewApprox, &yNewApprox,
		  &X0, &X1, &OldApprox, &Factor, ProcPtr);

  while ( !Found && !(*Error) && (*Iter < MaxIter) )
  {
	(*Iter)++;
	Muller_QuadraticFormula(Factor, OldApprox, &NewApprox);
	UserProcedure(NewApprox, &yNewApprox, ProcPtr);    /* Calculate a new yNewApprox */
	Found = Muller_TestForRoot(NewApprox, OldApprox, yNewApprox, Tol);
	X0 = X1;
	X1 = OldApprox;
	OldApprox = NewApprox;
	Factor.C = yNewApprox;
	MakeParabola(X0, X1, OldApprox, &Factor, Error, ProcPtr);
  }
  *Answer = NewApprox;
  *yAnswer = yNewApprox;
  if ( Found )
	*Error = 0;
  else
	if ( ((*Error) == 0) && ((*Iter) >= MaxIter) )
	  *Error = 1;
} /* function Muller */

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);
}
