// INTERPOLATION
// Program 9 : To impliment Stirling's interpolation formula
// ------------------------------------------------------------------------
//
// Table: I use only one 2D array to store all values of x, y and delta y's
// Logic ( i and j are the loop variables )
//	
//	  i\j |   0  |   1  |   2   |   3     |   4  |
//	  ----+------+------+-------+---------+------+...
//	   0  |  x1  |  y1  | y2-y1 | dy2-dy1 |      |
//	  ----+------+------+-------+---------+------+...
//	   1  |  x2  |  y2  | y3-y2 |         |      |
//	  ----+------+------+-------+---------+------+...
//	   2  |  x3  |  y3  |       |         |      !
//	  ----+------+------+-------+---------+...
//	      !      !      !       !         !      

#include <stdio.h>
#include <conio.h>

#define MAXITEMS 10

main()
{
	int i,j,i1,i2;	// for loop
	int n;		// no of terms
	// 2D array. 1st column for x, 2nd for y and other for difference
	double table[MAXITEMS][MAXITEMS+1];
	double h;			// constant difference
	double x, y = 0.0;	// y is f(x)
	double u;			// m is multiplier
	double m1 = 1.0, m2 = 1.0; 	// 1 for forward, 2 for backward
	int z;				// index of f(0)
	
	printf(" Enter no of terms : ");
	scanf("%d", &n);
	printf(" Enter first term : ");
	scanf("%lg", &table[0][0]);
	printf(" Enter difference between terms : ");
	scanf("%lg", &h);
	
	for( i=0; i<n; i++)
	{
		// calculate other x from first term and difference
		table[i][0] = table[0][0] + (double)(i) * h;
		printf(" Enter the value of f(%lg) : ", table[i][0]);
		scanf("%lg", &table[i][1]);
	}
	
	for( j=1; j<n; j++)	// loop to each column
		for( i=0; i<n-j; i++)	// loop to cell within a column
			table[i][j+1] = table[i+1][j] - table[i][j];
	
	printf("\n Enter the value of x for which you want to find f(x) : ");
	scanf("%lg", &x);
	
	//Finding the value of u for best case
	for( i=0; i<n; i++)
	{
		u = ( x - table[i][0] ) / h;
		if( -.25<=u && u<=0.25 )
			break;
	}
	if( i == n ) //search unsuccessful for best case
	{
		//Finding the value of u for general case
		for( i=0; i<n; i++)
		{
			u = ( x - table[i][0] ) / h;
			if( -.5<=u && u<=0.5 )
				break;
		}
		if( i == n ) //search unsuccessful
		{
			printf("\n Condition -.5<=u<=0.5 not satisfied");
			getch();
			return 1;//exit program
		}
		else
			z = i;
	}
	else
		z = i;
	
	printf("\n u = %g is found at f(%g)", u, table[z][0]);
	
	for( j=1, i1=0, i2=0; j<=n; j++)
	{
		// stirling is the mean of gauss forward & backward
		y += ( m1 * table[z-i1][j] + m2 * table[z-i2][j] ) / 2.0;
		
		if( j%2 == 1 )
		{
			m1 *= (u+double(i1)) / (double)(j);
			m2 *= (u-double(i2++)) / (double)(j);
		}
		else
		{
			m1 *= (u-double(++i1)) / (double)(j);
			m2 *= (u+double(i2)) / (double)(j);
		}
		
		if( (z-i1)<0 || (z-i1)>(n-j-1) || (z-i2)<0 || (z-i2)>(n-j-1) )
			break;	//checks validity of index in the table
	}
	
	printf("\n F(%g) = %g", x, y);	
	getch();
}