曲线拟合多项式

void gauss_solve(int n,double A[],double x[],double b[])
{
    int i,j,k,r;
    double max;
    for (k=0;k<n-1;k++)
    {
        max=fabs(A[k*n+k]); /*find maxmum*/
        r=k;
        for (i=k+1;i<n-1;i++)
        if (max<fabs(A[i*n+i]))
        {
            max=fabs(A[i*n+i]);
            r=i;
        }
        if (r!=k)
        for (i=0;i<n;i++) /*change array:A[k]&A[r] */
        {
            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;
        for (i=k+1;i<n;i++)
        {
            for (j=k+1;j<n;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];
        }
    }
    //for (i=n-1;i&gt;=0;x[i]/=A[i*n+i],i--) 
    for(i=n-1;i>=0;x[i]/=A[i*n+i],i--) 
        for(j=i+1,x[i]=b[i];j<n;j++)
            x[i]-=A[i*n+j]*x[j]; 
}

/*==================polyfit(n,x,y,poly_n,a)===================*/
/*=======拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n========*/
/*=====n是数据个数 xy是数据值 poly_n是多项式的项数======*/
/*===返回a0,a1,a2,……a[poly_n]，系数比项数多一（常数项）=====*/
void polyfit(int n,double x[],double y[],int poly_n,double a[])
{
    int i,j;
    double *tempx,*tempy,*sumxx,*sumxy,*ata;
    void gauss_solve(int n,double A[],double x[],double b[]);
    tempx=calloc(n,sizeof(double));if(tempx == NULL)     printf("err");
    
    tempx=malloc(n*sizeof(double));if(tempx == NULL)     printf("err");
    printf("6:%d\r\n",n*sizeof(double));
    
    sumxx=calloc(poly_n*2+1,sizeof(double));if(sumxx == NULL)     printf("err"); printf("%d\r\n",sizeof(sumxx));
    printf("2:%d\r\n",(poly_n*2+1)*sizeof(double));
    
    tempy=calloc(n,sizeof(double));if(tempy == NULL)     printf("err");
    printf("3:%d\r\n",(n)*sizeof(double));
    
    sumxy=calloc(poly_n+1,sizeof(double));if(sumxy == NULL)     printf("err");
    printf("4:%d\r\n",(poly_n+1)*sizeof(double));
    
    ata=calloc((poly_n+1)*(poly_n+1),sizeof(double));if(ata == NULL)     printf("err");
    printf("5:%d\r\n",((poly_n+1)*(poly_n+1))*sizeof(double));
    
    
    for (i=0;i<n;i++)
    {
        tempx[i]=1;
        tempy[i]=y[i];
    }
    for (i=0;i<2*poly_n+1;i++)
    for (sumxx[i]=0,j=0;j<n;j++)
    {
        sumxx[i]+=tempx[j];
        tempx[j]*=x[j];
    }
    for (i=0;i<poly_n+1;i++)
    for (sumxy[i]=0,j=0;j<n;j++)
    {
        sumxy[i]+=tempy[j];
        tempy[j]*=x[j];
    }
    for (i=0;i<poly_n+1;i++)
    for (j=0;j<poly_n+1;j++)
        ata[i*(poly_n+1)+j]=sumxx[i+j];
    gauss_solve(poly_n+1,ata,a,sumxy);
    free(tempx);
    free(sumxx);
    free(tempy);
    free(sumxy);
    free(ata);
}

//double x[7]={1,2,3,4,6,7,8},y[7]={2,3,6,7,5,3,2}; double a[4]; 
void PolynomialFitting(double *x,double *y,int n,int poly_n)
{
int i;//n=7,poly_n=3; 
double a[4]; 

polyfit(n,x,y,poly_n,a);
for (i=0;i<poly_n+1;i++)/*这里是升序排列，Matlab是降序排列*/      
printf("a[%d]=%g\r\n",i,a[i]);
}