/*
计算方法第五题第二部分
三次样条插值多项式  writed by xlzhao  2001-12-2
*/
# include <stdio.h>
# include <math.h>
# include <string.h>
#pragma hdrstop
void  log_file(char  log_str[256]);
void  view_spl_result(float X[99],float f[99],float S1[99],float S2[99],float S3[99],float S4[99]);
void  max_error(float f[99],float S1[99],float S2[99],float S3[99],float S4[99],float Es[4]);
void  view_error(float Es[4]);
void  view_sql1(float x1[6],float y1[6],float arr_res1[6],float d1[6]);
void  view_sql2(float x2[11],float y2[11],float arr_res2[11],float d2[11]);
void  view_sql3(float x3[16],float y3[16],float arr_res3[16],float d3[16]);
void  view_sql4(float x4[21],float y4[21],float arr_res4[21],float d4[21]);
void  TSS1(float arr_low1[5],float arr_up1[5],float d1[6],float arr_res1[6]);
void  TSS2(float arr_low2[10],float arr_up2[10],float d2[11],float arr_res2[11]);
void  TSS3(float arr_low3[15],float arr_up3[15],float d3[16],float arr_res3[16]);
void  TSS4(float arr_low4[20],float arr_up4[20],float d4[21],float arr_res4[21]);
void	val_spl1(float x1[6],float y1[6],float arr_res1[6],float X[99],float S1[99]);
void	val_spl2(float x2[11],float y2[11],float arr_res2[11],float X[99],float S2[99]);
void	val_spl3(float x3[16],float y3[16],float arr_res3[16],float X[99],float S3[99]);
void	val_spl4(float x4[99],float y4[99],float arr_res4[21],float X[99],float S4[99]);
int   init();
int   n = 5;
FILE  *f1;
char  log_str[8192];
#pragma argsused
int main(int argc, char* argv[])
{
	int i,k;
        char temp[127];
        float x1[6],x2[11],x3[16],x4[21],y1[6],y2[11],y3[16],y4[21],Es[4];
        float d1[6],d2[11],d3[16],d4[21],t1[5],t2[10],t3[15],t4[20];
        float arr_up1[5],arr_up2[10],arr_up3[15],arr_up4[20];
        float arr_low1[5],arr_low2[10],arr_low3[15],arr_low4[20];
        float arr_res1[6],arr_res2[11],arr_res3[16],arr_res4[21];
        float f[99],S1[99],S2[99],S3[99],S4[99],X[99];
        if((init())< 0) return -1 ;
	memset(&log_str,0,8192);
        memset(&temp,0,127);
	/*计算插值点函数值*/
	for(i = 0;i <= 5 ; i ++)
	{
		x1[i] = 0.4*i - 1 ;
		y1[i] = (float) (1/(1 + 25*x1[i]*x1[i]));
	}
	for(i = 0;i <= 10 ; i ++)
	{
		x2[i] = (float) -1 + 0.2*i ;
		y2[i] = (float) 1/(1 + 25*x2[i]*x2[i]);
	}
	for(i = 0;i <= 15 ; i ++)
	{
		x3[i] = (float) -1 + 0.4*i/3 ;
		y3[i] = (float) 1/(1 + 25*x3[i]*x3[i]);
	}
	for(i = 0;i <= 20 ; i ++)
	{
		x4[i] = (float) -1 + 0.1*i ;
		y4[i] = (float) 1/(1 + 25*x4[i]*x4[i]);
	}
	/*end of计算插值点函数值*/
	/*计算newton一阶差商*/
	for(i = 1;i <= 5 ; i ++)
	{
		t1[i-1] = (y1[i]-y1[i-1])/(x1[i]-x1[i-1]) ;
	}
	for(i = 1;i <= 10 ; i ++)
	{
		t2[i-1] = (y2[i]-y2[i-1])/(x2[i]-x2[i-1]) ;
	}
    	for(i = 1;i <= 15 ; i ++)
	{
		t3[i-1] = (y3[i]-y3[i-1])/(x3[i]-x3[i-1]) ;
	}
	for(i = 1;i <= 20 ; i ++)
	{
		t4[i-1] = (y4[i]-y4[i-1])/(x4[i]-x4[i-1]) ;
	}
	/*end of  计算newton一阶差商*/
	/*计算newton二阶差商*6即d 不包括d0,dn*/
	/*确定方程组矩阵的上下边带arr_low[n],arr_up[0]由边界确定*/
	/*由于是等分区间插值，所以上下边带都是0.5*/
	for(i = 0;i <= 3 ; i ++)
	{
		d1[i+1] = 6*(t1[i+1]-t1[i])/(x1[i+2]-x1[i]) ;
		arr_low1[i]= 0.5;
		arr_up1[i+1]= 0.5;
	}
	for(i = 0;i <= 8 ; i ++)
	{
		d2[i+1] = 6*(t2[i+1]-t2[i])/(x2[i+2]-x2[i]) ;
		arr_low2[i]= 0.5;
		arr_up2[i+1]= 0.5;
	}
    	for(i = 0;i <= 13 ; i ++)
	{
		d3[i+1] = 6*(t3[i+1]-t3[i])/(x3[i+2]-x3[i]) ;
		arr_low3[i]= 0.5;
		arr_up3[i+1]= 0.5;
	}
	for(i = 0;i <= 18 ; i ++)
	{
		d4[i+1] = 6*(t4[i+1]-t4[i])/(x4[i+2]-x4[i]) ;
		arr_low4[i]= 0.5;
		arr_up4[i+1]= 0.5;
	}
	/*end of  计算newton二阶差商*/
	/*边界条件利用已知函数求边界的导数即书上的的二种情况
	arr_up[0] = 1;
	arr_low[n] = 1;
	f(x)=1/(1+25*x*x),f'(x)=-50x/((1+25*x*x)(1+25*x*x))
	f'(-1) = 50/676,f'(1) = -50/676;
	d0,dn,arr_low[n],arr_up[0]*/
	arr_low1[4] = 1.0 ;
	arr_low2[9] = 1.0 ;
	arr_low3[14] = 1.0 ;
	arr_low4[19] = 1.0 ;
	arr_up1[0] = 1.0 ;
	arr_up2[0] = 1.0 ;
	arr_up3[0] = 1.0 ;
	arr_up4[0] = 1.0 ;
	d1[0]=6.0/(2.0/5)*(t1[0]-50.0/676);
	d2[0]=6.0/(2.0/10)*(t2[0]-50.0/676);
	d3[0]=6.0/(2.0/15)*(t3[0]-50.0/676);
	d4[0]=6.0/(2.0/20)*(t4[0]-50.0/676);
	d4[5]=6.0/(2.0/5)*(-50/676-t4[4]);
	d4[10]=6.0/(2.0/10)*(-50/676-t4[9]);
	d4[15]=6.0/(2.0/15)*(-50/676-t4[14]);
	d4[15]=6.0/(2.0/20)*(-50/676-t4[19]);
	/*利用TSS算法解方程组得到样条插值函数的插值点二阶导数值向量arr_res*/
	TSS1(arr_low1,arr_up1,d1, arr_res1);
	TSS2(arr_low2,arr_up2,d2, arr_res2);
	TSS3(arr_low3,arr_up3,d3, arr_res3);
	TSS4(arr_low4,arr_up4,d4, arr_res4);
	/*计算函数值*/
	for(k = 1 ; k < 100; k ++)
	{
		X[k-1] = 2.0/100.0*k - 1;
		f[k-1] = 1.0/(1.0 + 25.0*X[k-1]*X[k-1]);
	}/*end of计算函数值*/
	view_sql1(x1,y1,arr_res1,d1);
	view_sql2(x2,y2,arr_res2,d2);
	view_sql3(x3,y3,arr_res3,d3);
	view_sql4(x4,y4,arr_res4,d4);
	val_spl1(x1,y1,arr_res1,X,S1);
	val_spl2(x2,y2,arr_res2,X,S2);
	val_spl3(x3,y3,arr_res3,X,S3);
	val_spl4(x4,y4,arr_res4,X,S4);
	view_spl_result(X,f,S1,S2,S3,S4);
	max_error(f,S1,S2,S3,S4,Es);
	view_error(Es);
	fclose(f1);
        return 0;
}
int  init()
{
    if((f1 = fopen("reslt52.txt","w"))==NULL)
       {
        printf("can not create result text file");
        return -1;
        }
    fputs("*****三次样条插值多项式*****\n",f1);
    return 0;
}
void  log_file(char  log_str[8192])
{
    log_str[8191]='\0';
    printf("%s",log_str);
    fputs(log_str,f1);
    log_str[0]='\0';
    return;
}

