/********************************************************************
* Program 1.5: Newton-Raphson Method for the roots of f(x)=0.		*
*				Coded by Dilip Datta								*
*																	*
*********************************************************************/

# include "acm.h"

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

	void newt_raph(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 Newton Raphson 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)/1000;

	newt_raph(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 newt_raph(int nit, float dx, float eps)
/* Sub-routine for Newton Raphson Method for roots of f(x)=0 */
{
	int i;
	float x1, x2, f1, f2, df1, df2;
	float x1p, fprod, dprod, pvalue1, pvalue2, dvalue1, dvalue2;

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

		if((fprod = f1*f2) < 0)	pvalue1 = -1; /* Products of function values & */	/* Step-6 */
		else                    pvalue1 = 1;  /* derivatives and their natures */

		if((dprod = df1*df2) < 0)	dvalue1 = -1;
		else						dvalue1 = 1;

		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(fabs(f1) > 1/eps)												/* Step-9 */
				{ 
					root = FALSE;
					troot(x1, dx, eps, x2, f2, df2, &x1, &f1, &df1);
					i = nit+1;
				}
				else
				{
					if(f1 == 0)														/* Step-10 */
					{ 
						troot(x1, dx, eps, x2, f2, df2, &x1, &f1, &df1);
						i = nit+1;
					}
					else
					{
						if(df1 == 0)												/* Step-11 */
						{
							if(fabs(f1) <= eps) troot(x1, dx, eps, x2, f2, df2, &x1, &f1, &df1);/* Step-12 */
							else													/* Step-13 */
							{
								x1=x2;
								f1=f2;
								df1=df2;
							}				

							i = nit+1;
						}
						else														/* Step-14 */
						{
							x1p = x1;			/* Previous value of x1 */
							x1 -= f1/df1;

							if(fabs(x1-x1p) <= eps)									/* Step-15 */
							{
								if(fabs(f1) <= eps) troot(x1, dx, eps, x2, f2, df2, &x1, &f1, &df1);/* Step-16 */
								else												/* Step-17 */
								{
									if(fprod < 0)									/* Step-18 */
									{
										root = FALSE;
										troot(x1, dx, eps, x2, f2, df2, &x1, &f1, &df1);
									}
									else											/* Step-19 */
									{
										x1 = x2;
										f1 = f2;
										df1 = df2;
									}			
								}

								i = nit+1;
							}
							else													/* Step-20 */
							{ 
								if(x1 < x1p || x1 > x2) x1 = (x1p+x2)/2;			/* Step-21 */

								do													/* Step-22 */
								{
									rtfunct(x1); 									/* Step-23 */
									recover_rtfunct(&f1, &df1);

									if(f1*f2 < 0)	pvalue2 = -1;					/* Step-24 */
									else			pvalue2 = 1;

									if(df1*df2 < 0)	dvalue2 = -1;
									else			dvalue2 = 1;

									if(pvalue1 != pvalue2 && dvalue1 != dvalue2) x1 = (x1p+x1)/2; /* Step-25 */
								}while(pvalue1 != pvalue2 && dvalue1 != dvalue2);
							}		
						}
					}
				}

				if(i == nit && fprod < 0) 											/* Step-26 */
				{
					if(fabs(f1) <= 1.0)												/* Step-27 */
					{
// 						printf("\nINSUFFICIENT ITERATION IN FINDING REAL ROOT-%d\n", rcnt+1);
						troot(x1, dx, eps, x2, f2, df2, &x1, &f1, &df1);
					}
					else															/* Step-28 */
					{
// 						printf("\nINSUFFICIENT ITERATION IN FINDING INFINITE POINT-%d\n", icnt+1);
						root = FALSE;
						troot(x1, dx, eps, x2, f2, df2, &x1, &f1, &df1);
					}
				}
				else
				{
					if(i == nit && dprod < 0)										/* Step-29 */
					{
						if(fabs(f1) <= eps)											/* Step-30 */
						{
// 							printf("\nINSUFFICIENT ITERATION IN FINDING INFINITE POINT-%d\n", icnt+1);
							troot(x1, dx, eps, x2, f2, df2, &x1, &f1, &df1);
						}
						else														/* Step-31 */
						{
							x1 = x2;
							f1 = f2;
							df1 = df2;
						}
					}
				}

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