用C语言实现polyfit函数

用C语言实现polyfit多项式拟合,已知离散点上的数据集，即已知在点集上的函数值，构造一个解析函数（其图形为一曲线）使在原离散点上尽可能接近给定的值。

// shacha.cpp : Defines the entry point for the console application.
//

#include “stdafx.h“
//#include “system.h“
#include 

// polyfit.cpp : 定义控制台应用程序的入口点。


/**********************************************
 *Author  :wacs5
 *DATE  :20090408（YYYMMDD）
 *Functtion  :多项式拟合polyfit
 **********************************************/
#include 
#include 
#include 
#include 
void main（）
{
 int in=7poly_n=2;
 double x[7]y[7];
 double a[3];
 printf（“请输入%d个x“n）;
 printf（“\n“）;
 scanf（“%lf %lf %lf %lf %lf %lf %lf“ x x+1 x+2 x+3 x+4 x+5 x+6）;
 printf（“请输入%d个y“n）;
 printf（“\n“）;
 scanf（“%lf %lf %lf %lf %lf %lf %lf“ y y+1 y+2 y+3 y+4 y+5 y+6）;
 void polyfit（int ndouble x[]double y[]int poly_ndouble a[]）;
 system（“cls“）;
 polyfit（nxypoly_na）;
 
 for （i=0;i     printf（“a[%d]=%g\n“ia[i]）;
 getch（）;
}


/*==================polyfit（nxypoly_na）===================*/
/*=======拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n========*/
/*=====n是数据个数  xy是数据值  poly_n是多项式的项数======*/
/*===返回a0a1a2……a[poly_n]，系数比项数多一（常数项）=====*/
void polyfit（int ndouble x[]double y[]int poly_ndouble a[]）
{
 int ij;
 double *tempx*tempy*sumxx*sumxy*ata;
 void gauss_solve（int ndouble A[]double x[]double b[]）;
 tempx=（double *）calloc（nsizeof（double））;
 sumxx=（double *）calloc（poly_n*2+1sizeof（double））;
 tempy=（double *）calloc（nsizeof（double））;
 sumxy=（double *）calloc（poly_n+1sizeof（double））;
 ata=（double *）calloc（（poly_n+1）*（poly_n+1）sizeof（double））;
 for （i=0;i     {
      tempx[i]=1;
      tempy[i]=y[i];
     }
 for （i=0;i<2*poly_n+1;i++）
     for （sumxx[i]=0j=0;j   {
    sumxx[i]+=tempx[j];
    tempx[j]*=x[j];
   }
 for （i=0;i    for （sumxy[i]=0j=0;j   {
    sumxy[i]+=tempy[j];
    tempy[j]*=x[j];
   }
 for （i=0;i     for （j=0;j  ata[i*（poly_n+1）+j]=sumxx[i+j];
 gauss_solve（poly_n+1ataasumxy）;
 
 free（tempx）;
 free（sumxx）;
 free（tempy）;
 free（sumxy）;
 free（ata）;
}

void gauss_solve（int ndouble A[]double x[]double b[]）
{
 int ijkr;
 double max;
 //把ata拥有最大斜对角值的行换到第一行，从大到小依次排列，把b按大小依次排列
 for （k=0;k     {
      max=fabs（A[k*n+k]）; /*find maxmum，假设右斜线都是最大值*/
      r=k;
      for （i=k+1;i   if （max      {
       max=fabs（A[i*n+i]）;	//假设此时max小于其右边斜对角的哪个值，将那个最大值赋予max，并记下行号
       r=i;
      }
      if （r!=k）				//假设斜角上的值在当前行不是最大值的话
  for （i=0;i      {
       max=A[k*n+i];
       A[k*n+i]=A[r*n+i];	
       A[r*n+i]=max;		//把最大值那行的值与当前行互换
      }
      max=b[k];                    /*change array:b[k]&b[r]     */
      b[k]=b[r];		
      b[r]=max;				//ata的一行对应b的一个数，将变换后的ata与b的序号对应
      for （i=k+1;i   {
    for （j=k+1;j        A[i*n+j]-=A[i*n+k]*A[k*n+j]/A[k*n+k];
    b[i]-=A[i*n+k]*b[k]/A[k*n+k

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2013-05-23 22:13  polyfit\
     目录           0  2013-05-23 22:13  polyfit\debug\
     文件       40960  2013-05-22 22:33  polyfit\debug\polyfit.exe
     文件      329428  2013-05-22 22:33  polyfit\debug\polyfit.ilk
     文件      355328  2013-05-22 22:33  polyfit\debug\polyfit.pdb
     目录           0  2013-08-31 18:36  polyfit\polyfit\
     目录           0  2013-05-23 22:13  polyfit\polyfit\Debug\
     文件        8906  2013-05-22 22:33  polyfit\polyfit\Debug\BuildLog.htm
     文件          65  2013-05-22 22:33  polyfit\polyfit\Debug\mt.dep
     文件         403  2013-05-21 20:38  polyfit\polyfit\Debug\polyfit.exe.embed.manifest
     文件         468  2013-05-21 20:38  polyfit\polyfit\Debug\polyfit.exe.embed.manifest.res
     文件         385  2013-05-22 22:33  polyfit\polyfit\Debug\polyfit.exe.intermediate.manifest
     文件       11208  2013-05-22 22:33  polyfit\polyfit\Debug\polyfit.obj
     文件     1048576  2013-05-21 20:38  polyfit\polyfit\Debug\polyfit.pch
     文件       10613  2013-05-21 20:38  polyfit\polyfit\Debug\stdafx.obj
     文件       68608  2013-05-22 22:33  polyfit\polyfit\Debug\vc80.idb
     文件      110592  2013-05-22 22:33  polyfit\polyfit\Debug\vc80.pdb
     文件         950  2013-05-21 17:08  polyfit\polyfit\ReadMe.txt
     文件        3339  2013-05-22 22:29  polyfit\polyfit\polyfit - 副本.cpp
     文件        3392  2013-05-22 22:33  polyfit\polyfit\polyfit.cpp
     文件        4487  2013-05-21 17:08  polyfit\polyfit\polyfit.vcproj
     文件        1419  2013-05-22 22:41  polyfit\polyfit\polyfit.vcproj.liguojun-PC.liguojun.user
     文件        4885  2013-05-25 00:46  polyfit\polyfit\polyfit1.cpp
     文件         212  2013-05-21 17:08  polyfit\polyfit\stdafx.cpp
     文件         276  2013-05-21 17:08  polyfit\polyfit\stdafx.h
     文件      633856  2013-05-22 22:41  polyfit\polyfit.ncb
     文件         886  2013-05-21 17:08  polyfit\polyfit.sln
     文件        7680  2013-05-22 22:41  polyfit\polyfit.suo