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

# include "acm.h"

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

	void bisect(int n, 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 bisections desired                            : ");
	scanf("%d", &n);

	printf("\n\n");
	printf("Results from Bisection 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("Number of bisections desired : %d\n", n);
	printf("Converging factor            : %f\n", eps);

	dx = (xu-xl)/10;

	bisect(n, 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 bisect(int nb, float dx, float eps)
/* Sub-routine for Bisection Method for roots of f(x)=0 */
{
	int i;
	float x1, x2, xm, f1, f2, fm, df1, df2, dfm;
	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 <= nb; i++)											/* Step-8 */
			{
				xm = (x1+x2)/2;

				rtfunct(xm);
				recover_rtfunct(&fm, &dfm);

				if(fm == 0)														/* Step-9 */
				{ 
					troot(xm, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);

					i = nb+1;
				}
				else															/* Step-10 */
				{
					if(fabs(fm) > 1/eps)										/* Step-11 */
					{ 
						root=FALSE;
						troot(xm, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
						i = nb+1;
					}
					else														/* Step-12 */
					{
						if(f1*fm < 0 || df1*dfm < 0)							/* Step-13 */
						{
							x2 = xm;
							f2 = fm;
							df2 = dfm;
						}		
						else
						{
							x1 = xm;
							f1 = fm;
							df1 = dfm;
						}
					}
				}
	       
				if(i == nb && fprod < 0) 										/* Step-14 */
				{
					if(fabs(fm) <= 1.0)											/* Step-15 */
					{
// 						printf("\nINSUFFICIENT BISECTIONS IN FINDING REAL ROOT-%d\n", rcnt+1);
						troot(xm, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
					}
					else														/* Step-16 */
					{
// 						printf("\nINSUFFICIENT BISECTIONS IN FINDING INFINITE POINT-%d\n", icnt+1);
						root=FALSE;
						troot(xm, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
					}
				}
				else															/* Step-17 */
				{
					if(i == nb && dprod < 0)									/* Step-18 */
					{
						if(fabs(fm) <= eps)										/* Step-19 */
						{
// 							printf("\nINSUFFICIENT BISECTIONS IN FINDING INFINITE POINT-%d\n", icnt+1);
							troot(xm, dx, eps, x2o, f2o, df2o, &x1, &f1, &df1);
						}
						else													/* Step-20 */
						{
							x1 = x2o;
							f1 = f2o;
							df1 = df2o;
						}
					}
				}

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

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