/********************************************************************
* Program 1.3: Regula-Falsi Method for the roots of f(x)=0.			*
*				Coded by Dilip Datta								*
*																	*
*********************************************************************/

# include "acm.h"

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

	void rfalsi(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 Regula Falsi 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;

	rfalsi(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 rfalsi(int nit, float dx, float eps)
/* Sub-routine for Regula Falsi Method for roots of f(x)=0 */
{
	int i;
	float x1, x2, xf, f1, f2, ff, df1, df2, dff;
	float fprod, dprod;
	float 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;
		x2o = x2;		/* Original values of x2, f2 & df2 */
		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(f1*f2 < 0) 	xf = (x1*f2-x2*f1)/(f2-f1);						/* Step-9 */
				else			xf = (x1*df2-x2*df1)/(df2-df1);

				rtfunct(xf);													/* Step-10 */
				recover_rtfunct(&ff, &dff);

				if(ff == 0)														/* Step-11 */
				{ 
					troot(xf, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
					i = nit+1;
				}
				else															/* Step-12 */
				{
					if(fabs(ff) > 1/eps)										/* Step-13 */
					{ 
						root = FALSE;
						troot(xf, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
						i = nit+1;
					}
					else														/* Step-14 */
					{
						if(f1*ff < 0 || df1*dff < 0)							/* Step-15 */
						{
							x2=xf;
							f2=ff;
							df2=dff;
						}		
						else
						{
							x1=xf;
							f1=ff;
							df1=dff;
						}
					}
				}
	       
	        	if(i == nit && fprod < 0) 										/* Step-16 */
				{
					if(fabs(ff) <= 1.0)											/* Step-17 */
					{
// 						printf("\nINSUFFICIENT ITERATION IN FINDING REAL ROOT-%d\n", rcnt+1);
						troot(xf, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
					}
					else														/* Step-18 */
					{
// 						printf("\nINSUFFICIENT ITERATION IN FINDING INFINITE POINT-%d\n", icnt+1);
						root=FALSE;
						troot(xf, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
					}
				}
				else															/* Step-19 */
				{
					if(i == nit && dprod < 0)									/* Step-20 */
					{
						if(fabs(ff) <= eps)										/* Step-21 */
						{
// 							printf("\nINSUFFICIENT ITERATION IN FINDING INFINITE POINT-%d\n", icnt+1);
							troot(xf, 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;
		}
	}																			/* Step-26 */

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