void view_spl_result(float X[99],float f[99],float S1[99],float S2[99],float S3[99],float S4[99])
{
	int  i;
	char str_tmp[255];
	log_file("\nx         f(x)       S1(x)       S2(x)       S3(x)       S4(x)\n");
	for(i = 0;i < 99 ; i ++)
	{
		memset(&str_tmp,0,sizeof(str_tmp));
		sprintf(str_tmp,"%f ",X[i]);
		log_file(str_tmp);
		sprintf(str_tmp,"%f ",f[i]);
		log_file(str_tmp);
		sprintf(str_tmp,"%f ",S1[i]);
		log_file(str_tmp);
		sprintf(str_tmp,"%f ",S2[i]);
		log_file(str_tmp);
		sprintf(str_tmp,"%f ",S3[i]);
		log_file(str_tmp);
		sprintf(str_tmp,"%f\n",S4[i]);
		log_file(str_tmp);
	}
	return;
}
void  max_error(float f[99],float S1[99],float S2[99],float S3[99],float S4[99],float Es[4])
{
	int i,j;
	float f_bew[4];
	Es[0] = fabs(f[0]-S1[0]);
	Es[1] = fabs(f[0]-S2[0]);
	Es[2] = fabs(f[0]-S3[0]);
	Es[3] = fabs(f[0]-S4[0]);
	for(i = 1 ; i < 99 ; i ++)
	{
		f_bew[0] = f[i]-S1[i];
		f_bew[1] = f[i]-S2[i];
		f_bew[2] = f[i]-S3[i];
		f_bew[3] = f[i]-S4[i];
		for (j = 0 ; j < 4 ; j ++)
		if(fabs(f_bew[j])>fabs(Es[j])) Es[j] = fabs(f_bew[j]);
	}
	return ;
}

