#include   <stdio.h>
#include   <stdlib.h>  
#include   <malloc.h>  
#include   <math.h>   
void Smooth(double   *x,double   *y,double   *a,int   n,int   m,double   *dt1,double   *dt2,double   *dt3   );


void   main()   

{   
	int   i,n,m;   
	double   *x,*y,*a,dt1,dt2,dt3,b;   
	
	n = 20;   
	m = 6;   
	b = 0;   
	
	/*分別為x,y,a分配存貯空間*/   
	
	x = (double*)calloc(n,sizeof(double));   

	y = (double*)calloc(n,sizeof(double));   

	a = (double*)calloc(n,sizeof(double));   

	for(i = 1; i <= n; i++)   
	{   
		x[i - 1] = b + (i - 1) * 0.1;   
		/*每隔0.1取一個點，這樣連續取n個點*/   
		y[i - 1] = x[i - 1] - exp( -x[i-1]);   
		//y[i-1]   =   x[i-1]*x[i-1]*x[i-1];   
		//y[i-1]       =   x[i-1] ;   
		/*計算x[i-1]點對應的y值作為擬合已知值*/   
	}

	Smooth(x, y, a, n, m, &dt1, &dt2, &dt3);                                         /*調用擬合函數*/   

	for(i = 1; i <= m; i++)   

	{
		printf("a[%d]   =   %.10f\n", (i-1), a[i-1]);   
	}

	printf("擬合多項式與數據點偏差的平方和為：\n");   

	printf("%.10e\n",dt1);   

	printf("擬合多項式與數據點偏差的絕對值之和為：\n");   

	printf("%.10e\n",dt2);   

	printf("擬合多項式與數據點偏差的絕對值最大值為：\n");   

	printf("%.10e\n",dt3);   

	free(x);/*釋放存儲空間*/   

	free(y);/*釋放存儲空間*/   

	free(a);/*釋放存儲空間*/   

}   



void Smooth(double   *x,double   *y,double   *a,int   n,int   m,double   *dt1,double   *dt2,double   *dt3   )   
//   
//double   *x; /*實型一維數組，輸入參數，存放節點的xi值*/   

//double   *y; /*實型一維數組，輸入參數，存放節點的yi值*/   
//   
//double   *a; /*雙精度實型一維數組，長度為m。返回m - 1次擬合多項式的m個系數*/   
//   
//int   n;                                                                                               /*整型變量，輸入參數，給定數據點的個數*/   
//   
//int   m;                                                                                               /*整型變量，輸入參數，擬合多項式的項數*/   
//   
//double   *dt1; /*實型變量，輸出參數，擬合多項式與數據點偏差的平方和*/   
//   
//double   *dt2; /*實型變量，輸出參數，擬合多項式與數據點偏差的絕對值之和*/   
//   
//double   *dt3; /*實型變量，輸出參數，擬合多項式與數據點偏差的絕對值最大值*/   

{   
	int   i,j,k;   
	double   *s,*t,*b,z,d1,p,c,d2,g,q,dt;   

	/*分別為s,t,b分配存貯空間*/   
	s = (double*)calloc(n,sizeof(double));   
	t = (double*)calloc(n,sizeof(double));   
	b = (double*)calloc(n,sizeof(double));   
	z = 0;   
	for(i=1;i<=n;i++)   
	{
		z=z+x[i-1]/n;     /*z為各個x的平均值*/   
	}

	b[0]=1;   
	d1=n;   
	p=0;   
	c=0;   
	for(i=1;i<=n;i++)   
	{   
		p=p+x[i-1]-z;   
		c=c+y[i-1];   
	}   
	c=c/d1;   
	p=p/d1;   
	a[0]=c*b[0];   
	if(m>1)   
	{   
		t[1]=1;   
		t[0]=-p;   
		d2=0;   
		c=0;   
		g=0;   
		for(i=1;i<=n;i++)   
		{   
			q=x[i-1]-z-p;   
			d2=d2+q*q;   
			c=y[i-1]*q+c;   
			g=(x[i-1]-z)*q*q+g;   
		}   
		c=c/d2;   
		p=g/d2;   
		q=d2/d1;   
		d1=d2;   
		a[1]=c*t[1];   
		a[0]=c*t[0]+a[0];   
	}   

	for(j=3;j<=m;j++)   
	{   
		s[j-1]=t[j-2];   
		s[j-2]=-p*t[j-2]+t[j-3];   
		if(j>=4)   
			for(k=j-2;k>=2;k--)   
				s[k-1]=-p*t[k-1]+t[k-2]-q*b[k-1];   
		s[0]=-p*t[0]-q*b[0];   
		d2=0;   
		c=0;   
		g=0;   
		for(i=1;i<=n;i++)   
		{   
			q=s[j-1];   
			for(k=j-1;k>=1;k--)   
				q=q*(x[i-1]-z)+s[k-1];   
			d2=d2+q*q;   
			c=y[i-1]*q+c;   
			g=(x[i-1]-z)*q*q+g;   
		}   
		c=c/d2;   
		p=g/d2;   
		q=d2/d1;   
		d1=d2;  
		a[j-1]=c*s[j-1];   
		t[j-1]=s[j-1];   
		for(k=j-1;k>=1;k--)   
		{   
			a[k-1]=c*s[k-1]+a[k-1];   
			b[k-1]=t[k-1];   
			t[k-1]=s[k-1];   
		}   
	}   
	*dt1=0;   
	*dt2=0;   
	*dt3=0;   
	for(i=1;i<=n;i++)   
	{   
		q=a[m-1];   
		for(k=m-1;k>=1;k--)   
			q=q*(x[i-1]-z)+a[k-1];   
		dt=q-y[i-1];   
		if(fabs(dt)>*dt3)   
			*dt3=fabs(dt);   
		*dt1=*dt1+dt*dt;   
		*dt2=*dt2+fabs(dt);   
	}   

	/*釋放存儲空間*/   
	free(s);   
	free(t);   
	free(b);   
}