/********************************************************************
* Program 1.4: Secant Method for the roots of f(x)=0.				*
*				Coded by Dilip Datta								*
*																	*
*********************************************************************/

# include "acm.h"

int main()
{
	int i, maxit;
	float dx, eps;

	void secant(int maxit, float dx, float eps);

	printf("\nProblem ID in function.c                                : ");
	scanf("%d", &problem_id);
	printf("\nLower and upper limits of root-searching                : ");			/* Step-1 */
	scanf("%f %f", &xl, &xu);
	printf("\nA small positive number as converging factor            : ");
	scanf("%f", &eps);
	printf("\nNumber of iterations desired                            : ");
	scanf("%d", &maxit);

	printf("\n\n");
	printf("Results from Secant Method for roots of f(x)=0\n");
	printf("----------------------------------------------\n");
	printf("Problem ID                   : %d\n", problem_id);
	printf("Limits of search             : (%f, %f)\n", xl, xu);
	printf("Converging factor            : %f\n", eps);
	printf("Maximum iterations allowed   : %d\n", maxit);

	dx = (xu-xl)/10;

	secant(maxit, dx, eps);

	if(rcnt > 0)
	{
		if(rcnt == 1)	printf("\nRequired real root           :");
		else			printf("\nReal roots (total=%d)         :", rcnt);

		for(i = 1; i <= rcnt-1; i++)  printf(" %f,", rx[i]);
		printf(" %f\n", rx[rcnt]);
	}
	else	printf("\nThere is no real root in the given limit\n");

	if(icnt > 0)
	{
		printf("\nPoints of infinity (total=%d) :", icnt);

		for(i = 1; i <= icnt-1; i++)  printf(" %f,", ix[i]);
		printf(" %f\n", ix[icnt]);
	}

	printf("\nNo. of function evaluations  : %d\n\n", nfun);

	return(0);
}
/********************************************************************/

void secant(int nit, float dx, float eps)
/* Sub-routine for Secant Method for roots of f(x)=0 */
{
	int i;
	float x1, x2, xs, f1, f2, fs, df1, df2, dfs;
	float fprod, dprod;
	float x1o, x2o, f2o, df2o;

	void rtfunct(float x);

	x1 = xl;																		/* Step-2 */

	rtfunct(x1);
	recover_rtfunct(&f1, &df1);

	for(; ;)																		/* Step-3 */
	{
		iroot(&x1, &x2, &f1, &f2, &df1, &df2, dx, eps);								/* Step-4 */

		if(end == TRUE) return;														/* Step-5 */

		fprod = f1*f2;	/* Products of function values & derivatives */				/* Step-6 */
		dprod = df1*df2;
		x1o = x1;		/* Original values of x1, x2, f2 & df2 */
		x2o = x2;
		f2o = f2;
		df2o = df2;
       
		if(fprod < 0 || dprod < 0)													/* Step-7 */
		{
			fvalue = FALSE;  /* Function values at x1 to be copied from at x2 */

			for(i = 1; i <= nit; i++)												/* Step-8 */
			{
				if(fprod < 0 && f1 == f2)											/* Step-9 */
				{
					xs = (x1+x2)/2;
				}
				else																/* Step-10 */
				{
					if(fprod < 0 && f1 != f2)
					{
						xs = (x1*f2-x2*f1)/(f2-f1);
					}
					else															/* Step-11 */
					{
						if(dprod < 0 && df1 == df2)
						{
							xs=(x1+x2)/2;
						}
						else if(dprod < 0 && df1 != df2)	xs = (x1*df2-x2*df1)/(df2-df1);	/* Step-12 */
					}
				}
	
				if(xs < x1o || xs > x2o)  xs = (x1o+x2o)/2;							/* Step-13 */

				rtfunct(xs);
				recover_rtfunct(&fs, &dfs);

				if(fabs(fs) <= eps)													/* Step-14 */
				{ 
					troot(xs, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
					i = nit+1;
				}
				else
				{
					if(fabs(fs) > 1/eps)											/* Step-15 */
					{ 
						root = FALSE;
						troot(xs, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
						i = nit+1;
					}
					else															/* Step-16 */
					{
						x1=x2;
						f1=f2;
						df1=df2;
						x2=xs;
						f2=fs;
						df2=dfs;
					}
				}
	       
				if(i == nit && fprod < 0)											/* Step-17 */
				{
					if(fabs(fs) <= 1.0)												/* Step-18 */
					{
// 						printf("\nINSUFFICIENT ITERATION IN FINDING REAL ROOT-%d\n", rcnt+1);
						troot(xs, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
					}
					else															/* Step-19 */
					{
// 						printf("\nINSUFFICIENT ITERATION IN FINDING INFINITE POINT-%d\n", icnt+1);
						root = FALSE;
						troot(xs, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
					}
				}
				else
				{
					if(i == nit && dprod < 0)										/* Step-20 */
					{
						if(fabs(fs)<=eps)											/* Step-21 */
						{
// 							printf("\nINSUFFICIENT ITERATION IN FINDING INFINITE POINT-%d\n", icnt+1);
							troot(xs, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
						}
						else														/* Step-22 */
						{
							x1=x2o;
							f1=f2o;
							df1=df2o;
						}				
		   			}
				}

				if(end == TRUE) return;												/* Step-23 */
			}		/* Ending of for(i=1; i<=nit; i++) loop */						/* Step-24 */
		}
		else																		/* Step-25 */
		{
			x1=x2;
			f1=f2;
			df1=df2;
		}						
	}
	return;																			/* Step-26 */
}
/********************************************************************/