//#include<stdio.h>
#include<iostream>
//#include<cstdlib>
#include<cmath>
///#include<string>
//#include<ctime>
using namespace std;
#define pi 3.1415926
//#define float double  //----這個當N>90時，浮點行運算溢出，無法得到結果。所以將浮點轉成雙精度。



float Fun_x(float x){
	float y;
	y=1/(1+x*x);
	return y;
}//F(x)--原函數


float Lx_1(float x,int N)
{
	float *x1=new float[N+1],*y1=new float[N+1],LX1=0.0,tmp;
	float n;
	int i,j;
	for(i=0;i<=N;i++)
	{
		n=N;
		i=float (i);
		x1[i]=-5+10*i/n;
		i=int (i);
		y1[i]=Fun_x(x1[i]);
	}
	for(i=0;i<=N;i++)
	{
		tmp=1.0;
		for(j=0;j<i;j++)
			tmp=tmp*(x-x1[j])/(x1[i]-x1[j]);

		for(j=i+1;j<=N;j++)
				tmp=tmp*(x-x1[j])/(x1[i]-x1[j]);
		LX1=LX1+tmp*y1[i];
	}
	return LX1;
}//計算第一種插值算法求得的LX;

float Lx_2(float x,int N)
{
	float n;
	int i,j;
	float *x2=new float[N+1],*y2=new float[N+1],LX2=0.0,tmp,k;
	for(i=0;i<=N;i++)
	{
		i=float (i);
		n=N;
		k=(2*i+1)/(2*n+2);
		i=int (i);
		x2[i]=(-5)*cos(k*pi);
		y2[i]=Fun_x(x2[i]);
	}
	for(i=0;i<=N;i++)
	{
		tmp=1.0;
		for(j=0;j<i;j++)
			tmp=tmp*(x-x2[j])/(x2[i]-x2[j]);
		for(j=i+1;j<=N;j++)
			tmp=tmp*(x-x2[j])/(x2[i]-x2[j]);
		LX2=LX2+tmp*y2[i];
	}
	return LX2;
}//計算第二種插值算法求得的LX；


void main(){
	float X[101],error1,error2,temp1,temp2;//X[]用于存放測試誤差的數yi；error1/2分別是LX1和LX2的最大誤差。
	int i,j,N ,n[4]={5,10,20,40};
	FILE *fs;
	fs=fopen("D:\\lab02.txt","w");
	for(j=0;j<=3;j++)
	{
		N=n[j];
		error1=0.0,error2=0.0,temp1=0.0,temp2=0.0;
		for(i=70;i<=100;i++)
		{
			i=float (i);//用于計算X[i]，故將類型轉換。
			X[i]=-5+0.1*i;
			i=int (i);
			temp1=fabs(Fun_x(X[i])-Lx_1(X[i],N));
			if(error1<temp1)error1=temp1;//求最大的誤差
			temp2=fabs(Fun_x(X[i])-Lx_2(X[i],N));
			if(error2<temp2)error2=temp2;//同上
		}
		cout<<"N="<<N<<endl;
		cout<<"Max Error of grid (1) : "<<error1<<endl;
		cout<<"Max Error of grid (2) : "<<error2<<endl;//輸出
		fprintf(fs,"N= %d \n",N);
		fprintf(fs,"Max Error of grid (1) : %2.8f\n",error1);
		fprintf(fs,"Max Error of grid (2) : %2.8f\n",error2);//將結果輸入到文檔
	}

}