/********************************************************************
* Program 1.6: Newton's Second Order Method for the roots of f(x)=0.*
*				Coded by Dilip Datta								*
*																	*
*********************************************************************/

# include "acm.h"

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

	void newton_order2(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's Second Order 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)/100;

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

	void rtfunct(float x);
	void gradient(float x, float *df);

	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																/* Step-11 */
					{
						gradient(x1, &d2f1);
						denom = df1-f1*d2f1/(2*df1);

						if(denom == 0)													/* Step-12 */
						{
							if(fabs(f1) <= eps)	troot(x1, dx, eps, x2, f2, df2, &x1, &f1, &df1);
							else
							{
								x1 = x2;
								f1 = f2;
								df1 = df2;
							}

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

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

								i = nit+1;
							}
							else														/* Step-17 */
							{
								if(x1 < x1p || x1 > x2) x1 = (x1p+x2)/2;				/* Step-18 */
			    				do														/* Step-19 */
								{
									rtfunct(x1); 										/* Step-20 */
									recover_rtfunct(&f1, &df1);

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

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

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

				if(i == nit && fprod < 0) 												/* Step-23 */
				{
					if(fabs(f1) <= 1.0)													/* Step-24 */
					{
// 						printf("\nINSUFFICIENT ITERATION IN FINDING REAL ROOT-%d\n", rcnt+1);
						troot(x1, dx, eps, x2, f2, df2, &x1, &f1, &df1);
					}
					else																/* Step-25 */
					{
// 						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-26 */
					{
						if(fabs(f1)<=eps)												/* Step-27 */
						{
// 							printf("\nINSUFFICIENT ITERATION IN FINDING INFINITE POINT-%d\n", icnt+1);
							troot(x1, dx, eps, x2, f2, df2, &x1, &f1, &df1);
						}
						else															/* Step-28 */
						{
							x1 = x2;
							f1 = f2;
							df1 = df2;
						}				
					}
				}
		
				if(end == TRUE) return;													/* Step-29 */
			}		/* Ending of for(i=1; i<=nit; i++) loop */							/* Step-30 */
		}
        else																			/* Step-31 */
		{
			x1 = x2;
			f1 = f2;
			df1 = df2;
		}
	}

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

void gradient(float x, float *df)
/* Sub-routine to calculate the gradient of a function for given x */
{
	float dx, fx, dfx1, dfx2;

	void rtfunct(float x);

	if(x < 0.01) 	dx = 0.01;			/* dx for finite difference method */
	else			dx = 0.01*x;

	x += dx;							/* Store function values at (x+dx) */
	rtfunct(x);
	recover_rtfunct(&fx, &dfx1);

	x += -2*dx;							/* Store function values at (x-dx) */
	rtfunct(x);
	recover_rtfunct(&fx, &dfx2);

	*df = (dfx1-dfx2)/(2*dx);			/* Derivative of the function */

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