//C3
//Interpolation,using Newton Methom

#include <iostream.h>

const int N=5;

//Divied Difference
void fd(double *ax,double *ay,int n,double f[N][N])
{
	int i,j;
	for (i=0;i<n;i++)
		f[i][0]=ay[i];
	for (i=1;i<n;i++)
	{
		for (j=1;j<=i;j++)
		{
			f[i][j]=(f[i][j-1]-f[i-1][j-1])/(ax[i]-ax[i-j]);
		}
	}
}
	   
//Newton polynomial
double NP(double x,double f[N][N],double *ax,int n)
{
	double sum=0,Tsum=1;
	int i,j;
	for (i=0;i<n;i++)
	{
		for (j=0;j<i;j++)
		{
			Tsum=Tsum*(x-ax[j]);
		}
		sum=sum+Tsum*f[i][i];
		Tsum=1;
	}
	return sum;
}

void main()
{
	double ax[]={0.4,0.55,0.65,0.8,0.9},ay[]={0.41075,0.57815,0.69675,0.88811,1.02652},x1=0.596,x2=0.895;
	double func[N][N]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
	int i;
	fd(ax,ay,N,func);
	cout<<"The value at x1 is: "<<NP(x1,func,ax,N)<<endl;
	cout<<"The value at x2 is: "<<NP(x2,func,ax,N)<<endl;
	cin>>i;
}
//F(x1)=0.631918,F(x2)=1.01937.Newton法運(yùn)算時(shí)間與Lagrange基本一致
//但可變更性遠(yuǎn)強(qiáng)于前者