/********************************************************************
*	acm.h: Header file to be used in programs given in the book		*
*			"Applied Computational Methods in C - Dilip Datta"		*
*						Coded by Dilip Datta						*
*																	*
*********************************************************************/

# include <stdio.h>
# include <math.h>
# include <malloc.h>
/********************************************************************/

int problem_id;			/* Problem ID in function.c */
# include "function.c"
/********************************************************************/

# define TRUE 1
# define FALSE 0
/********************************************************************/

int max_nrt = 101;		/* Maximum number of roots */
int rcnt = 0, icnt = 0;	/* Counters for roots & points of infinity */
int root = TRUE;		/* Real root = TRUE (default), infinite point = FALSE */
int fvalue = TRUE;		/* Function values to be calculated = TRUE (default) */
int nfun = 0;			/* Number of function evaluations */
int end = FALSE;		/* TRUE if the program is to be terminated */

float xl, xu;			/* Range of search for the roots */
float fx, dfx;			/* Values of function and its derivative at given point x */
float *rx, *ix;			/* Pointers to store roots & points of infinity 
			   			   or real & imaginary parts of roots           */
/********************************************************************/

void iroot(float *a, float *b, float *fa, float *fb, float *dfa, float *dfb, float dx, float eps);
void qdroot(float a, float b, float c, float real[3], float imag[3]);
void recover_rtfunct(float *fv, float *dv);
void troot(float xt, float dx, float eps, float x, float f2, float df2, float *x1, float *f1, float *df1);
/********************************************************************/

void iroot(float *a, float *b, float *fa, float *fb, float *dfa, float *dfb, float dx, float eps)
/* Function to check the roots for initial guesses */ 
{
	float x1, x2, f1, f2, df1, df2;

	void rtfunct(float x);
	void recover_rtfunct(float *fv, float *dv);
  
	x1 = *a;
	x2 = *b;
	f1 = *fa;
	f2 = *fb;
	df1 = *dfa;
	df2 = *dfb;

	do																		/* Step-1 */
	{
		do																	/* Step-2 */
		{
			if(f1 == 0)	troot(x1, dx, eps, x1, f1, df1, &x1, &f1, &df1);	/* Step-3 */
			else															/* Step-4 */
			{
				if(fabs(f1) > 1/eps)
				{
					root = FALSE;
					troot(x1, dx, eps, x1, f1, df1, &x1, &f1, &df1);
				}
			}

			if(end == TRUE)	return; 										/* Step-5 */
		}while(f1 == 0 || fabs(f1) >= 1/eps);								/* Step-6 */

		if(x1 >= xu)														/* Step-7 */
		{
			end = TRUE;
			return;
		}

		if((x2 = x1+dx) > xu)  x2 = xu;										/* Step-8 */
		rtfunct(x2);														/* Step-9 */
		recover_rtfunct(&f2, &df2);

		if(f2 == 0)	troot(x2, dx, eps, x2, f2, df2, &x1, &f1, &df1);		/* Step-10 */
		else																/* Step-11 */
		{
			if(fabs(f2) > 1/eps)
			{ 
				root = FALSE;
				troot(x2, dx, eps, x2, f2, df2, &x1, &f1, &df1);
			}
		}

		if(end == TRUE) return; 											/* Step-12 */
	}while(f1 == 0 || fabs(f1) > 1/eps);									/* Step-13 */

	*a = x1;
	*b = x2;
	*fa = f1;
	*fb = f2;
	*dfa = df1;
	*dfb = df2;

	return;
}
/********************************************************************/

void qdroot(float a, float b, float c, float real[3], float imag[3])
/* Sub-routine for quadratic roots of ax^2 + bx + c = 0 */
{
	float disc;

	disc = b*b - 4*a*c;														/* Step-1 */
	if(disc >= 0)															/* Step-2 */
	{
		real[1] = (-b+sqrt(disc))/(2*a);
		real[2] = (-b-sqrt(disc))/(2*a);
		imag[1] = 0;
		imag[2] = 0;
	}
	else																	/* Step-3 */
	{
		real[1] = -b/(2*a);
		real[2] = real[1];
		imag[1] = sqrt(-disc)/(2*a);
		imag[2] = -imag[1];
	}

	return;
}
/********************************************************************/

void recover_rtfunct(float *fv, float *dv)
/* Assign the values of function and its derivative (fx & dfx) to fv & dv 
   and also store the number of function evaluation */
{
	*fv = fx;
	*dv = dfx;
	nfun += 1;

	return;		
}
/********************************************************************/

void troot(float xt, float dx, float eps, float x2, float f2, float df2, float *x1, float *f1, float *df1)
/* Function to store the roots of trascendental equation and to
   compute new point and function value for next root search    */
{
	float x, fv, dfv;

	void rtfunct(float x);
	void recover_rtfunct(float *fv, float *dv);

	if(xt > xu)					/* xt is outside of search range */			/* Step-1 */
	{
		end = TRUE;
		return;
	} 

	if(root == TRUE)			/* Root found */							/* Step-2 */
	{
		if(rcnt == 0) 		rx = (float *)malloc(max_nrt*sizeof(float));
		if(rcnt >= max_nrt)	rx = (float *)realloc(rx, sizeof(float));

		rx[++rcnt] = xt;
	}
	else						/* Function approaches infinity */			/* Step-3 */
	{
		if(icnt == 0) 		ix = (float *)malloc(max_nrt*sizeof(float));
		if(icnt >= max_nrt)	ix = (float *)realloc(ix, sizeof(float));

        ix[++icnt] = xt;
	}

	if(fvalue == FALSE)			/* Values copied */    						/* Step-4 */
	{
		x = x2;
		fv = f2;
		dfv = df2;
	}	
	else						/* Function values to be calculated */		/* Step-5 */
	{
		if((x = xt+dx) > xu)												/* Step-6 */
		{
			end = TRUE;
			return;
		}

		rtfunct(x);															/* Step-7 */
		recover_rtfunct(&fv, &dfv);
    
		while((root == TRUE && fv == 0) || (root == FALSE && fabs(fv) > 1/eps))	/* Step-8 */
		{
			if((x += dx) >= xu)												/* Step-9 */
			{
				end=TRUE;
				return;
			}

			rtfunct(x);														/* Step-10 */ 
			recover_rtfunct(&fv, &dfv); 
		}
	}

	root = TRUE;			/* Default value of root */						/* Step-11 */
	fvalue = TRUE;			/* Default value of fvalue */

	*x1 = x;				/* Return value of x1 and its */
	*f1 = fv;				/* corresponding function & */
	*df1 = dfv;				/* derivative values */

	return;
}
/********************************************************************/