//曲線擬合,例題4.4

#include "stdio.h"
#include "math.h"
#define num 10


float neiji(float b[num],float c[num])      //內積函數
{  int p;
   float nj=0;
    for (p=1;p<num;p++)
          nj+=c[p]*b[p];
    return nj;
}
float s[num],x[num],y[num],fai[num][num],afa[num],beida[num],a[num],xfai[num],yd[num],max,pcpfh;


void main()
{ int i,j,k,n,index,flag;
  char conti;
  conti=' ';
 

  
    printf("請輸入已知點的個數n=\n");
    scanf("%d",&n);

	printf("請輸入x和y:");
   for(i=1;i<=n;i++)
   {   printf("x[%d]=",i);
	   scanf("%f",&x[i]);
       printf("y[%d]=",i);
	   scanf("%f",&y[i]);
   }   
while(conti==' ')
{	 printf("請輸入擬和次數=");
    scanf("%d",&index);  
    pcpfh=0;
   afa[1]=0;
   a[0]=0;
   for(i=1;i<=n;i++)
   {afa[1]+=x[i];
    a[0]+=y[i];
	fai[0][i]=1;
   }
   afa[1]=afa[1]/n;
   a[0]=a[0]/n;   
   
   for (i=1;i<=n;i++)
   {
    fai[1][i]=x[i]-afa[1];
   }    

   a[1]=neiji(fai[1],y)/neiji(fai[1],fai[1]);
 
   for(k=1;k<index;k++)
   {   for(i=1;i<=n;i++)
           xfai[i]=x[i]*fai[k][i];

	   afa[k+1]=neiji(fai[k],xfai)/neiji(fai[k],fai[k]);
        beida[k]=neiji(fai[k],fai[k])/neiji(fai[k-1],fai[k-1]);
	   for(j=1;j<=n;j++)
	        fai[k+1][j]=(x[j]-afa[k+1])*fai[k][j]-beida[k]*fai[k-1][j];
      a[k+1]=neiji(fai[k+1],y)/neiji(fai[k+1],fai[k+1]);
	 
   }
   
   printf("%d次擬和結果為\n",index);
   for(i=0;i<=index;i++)   
    printf("a[%d]=%f\n",i,a[i]);
   for(i=1;i<=index;i++)
	   printf("afa[%d]=%f\n",i,afa[i]);
   for(i=1;i<index;i++)
       printf("beida[%d]=%f\n",i,beida[i]);   

   for(i=1;i<=n;i++)
   {	for(k=0;k<=index;k++)
	      s[i]+=a[k]*fai[k][i];
       yd[i]=float(fabs(y[i]-s[i]));
	   pcpfh+=yd[i]*yd[i];
       s[i]=0;
   }

   max=0;
   for(i=1;i<=n;i++)
	   if(yd[i]>max)
	   {max=yd[i];
	    flag=i;
	   }

   
   printf("當x=%f時,偏差最大=%f,偏差平方和為%f\n",x[flag],max,pcpfh);

   printf("繼續擬和請按space,按其他鍵退出");

   conti=getchar();
   conti=getchar();

  }
}