/*
	Title: Newton divided difference interpolating polynomial,
	given in Section 5.2.

	This routine will do interpolation of  fcn(x), defined below.
	The program reads the degree  n  and an initial point  x(0).
	Then the divided differences
         	f[x(0),x(1),...,x(j)] ,  1 <= j <= n
	are constructed for use in evaluating the interpolating
	polynomial.  The interpolation nodes x(j) are evenly spaced.
	
	Following the above, an x at which the interpolation is to be
	done is requested from the user.  For this x, the
	interpolating polynomials of degree <= n are
	constructed, along with the error in the interpolation.
	The interpolation is carried out with the newton divided
	difference form of the interpolating polynomial.
*/

#include <stdio.h>
#include <math.h>

#define MAX_N	20
#define fcn(x)	exp(x)

void divdif(int n,float x[],float f[]) ;


main()
{
	char q;
	float x[MAX_N+1], f[MAX_N+1], df[MAX_N+1];
	FILE *fp;
	int i, j, n;
	float error, h, poly, z;

	/*  Input problem parameters  */

	printf("\n Do you want the output on the terminal ?  (y/n)");
	printf("\n If not, then it is stored on file interp.output : ");
	fflush(stdin);
	q=getchar();
	if ((q == 'N') || (q == 'n'))  
	{
		if ((fp = fopen("interp.output", "w")) == NULL)  
		{
			printf("\n Cannot open file interp.output\n");
			return(1);
		}
	}
	else
		fp = stdout;
      
	printf("\n Degree, x(0), h = ? : ");
	scanf("%d %f %f", &n, &x[0], &h);

	/*  Set up node and function values  */
      
	for (i=0; i <= n; i++)  {
		x[i] = x[0] + i*h;
		f[i] = fcn(x[i]);
		df[i] = f[i];
	}

	/*  Calculate divided differences and print them  */

	divdif(n,x,df);
	fprintf(fp, "\n  i    x[i]      f[i]        df[i]\n");
	for (i=0; i <= n; i++)
		fprintf(fp, "\n%3d    %3.1f    %8.6f    %14.7e", i, x[i], f[i], df[i]);

	/*  Begin interpolation loop  */

	while (1)  
	{
		printf("\n\n Give interpolation point z : ");
		scanf("%f", &z);
		fprintf(fp, "\n\n\n x = %9.7f\n", z);

		/*  Construct interpolation polynomial at z, for
			all degrees i <= n  */
      
		for(i=1; i <= n; i++)  
		{
			poly = df[i];
			for (j=i-1; j >= 0; j--)
				poly = (z - x[j])*poly + df[j];
			error = fcn(z) - poly;
			fprintf(fp, "\n Degree = %2d    Interpolant = %10.7f", i, poly);
			fprintf(fp, "    Error = %12.3e", error);
		}

		printf("\n\n Another interpolation point x? (y/n) : ");
		scanf("\n%c", &q);
		if ((q == 'N') || (q == 'n')) 
		{
			if (fp != stdout)
				fclose(fp);
			break;
		}
	}
}
      


void divdif(int n,float x[],float f[]) 
{

/*
	Input: (1) 'n' denotes a positive integer.
             	(2) 'x' denotes a vector of points x(0),...,x(n).
             	(3) 'f' denotes a vector of values of some function
                 	evaluated at the nodes in x.
      	Output: The vector f is converted to a vector of divided
                 	differences: f(i)=f[x(0),...,x(i)], i=0,...,n     
*/

	int i, j;

	for (i=1; i <= n; i++)
		for (j=n; j >= i; j--)
			f[j] = (f[j] - f[j-1])/(x[j] - x[j-i]);

}