// *Hermite's Polynominal Interpolation with Differentials* by Krzysztof Bosak, 1999-01-07

#include <stdio.h>
#include <memory.h>

FILE *plik=fopen("out.txt", "wt");

void plot(int x, int y)
{
	static int ox=65536, oy=65536;
	if(ox!=x || oy!=y)
	{
		fprintf(plik, "%d = X\n%d =Y\n", x, y);
		ox=x;
		oy=y;
	}
}

void cross(int x, int y)
{
	for(int j=-4; j<=4; j++)
	{
		plot(x+j, y+j);
		plot(x-j, y+j);
	}
}

int FindDegree(double *xarray, double terminator)
{
	int d=0;
	while(xarray[d] != terminator)
		d++;
	return d-1;
}

void PlotNodes(int n, double xscale, double yscale, double *xarray, double *yarray)
{
	for(int i=0; i<=n; i++)
		cross(static_cast<int>(xarray[i]*xscale), static_cast<int>(yarray[i]*yscale));
}

void HermiteInterpolation(int deg, double *xarray, double *yarray, double *result)
{	
	memcpy(result, yarray, (deg+1)*sizeof(double));
	double val;
	for(int k=1; k<=deg; k++)
	{		
		val=result[deg];
		for(int l=deg; l>=k; l--)
		{
			if(xarray[l] != xarray[l-k])
			{
				result[l]=(val-result[l-1])/(xarray[l]-xarray[l-k]);
				val=result[l-1];
			}
			else 
				result[l]=result[l-1]/k;
		}
	}
};

double HermiteValue(int deg, double x, double *xarray, double *yarray)
{
	yarray+=deg;
	deg--;
	xarray+=deg;
	double result=*yarray;
	for(; deg>=0; deg--)
	{		
		yarray--;
		result*=x-*xarray;
		xarray--;
		result+=*yarray;
	}
	return result;
};

int main()
{
	const double EndOfList=1E99;				//terminująca, zastrzeżona wartość x
	double xtab[]={0, 0, 0, 40, 100, EndOfList};		//dane OX
	double ytab[]={0, 0, 0, 40, 0};					//dane OY i pochodne
	double result[256];
	const double xscale=1;						//skalowanie OX
	const double yscale=1;						//skalowanie OY
	const double x1=-50, x2=150, step=0.01;		//początek, koniec rysowania i krok
	plot(250, 70);	//set axes mode
	plot(0, 0);		//x: connect points
	int N=FindDegree(xtab, EndOfList);
	PlotNodes(N, xscale, yscale, xtab, ytab);

	HermiteInterpolation(N, xtab, ytab, result);
	for(double x=x1; x<=x2; x+=step)
		plot(static_cast<int>(x*xscale), static_cast<int>(HermiteValue(N, x, xtab, result)*yscale));
	return 0;
}