void  view_error(float Es[4])
{
	int i;
	char str_tmp[255];
	log_file("\n****三次样条插值函数误差：****\n");
	for(i = 0;i < 4 ; i ++)
	{
		memset(&str_tmp,0,sizeof(str_tmp));
		sprintf(str_tmp,"第%d次三次样条最大误差：%f \n",i,Es[i]);
		log_file(str_tmp);
	}
	return;
}
void  TSS1(float arr_low1[5],float arr_up1[5],float d1[6],float arr_res1[6])
{
	int k;
        int arr_diag[6] = {2,2,2,2,2,2};
        float arr_q[6] ;
        float arr_y[6];
        float arr_p[5];
        arr_q[0]= arr_diag[0];
        arr_y[0]= d1[0];
        for(k = 2;k <= 6;k ++)
        {
        arr_p[k-2]= arr_low1[k-2]/arr_q[k-2];
        arr_q[k-1]= arr_diag[k-1]-arr_p[k-2]*arr_up1[k-2];
        arr_y[k-1]= d1[k-1]-arr_p[k-2]*arr_y[k-2];
        }
        arr_res1[5]= arr_y[5]/arr_q[5];
        for(k = 6 - 2 ;k >= 0;k --)
        {
        arr_res1[k]=(arr_y[k]-arr_up1[k]*arr_res1[k+1])/arr_q[k];
        }
        return;
}
void  TSS2(float arr_low2[10],float arr_up2[10],float d2[11],float arr_res2[11])
{
	int k;
        int arr_diag[11] = {2,2,2,2,2,2,2,2,2,2,2};
        float arr_q[11] ;
        float arr_y[11];
        float arr_p[10];
        arr_q[0]= arr_diag[0];
        arr_y[0]= d2[0];
        for(k = 2;k <= 11;k ++)
        {
        arr_p[k-2]= arr_low2[k-2]/arr_q[k-2];
        arr_q[k-1]= arr_diag[k-1]-arr_p[k-2]*arr_up2[k-2];
        arr_y[k-1]= d2[k-1]-arr_p[k-2]*arr_y[k-2];
        }
        arr_res2[10]= arr_y[10]/arr_q[10];
        for(k = 11 - 2 ;k >= 0;k --)
        {
        arr_res2[k]=(arr_y[k]-arr_up2[k]*arr_res2[k+1])/arr_q[k];
        }
        return;
}
void  TSS3(float arr_low3[15],float arr_up3[15],float d3[16],float arr_res3[16])
{
	int k;
        int arr_diag[16] = {2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2};
        float arr_q[16] ;
        float arr_y[16];
        float arr_p[15];
        arr_q[0]= arr_diag[0];
        arr_y[0]= d3[0];
        for(k = 2;k <= 16;k ++)
        {
        arr_p[k-2]= arr_low3[k-2]/arr_q[k-2];
        arr_q[k-1]= arr_diag[k-1]-arr_p[k-2]*arr_up3[k-2];
        arr_y[k-1]= d3[k-1]-arr_p[k-2]*arr_y[k-2];
        }
        arr_res3[15]= arr_y[15]/arr_q[15];
        for(k = 16 - 2 ;k >= 0;k --)
        {
        arr_res3[k]=(arr_y[k]-arr_up3[k]*arr_res3[k+1])/arr_q[k];
        }
        return;
}
void  TSS4(float arr_low4[20],float arr_up4[20],float d4[21],float arr_res4[21])
{
	int k;
        int arr_diag[21] = {2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2};
        float arr_q[21] ;
        float arr_y[21];
        float arr_p[20];
        arr_q[0]= arr_diag[0];
        arr_y[0]= d4[0];
        for(k = 2;k <= 21;k ++)
        {
        arr_p[k-2]= arr_low4[k-2]/arr_q[k-2];
        arr_q[k-1]= arr_diag[k-1]-arr_p[k-2]*arr_up4[k-2];
        arr_y[k-1]= d4[k-1]-arr_p[k-2]*arr_y[k-2];
        }
        arr_res4[20]= arr_y[20]/arr_q[20];
        for(k = 21 - 2 ;k >= 0;k --)
        {
        arr_res4[k]=(arr_y[k]-arr_up4[k]*arr_res4[k+1])/arr_q[k];
        }
        return;
}
void	val_spl1(float x1[6],float y1[6],float arr_res1[6],float X[99],float S1[99])
{
	int i,k;
	float h = 2.0/5;
	for(i = 0;i < 99 ; i ++)
	{
		for(k = 0 ;k < 6 ; k ++)
		{
		if((X[i]<=x1[k])&&(X[i]>=x1[k-1])) break;
		}
		S1[i]=(arr_res1[k-1]*(x1[k]-X[i])*(x1[k]-X[i])*(x1[k]-X[i])/6
		     + arr_res1[k]*(X[i]-x1[k-1])*(X[i]-x1[k-1])*(X[i]-x1[k-1])/6
		     + (y1[k-1]-arr_res1[k-1]*h*h/6)*(x1[k]-X[i])
		     + (y1[k]-arr_res1[k]*h*h/6)*(X[i]-x1[k-1]))/h;
	}
	return;
}
void	val_spl2(float x2[11],float y2[11],float arr_res2[11],float X[99],float S2[99])
{
	int i,k;
	float h = 2.0/10;
	for(i = 0;i < 99 ; i ++)
	{
		for(k = 0 ;k < 11 ; k ++)
		{
		if((X[i]<=x2[k])&&(X[i]>=x2[k-1])) break;
		}
		S2[i]=(arr_res2[k-1]*(x2[k]-X[i])*(x2[k]-X[i])*(x2[k]-X[i])/6
		     + arr_res2[k]*(X[i]-x2[k-1])*(X[i]-x2[k-1])*(X[i]-x2[k-1])/6
		     + (y2[k-1]-arr_res2[k-1]*h*h/6)*(x2[k]-X[i])
		     + (y2[k]-arr_res2[k]*h*h/6)*(X[i]-x2[k-1]))/h;
	}
	return;
}
void	val_spl3(float x3[16],float y3[16],float arr_res3[16],float X[99],float S3[99])
{
	int i,k;
	float h = 2.0/15;
	for(i = 0;i < 99 ; i ++)
	{
		for(k = 0 ;k < 16 ; k ++)
		{
		if((X[i]<=x3[k])&&(X[i]>=x3[k-1])) break;
		}
		S3[i]=(arr_res3[k-1]*(x3[k]-X[i])*(x3[k]-X[i])*(x3[k]-X[i])/6
		     + arr_res3[k]*(X[i]-x3[k-1])*(X[i]-x3[k-1])*(X[i]-x3[k-1])/6
		     + (y3[k-1]-arr_res3[k-1]*h*h/6)*(x3[k]-X[i])
		     + (y3[k]-arr_res3[k]*h*h/6)*(X[i]-x3[k-1]))/h;
	}
	return;
}
void	val_spl4(float x4[99],float y4[99],float arr_res4[21],float X[99],float S4[99])
{
	int i,k;
	float h = 2.0/20;
	for(i = 0;i < 99 ; i ++)
	{
		for(k = 0 ;k < 21 ; k ++)
		{
		if((X[i]<=x4[k])&&(X[i]>=x4[k-1])) break;
		}
		S4[i]=(arr_res4[k-1]*(x4[k]-X[i])*(x4[k]-X[i])*(x4[k]-X[i])/6
		     + arr_res4[k]*(X[i]-x4[k-1])*(X[i]-x4[k-1])*(X[i]-x4[k-1])/6
		     + (y4[k-1]-arr_res4[k-1]*h*h/6)*(x4[k]-X[i])
		     + (y4[k]-arr_res4[k]*h*h/6)*(X[i]-x4[k-1]))/h;
	}
	return;
}

