#include <iostream.h>

double X[21]={0.0};
double Y[21]={0.0};

double Lagrange2(double x,int i,int j,int k)
{
	double l0=(x-X[j])*(x-X[k])/((X[i]-X[j])*(X[i]-X[k]));
	double l1=(x-X[i])*(x-X[k])/((X[j]-X[i])*(X[j]-X[k]));
	double l2=(x-X[i])*(x-X[j])/((X[k]-X[i])*(X[k]-X[j]));
	return Y[i]*l0+Y[j]*l1+Y[k]*l2;
}

double Lagrange(double x, int i, int j)
{
	return Y[i]*(x-X[j])/(X[i]-X[j]) + Y[j]*(x-X[i])/(X[j]-X[i]);
}

double p(double x)
{
	return -50*x/((1+25*x*x)*(1+25*x*x));
}

void main(){

	cout<<endl;
	cout<<"當n＝20 各結點為："<<endl;

	for (int i=0;i<21;i++)
	{
		X[i] = -1+i*(1.0+1.0)/20;
		Y[i]= 1/(1+25*X[i]*X[i]);
		cout<<"x"<<i<<" = "<<X[i]<<"   "<<"y"<<i<<" = "<<Y[i]<<endl;
	}

	cout<<endl;
	cout<<"所求點的真值："<<endl;
	for (i=0;i<21;i++)
	{
		double x=(X[i]+X[i+1])/2;
		double y= 1/(1+25*x*x);
		cout<<"x"<<i<<" = "<<x<<"   "<<"y"<<i<<" = "<<y<<endl;
	}

	cout<<endl;
	cout<<"n＝10, 二次Lagrange 插值法所得結果:"<<endl;
	for( i=0;i+1<11;i++){
		double x=X[2*i+1];
		double y= Lagrange2(x, 0,10,20);
		cout<<"x= "<<x<<"   "<<"y = "<<y<<endl;
		}

	cout<<endl;
	cout<<"n＝20, 二次Lagrange 插值法所得結果:"<<endl;
	for( i=0;i+1<21;i++){
		double x=(X[i]+X[i+1])/2;
		double y= Lagrange2(x, 0,10, 20);
		cout<<"x= "<<x<<"   "<<"y = "<<y<<endl;
	}

	cout<<endl;
	cout<<"n＝10, 分段線性插值法所得結果:"<<endl;
	for( i=0;i+2<21;){
		double x=X[i+1];
		double y= Lagrange(x, i, i+2);
		cout<<"x= "<<x<<"   "<<"y = "<<y<<endl;
		i=i+2;
	}

	cout<<endl;
	cout<<"n＝20, 分段線性插值法所得結果:"<<endl;
	for( i=0;i+1<21;i++){
		double x=(X[i]+X[i+1])/2;
		double y= Lagrange(x, i, i+1);
		cout<<"x= "<<x<<"   "<<"y = "<<y<<endl;
	}

	cout<<endl;
	cout<<"n＝10, I型三次Spline（樣條）插值所得結果:"<<endl;	
	double det=(1.0+1.0)/10;
	double C[11], G[11];
	double A[11]={0.0};
	C[0]=1;
	G[0]=3*(Y[2]-Y[0])/det-p(X[0])*det/2;
	for(i=1;i+1<11;i++)
	{
		A[i]=det/(det+det);
		C[i]=1-A[i];
		G[i]=3*(C[i]*(Y[2*(i+1)]-Y[2*i])/det+A[i]*(Y[2*i]-Y[2*(i-1)])/det);
	}
	A[10]=1;
	G[10]=3*(Y[20]-Y[18])/det+p(Y[20])*det/2;

	double B1[11], B2[11];
	B1[0]=2;
	B2[0]=C[0]/B1[0];
	for(i=1;i<11;i++){
		B1[i]=2-A[i]*B2[i-1];
		B2[i]=C[i]/B1[i];
	}
	G[0]=G[0]/2;
	for(i=1;i<11;i++)
		G[i]=((G[i]-A[i])*G[i-1])/B1[i];
	double M[11];
	M[10]=G[10];
	for(i=9;i>=0;i--)
		M[i]=G[i]-B2[i]*M[i+1];

	double S[11];
	for(i=0;i+1<11;i++)
	{
		S[i]=0;
		double x=X[2*i+1];
		S[i]=S[i]+(det+2*(x-X[2*i]))*(x-X[2*(i+1)])*(x-X[2*(i+1)])*Y[2*i]/(det*det*det);
		S[i]=S[i]+(det-2*(x-X[2*(i+1)]))*(x-X[2*i])*(x-X[2*i])*Y[2*(i+1)]/(det*det*det);
		S[i]=S[i]+(x-X[2*i])*(x-X[2*(i+1)])*(x-X[2*(i+1)])*M[i]/(det*det);
		S[i]=S[i]+(x-X[2*(i+1)])*(x-X[2*i])*(x-X[2*i])*M[i+1]/(det*det);
		cout<<"x= "<<x<<"   "<<"y = "<<S[i]<<endl;
	}
	
	cout<<endl;
	cout<<"n＝20, I型三次Spline（樣條）插值所得結果:"<<endl;	
	det=(1.0+1.0)/20;
	double c[21], g[21];
	double a[21]={0.0};
	c[0]=1;
	g[0]=3*(Y[1]-Y[0])/det-p(X[0])*det/2;
	for(i=1;i+1<21;i++)
	{
		a[i]=det/(det+det);
		c[i]=1-a[i];
		g[i]=3*(c[i]*(Y[i+1]-Y[i])/det+a[i]*(Y[i]-Y[i-1])/det);
	}
	a[20]=1;
	g[20]=3*(Y[20]-Y[19])/det+p(X[20])*det/2;

	double b1[21], b2[21];
	b1[0]=2;
	b2[0]=c[0]/b1[0];
	for(i=1;i<21;i++){
		b1[i]=2-a[i]*b2[i-1];
		b2[i]=c[i]/b1[i];
	}
	g[0]=g[0]/2;
	for(i=1;i<21;i++)
		g[i]=((g[i]-a[i])*g[i-1])/b1[i];
	double m[21];
	m[20]=g[20];
	for(i=19;i>=0;i--)
		m[i]=g[i]-b2[i]*m[i+1];

	double s[21];
	for(i=0;i+1<21;i++)
	{
		s[i]=0;
		double x=(X[i]+X[i+1])/2;
		s[i]=s[i]+(det+2*(x-X[i]))*(x-X[i+1])*(x-X[i+1])*Y[i]/(det*det*det);
		s[i]=s[i]+(det-2*(x-X[i+1]))*(x-X[i])*(x-X[i])*Y[i+1]/(det*det*det);
		s[i]=s[i]+(x-X[i])*(x-X[i+1])*(x-X[i+1])*m[i]/(det*det);
		s[i]=s[i]+(x-X[i+1])*(x-X[i])*(x-X[i])*m[i+1]/(det*det);
		cout<<"x= "<<x<<"   "<<"y = "<<s[i]<<endl;
	}
}