//xi		0.5			1.0			1.5			2.0			2.5			3.0
//f(xi)	0.541688	-0.042060	0.375089	0.624706	-0.887597	0.099809
//輸入節點數及節點數據，計算并輸出差商表；
//輸入待使用的插值多項式次數，待求的自變量 x，輸出內插需要選擇的節點及 Newton 插值多項式；
//用 3 次 Newton 插值多項式計算并輸出 f (1.25)；
//用 5 次 Newton 插值多項式計算并輸出 f (1.987)。
//
//

#include <stdio.h>
#include <malloc.h>
#include <conio.h>
#include <stdlib.h>
typedef struct vex
{ 
	double x;
	double y;
}vex,*pVex;
void function(pVex dvex,double F[5][5],int style,double x,double vex[],char* Newton,double result);
void main()
{
	int N=6;
	pVex dvex=(pVex)malloc(N*sizeof(vex));
	dvex[0].x=0.5;dvex[0].y=0.541688;
	dvex[1].x=1.0;dvex[1].y=-0.042060;
	dvex[2].x=1.5;dvex[2].y=0.375089;
	dvex[3].x=2.0;dvex[3].y=0.624706;
	dvex[4].x=2.5;dvex[4].y=-0.887597;
	dvex[5].x=3.0;dvex[5].y=0.099809;
	//////////////////////////
	int i,j;
	double F[5][5];
	///////////////////////
	//差商表的形成過程
	for(i=0;i<N-1;i++)
		F[i][0]=(dvex[i+1].y-dvex[i].y)/(dvex[i+1].x-dvex[i].x);
	for(j=1;j<N-1;j++)
		for(i=j;i<N-1;i++)
			F[i][j]=(F[i][j-1]-F[i-1][j-1])/(dvex[i+1].x-dvex[i-j].x);
	/////////////////////////
	printf("下面輸出Newton差商表:\n");
	printf("x值\ty值\t\t一階差商\t二階差商\t三階差商\t四階差商\t五階差商\n");
	printf("%.1f\t%f\n",dvex[0].x,dvex[0].y);
	printf("%.1f\t%f\t%f\n",dvex[1].x,dvex[1].y,F[0][0]);
	printf("%.1f\t%f\t%f\t%f\n",dvex[2].x,dvex[2].y,F[1][0],F[1][1]);
	printf("%.1f\t%f\t%f\t%f\t%f\n",dvex[3].x,dvex[3].y,F[2][0],F[2][1],F[2][2]);
	printf("%.1f\t%f\t%f\t%f\t%f\t%f\n",dvex[4].x,dvex[4].y,F[3][0],F[3][1],F[3][2],F[3][3]);
	printf("%.1f\t%f\t%f\t%f\t%f\t%f\t%f\n",dvex[5].x,dvex[5].y,F[4][0],F[4][1],F[4][2],F[4][3],F[4][4]);
	/////////////////////////
	double pVex[5]={0,0,0,0,0};
	char pNewton[200]="";
	double result=0;
	////////////////////////////////
	printf("------------------------------------------------------------\n");
	printf("下面用3次Newton插值多項式計算并輸出f(1.125)\n");
	function(dvex,F,3,1.125,pVex,pNewton,result);
	printf("------------------------------------------------------------\n");
	printf("下面用5次Newton插值多項式計算并輸出f(1.987)");
	function(dvex,F,5,1.987,pVex,pNewton,result);
	////////////////////////////////
	printf("------------------------------------------------------------\n");
	printf("下面，請輸入待使用的插值多項式次數，待求的自變量 x\n");
	int nN;float dX;
	for(;;)
	{
		printf("\n待使用的插值多項式次數(輸入0退出程序):");
		scanf("%d",&nN);
		if(nN==0)exit(0);
		while(nN<1||nN>5)
		{
			printf("輸入的次數建議在1~5之間\n");
			scanf("%d",&nN);
		}
		printf("待求的自變量 x:");
		scanf("%f",&dX);
		while(dX<dvex[0].x||dX>dvex[5].x)
		{
			printf("輸入的自變量不在%f和%f之間\n",dvex[0].x,dvex[5].x);
			scanf("%f",&dX);
		}
		printf("下面開始顯示結果\n");

		function(dvex,F,nN,dX,pVex,pNewton,result);
	}
	getch();	
}
void function(pVex dvex,double F[5][5],int style,double x,double vex[],char* Newton,double result)
{
	int i=0;
	int vex_i;
	int j=0;
	switch(style)
	{
	case 1:
		{
			while(x>dvex[i].x)i++;
			printf("\n所選節點為:\n");
			for(vex_i=i-1,j=0;vex_i<i+1;vex_i++,j++)printf("%1.1f\t",vex[j]=dvex[vex_i].x);
			sprintf(
					Newton,
					"f(x)=\n\t\t%f+\n\t\t%f*(x-%1.1f)",
					dvex[i-1].y,
					F[i-1][0],dvex[i-1].x
					);
			printf("\n表達式為:%s",Newton);
			result=dvex[i-1].y+F[i-1][0]*(x-dvex[i-1].x);
			printf("\n其f(%.6f)結果為%2.6f\n",x,result);
			break;
		}
	case 2:
		{
			while(x>dvex[i].x)i++;
			if(i==5)
			{
				printf("\n所選節點為:\n%1.1f\t%1.1f\t%1.1f\n",dvex[3].x,dvex[4].x,dvex[5].x);
				sprintf(
					Newton,
					"f(x)=\n\t\t%f+\n\t\t%f*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)",
					dvex[3].y,
					F[3][0],dvex[3].x,
					F[4][1],dvex[3].x,dvex[4].x
					);
				printf("\n表達式為:%s",Newton);
				result=dvex[3].y+F[3][0]*(x-dvex[3].x)+F[4][1]*(x-dvex[4].x);
			}
			else
			{
				printf("\n所選節點為:\n");
				for(vex_i=i-1,j=0;vex_i<i+2;vex_i++,j++)printf("%1.1f\t",vex[j]=dvex[vex_i].x);
				sprintf(
					Newton,
					"f(x)=\n\t\t%f+\n\t\t%f*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)",
					dvex[i-1].y,
					F[i-1][0],dvex[i-1].x,
					F[i][1],dvex[i-1].x,dvex[i].x
					);
				printf("\n表達式為:%s",Newton);
				result=dvex[i-1].y+F[i-1][0]*(x-dvex[i-1].x)+F[i][1]*(x-dvex[i-1].x)*(x-dvex[i].x);
			}
			printf("\n其f(%.6f)結果為%2.6f\n",x,result);
			break;
		}
	case 3:
		{
			while(x>dvex[i].x)i++;
			//x[0]	y[0]										F[i][j]
			//x[1]	y[1]	F[0][0]
			//x[2]	y[2]	F[1][0]		F[1][1]
			//x[3]	y[3]	F[2][0]		F[2][1]		F[2][2]
			//x[4]	y[4]	F[3][0]		F[3][1]		F[3][2]		F[3][3]
			//x[5]	y[5]	F[4][0]		F[4][1]		F[4][2]		F[4][3]		F[4][4]
			if(i==5||x==4)
			{
				printf("\n所選節點為:\n%1.1f\t%1.1f\t%1.1f\t%1.1f\n",dvex[2].x,dvex[3].x,dvex[4].x,dvex[5].x);
				sprintf(
					Newton,
					"f(x)=\n\t\t%f+\n\t\t%f*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)",
					dvex[2].y,
					F[2][0],dvex[2].x,
					F[3][1],dvex[2].x,dvex[3].x,
					F[4][2],dvex[2].x,dvex[3].x,dvex[4].x
					);
				printf("\n表達式為:%s",Newton);
				result=dvex[2].y+F[2][0]*(x-dvex[2].x)+F[3][1]*(x-dvex[2].x)*(x-dvex[3].x)+F[4][2]*(x-dvex[2].x)*(x-dvex[3].x)*(x-dvex[4].x);
			}
			else
			{
				printf("\n所選節點為:\n");
				for(vex_i=i-1,j=0;vex_i<i+3;vex_i++,j++)printf("%1.1f\t",vex[j]=dvex[vex_i].x);
				sprintf(
					Newton,
					"f(x)=\n\t\t%f+\n\t\t%f*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)",
					dvex[i-1].y,
					F[i-1][0],dvex[i-1].x,
					F[i][1],dvex[i-1].x,dvex[i].x,
					F[i+1][2],dvex[i-1].x,dvex[i].x,dvex[i+1].x
					);
				printf("\n表達式為:%s",Newton);
				result=dvex[i-1].y+F[i-1][0]*(x-dvex[i-1].x)+F[i][1]*(x-dvex[i-1].x)*(x-dvex[i].x)+F[i+1][2]*(x-dvex[i-1].x)*(x-dvex[i].x)*(x-dvex[i+1].x);
			}
			
			printf("\n其f(%.6f)結果為%2.6f\n",x,result);
			break;
		}
	case 4:
		{
			while(x>dvex[i].x)i++;
			if(i==5||i==4||x==3)
			{
				printf("\n所選節點為:\n%1.1f\t%1.1f\t%1.1f\t%1.1f\t%1.1f\n",dvex[1].x,dvex[2].x,dvex[3].x,dvex[4].x,dvex[5].x);
				sprintf(
					Newton,
					"f(x)=\n\t\t%f+\n\t\t%f*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)",
					dvex[1].y,
					F[1][0],dvex[1].x,
					F[2][1],dvex[1].x,dvex[2].x,
					F[3][2],dvex[1].x,dvex[2].x,dvex[3].x,
					F[4][3],dvex[1].x,dvex[3].x,dvex[3].x,dvex[4].x
					);
				printf("\n表達式為:%s",Newton);
				result=dvex[2].y+F[2][0]*(x-dvex[2].x)+F[3][1]*(x-dvex[2].x)*(x-dvex[3].x)+F[4][1]*(x-dvex[2].x)*(x-dvex[3].x)*(x-dvex[4].x);
			}
			else
			{
				printf("\n所選節點為:\n");
				for(vex_i=i-1;vex_i<i+4;vex_i++,j++)printf("%1.1f\t",vex[j]=dvex[vex_i].x);
				sprintf(
					Newton,
					"f(x)=\n\t\t%f+\n\t\t%f*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)",
					dvex[i-1].y,
					F[i-1][0],	dvex[i-1].x,
					F[i][1],	dvex[i-1].x,	dvex[i].x,
					F[i+1][2],	dvex[i-1].x,	dvex[i].x,	dvex[i+1].x,
					F[i+2][3],	dvex[i-1].x,	dvex[i].x,	dvex[i+1].x,	dvex[i+2].x
					);
				printf("\n表達式為:%s",Newton);
				result=dvex[i-1].y+F[i-1][0]*(x-dvex[i-1].x)+F[i][1]*(x-dvex[i-1].x)*(x-dvex[i].x)+
					F[i+1][2]*(x-dvex[i-1].x)*(x-dvex[i].x)*(x-dvex[i+1].x)+F[i+2][3]*(x-dvex[i-1].x)*(x-dvex[i].x)*(x-dvex[i+1].x)*(x-dvex[i+2].x);
			}
			printf("\n其f(%.6f)結果為%2.6f\n",x,result);
			break;
		}
	case 5:
		{
			printf("\n所選節點為:\n");
			for(j=0;j<6;j++)printf("%1.1f\t",dvex[j].x);
			for(vex_i=0;vex_i<6;vex_i++)vex[i]=dvex[i].x;
			sprintf(
				Newton,
				"f(x)=\n\t\t%f+\n\t\t%f*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)+\n\t\t%f*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)*(x-%1.1f)",
				dvex[0].y,
				F[0][0],dvex[0].x,
				F[1][1],dvex[0].x,dvex[1].x,
				F[2][2],dvex[0].x,dvex[1].x,dvex[2].x,
				F[3][3],dvex[0].x,dvex[1].x,dvex[2].x,dvex[3].x,
				F[4][4],dvex[0].x,dvex[1].x,dvex[2].x,dvex[3].x,dvex[4].x
				);
			printf("\n表達式為:%s",Newton);
			result=dvex[0].y+
				F[0][0]*(x-dvex[0].x)+
				F[1][1]*(x-dvex[0].x)*(x-dvex[1].x)+
				F[2][2]*(x-dvex[0].x)*(x-dvex[1].x)*(x-dvex[2].x)+
				F[3][3]*(x-dvex[0].x)*(x-dvex[1].x)*(x-dvex[2].x)*(x-dvex[3].x),
				F[4][4]*(x-dvex[0].x)*(x-dvex[1].x)*(x-dvex[2].x)*(x-dvex[3].x)*(x-dvex[4].x);
			printf("\n其f(%.6f)結果為%2.6f\n",x,result);
			break;
		}
	}
}