void  view_sql1(float x1[6],float y1[6],float arr_res1[6],float d1[6])
{
	int i;
	char tmp[32];
	log_file("\n*****当n为5时*****\n自变量的值:\n");
	for(i = 0;i < 6 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",x1[i]);
		log_file(tmp);
	}
	log_file("\n函数值:\n");
	for(i = 0;i < 6 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",y1[i]);
		log_file(tmp);
	}
	log_file("\n对应插值方程的d:\n");
	for(i = 0;i < 6 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",d1[i]);
		log_file(tmp);
	}
	log_file("\n样条函数对应点的二阶导数值:\n");
	for(i = 0;i < 6 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",arr_res1[i]);
		log_file(tmp);
	}
}
void  view_sql2(float x2[11],float y2[11],float arr_res2[11],float d2[11])
{
	int i;
	char tmp[32];
	log_file("\n*****当n为10时*****\n自变量的值:\n");
	for(i = 0;i < 11 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",x2[i]);
		log_file(tmp);
	}
	log_file("\n函数值:\n");
	for(i = 0;i < 11 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",y2[i]);
		log_file(tmp);
	}
	log_file("\n对应插值方程的d:\n");
	for(i = 0;i < 11 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",d2[i]);
		log_file(tmp);
	}
	log_file("\n样条函数对应点的二阶导数值:\n");
	for(i = 0;i < 11 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",arr_res2[i]);
		log_file(tmp);
	}
}
void  view_sql3(float x3[16],float y3[16],float arr_res3[16],float d3[16])
{
	int i;
	char tmp[32];
	log_file("\n*****当n为15时*****\n自变量的值:\n");
	for(i = 0;i < 16 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",x3[i]);
		log_file(tmp);
	}
	log_file("\n函数值:\n");
	for(i = 0;i < 16 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",y3[i]);
		log_file(tmp);
	}
	log_file("\n对应插值方程的d:\n");
	for(i = 0;i < 16 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",d3[i]);
		log_file(tmp);
	}
	log_file("\n样条函数对应点的二阶导数值:\n");
	for(i = 0;i < 16 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",arr_res3[i]);
		log_file(tmp);
	}
}
void  view_sql4(float x4[21],float y4[21],float arr_res4[21],float d4[21])
{
	int i;
	char tmp[32];
	log_file("\n*****当n为20时*****\n自变量的值:\n");
	for(i = 0;i < 21 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",x4[i]);
		log_file(tmp);
	}
	log_file("\n函数值:\n");
	for(i = 0;i < 21 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",y4[i]);
		log_file(tmp);
	}
	log_file("\n对应插值方程的d:\n");
	for(i = 0;i < 21 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",d4[i]);
		log_file(tmp);
	}
	log_file("\n样条函数对应点的二阶导数值:\n");
	for(i = 0;i < 21 ;i ++)
	{
		memset(&tmp,0,sizeof(tmp));
		sprintf(tmp,"%f  ",arr_res4[i]);
		log_file(tmp);
	}
}