Полиномиальная аппроксимация функции методом наименьших квадратов

#include<stdio.h>
#include <vcl.h>
#pragma hdrstop
#include "Unit1.h"
#include <math.h>
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double X[100];
double Y[100];
int K,n=2;
bool Apr(double X[100], double Y[100], double a[100], int K, int n)
{
        bool Result=true;
        int i, j, k;
        int Order=K+1;
        double *b = new double[Order];
        memset(b, 0x00,sizeof(double)*(Order));
        double **sums = new double*[Order];
        for(i=0; i<Order; i++)
        {
                sums[i] = new double[Order];
        }
        for(i=0; i<Order; i++)
        {
                for(j=0; j<Order; j++)
                {
                        sums[i][j] = 0;
                        for(k=0; k<n; k++)
                        {
                                sums[i][j] += pow(X[k], i+j);
                        }
                }
        }
        for(i=0; i<Order; i++)
        {
                for(j=0; j<n; j++)
                {
                        b[i] += pow(X[j], i) * Y[j];
                }
        }
        double temp=0;
        for(i=0; i<Order; i++)
        {
                if(sums[i][i]==0)
                {
                        for(j=0; j<Order; j++)
                        {
                                if(j==i)
                                        continue;
                                if(sums[j][i] !=0 && sums[i][j]!=0)
                                {
                                        for(k=0; k<Order; k++)
                                        {
                                                temp = sums[j][k];
                                                sums[j][k] = sums[i][k];
                                                sums[i][k] = temp;
                                        }
                                        temp = b[j];
                                        b[j] = b[i];
                                        b[i] = temp;
                                        break;
                                }
                        }
                }
        }
        for(k=0; k<Order; k++)
        {
                for(i=k+1; i<Order; i++)
                {
                        if(sums[k][k]==0)
                        {
                               Result = false;
                                goto CleanUp;
                        }
                        double M = sums[i][k] / sums[k][k];
                        for(j=k; j<Order; j++)
                        {
                                sums[i][j] -= M * sums[k][j];
                        }
                        b[i] -= M*b[k];
                 }
        }
 
        for(i=(Order)-1; i>=0; i--)
        {
                double s = 0;
                for(j = i; j<Order; j++)
                {
                        s = s + sums[i][j]*a[j];
                }
                a[i] = (b[i] - s) / sums[i][i];
        }
 
        CleanUp:
        delete [] b;
        for(i=0; i<Order; i++)
        {
                 delete [] sums[i];
        }
        delete [] sums;
        return Result;
 
}
 
 
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
        : TForm(Owner)
{
}
//---------------------------------------------------------------------------
 
 
void __fastcall TForm1::Button1Click(TObject *Sender)
{
     for(int i=0;i<n;i++)
        if(StringGrid1->Cells[i][0]!="")
           if(StringGrid1->Cells[i][1]!="")
             {
                X[i]=StrToFloat(StringGrid1->Cells[i+1][0]);
                Y[i]=StrToFloat(StringGrid1->Cells[i+1][1]);
             }
 
    for(int i=0;i<n;i++)
       Series1->AddXY(X[i],Y[i]);
 
     double l=StrToFloat(StringGrid1->Cells[1][0]);
     double m=StrToFloat(StringGrid1->Cells[n][0]);
 
        K=UpDown2->Position;
        double * a = new double[K+1];
        memset(a, 0,sizeof(double)*(K+1));
 
    bool Result=Apr(X,Y,a,K,n);
 
        if(Result)
        {
                for(double X=l;X<=m;X++)
                {
                        double Y=0;
                        for(int i=K;i>=1;i--)
                        Y+= a[i]*pow(X,i);
                        Y+=a[0];
                        Series3->AddXY(X,Y);
                }
        }
        delete [] a;
 
}
//---------------------------------------------------------------------------