#include<stdio.h>
#include<math.h>
#define MAX_N 20             //定義（x_i,y_i）的最大維數(shù)
typedef struct tagPOINT      //定義點(diǎn)結(jié)構(gòu)
{double x;
 double y;
}POINT;
int main()
{
 int n;
 int i,j,k,t,i_k,q;                      //q也是循環(huán)變量
 double s,a,b;
 POINT points[MAX_N+1];
 double l[MAX_N+1][MAX_N+1];             //設(shè)頂定初始值為零
 double m[MAX_N+1];                      //存儲(chǔ)正規(guī)方程組的右端項(xiàng)
 double p[2*MAX_N+1]={};                    //用p[k]存儲(chǔ)x_i的k次冪的和
 printf("請(qǐng)輸入被擬合點(diǎn)的數(shù)目,n∈(1,%d)：\n",MAX_N);
 scanf("%d",&n);
 printf("請(qǐng)輸入被擬合點(diǎn)（x_i,y_i），i=0,1,2,...,%d:\n",n);
 for(i=0;i<n;i++)
  scanf("%lf%lf",&points[i].x,&points[i].y);
 printf("請(qǐng)輸入擬合多項(xiàng)式的次數(shù)k(k∈[1,n]):\n");
 scanf("%d",&k);
 for(j=0;j<=2*k;j++)                  //計(jì)算x_i的k次冪的和
 {
  p[j]=0;
  for(i=0;i<n;i++)
  {
   p[j]=p[j]+pow(points[i].x,j);
  }
//  printf("%lf\n",p[j]);
 }
 for(j=0;j<=k;j++)                    //計(jì)算正規(guī)方程組的右端項(xiàng)
 {
  m[j]=0;
  for(i=0;i<n;i++)
  {
   m[j]=m[j]+points[i].y*pow(points[i].x,j);
  }
 }
 for(i=0;i<=k;i++)                    //計(jì)算正規(guī)方程組的系數(shù)矩陣
 {
  for(j=0;j<=k;j++)
  {
   t=i+j;
   l[i][j]=p[t];
  }
 }
 for(i=0;i<=k;i++)                   //輸出正規(guī)矩陣方程組
 {
  printf("\n");
  for(j=0;j<=k;j++)
  printf("%lf ",l[i][j]);
  printf(" %lf ",m[i]);
 }

//用高斯列主元消去法解正規(guī)方程組 ，即求擬合多項(xiàng)式的系數(shù)

 for(q=0;q<=k;q++)
 {
  s=0;                      //開(kāi)始按列選主元
  i_k=q;
  for(i=q;i<=k;i++)
  {
   if(fabs(l[i][q])>s)
   {
    s=fabs(l[i][q]);
    i_k=i;                  //記錄主元所在行的指標(biāo)
   } //if
  }  //for i
  if(s==0)
  {
   printf("系數(shù)矩陣奇異！\ndetA=0\n");
   break;
  } //if
  else if(i_k!=q)
  {
   for(j=q;j<=k;j++)       //換行     書(shū)上的算法這里錯(cuò)了！！！
   {
    a=l[i_k][j];
    l[i_k][j]=l[q][j];
    l[q][j]=a;
   } //for j
   b=m[i_k];
   m[i_k]=m[q];
   m[q]=b;
  } //else if
  for(i=q+1;i<=k;i++)           //計(jì)算乘子
  {
   l[i][q]=l[i][q]/l[q][q];
  } //for i
  for(i=q+1;i<=k;i++)           //消元
  {
   for(j=q+1;j<=k;j++)
   {
    l[i][j]=l[i][j]-l[i][q]*l[q][j];
   } //for j
   m[i]=m[i]-l[i][q]*m[q];
  }  //for i
 } //for q
  for(i=k;i>=0;i--)          //回代；臨時(shí)用b來(lái)存儲(chǔ)這里的和
  {
   b=0;
   for(j=i+1;j<=k;j++)
   {b+=l[i][j]*m[j];}
   m[i]=(m[i]-b)/l[i][i];
  } //for i
 printf("\n解為x_i=\n");     //輸出正規(guī)方程組的解，即擬合多項(xiàng)式的系數(shù)
 for(i=0;i<=k;i++)
 printf("%lf\n",m[i]);
 return 0;
}