аппроксимация

Нужна программа для аппроксимации методом наименьших квадратов, подскажите где взять.

Вот, из старых запасов:

Program Ex1ForEGA;
  Uses Graph,PGraph;
  Const
    Nmax=50;
  Var
    a:array [0..1] of real;
    b:array [0..2] of real;
    c:array [0..3] of real;
    x,y:array [1..Nmax] of real;
    f:text;
    n,i,j:integer;
    u:array [1..4,1..4] of real;
    r:array [1..4] of real;
    MaxX,MinX,KX,MaxY,MinY,KY:real;
    St1,St2,St3:string;
  Function f1(x:real):real;
    Begin
      f1:=a[0]+a[1]*x;
    End;
  Function f2(x:real):real;
    Begin
      f2:=b[0]+b[1]*x+b[2]*x*x;
    End;
  Function f3(x:real):real;
    Begin
      f3:=c[0]+c[1]*x+c[2]*x*x+c[3]*x*x*x;
    End;
  Procedure PStep(n1,n2:byte);
    Var
      i,j:byte;
    Begin
      if n1=n2 then
       begin
    r[n1]:=r[n1]/u[n1,n1];
    u[n1,n1]:=1;
    Exit;
       end;
      r[n1]:=r[n1]/u[n1,n1];
      For i:=n1+1 to n2 do
       u[n1,i]:=u[n1,i]/u[n1,n1];
      u[n1,n1]:=1;
      For i:=n1+1 to n2 do
       begin
    For j:=n1+1 to n2 do
     u[i,j]:=u[i,j]-u[i,n1]*u[n1,j];
    r[i]:=r[i]-u[i,n1]*r[n1];
    u[i,n1]:=0;
       end;
      PStep(n1+1,n2);
    End;
  Procedure OStep(n:byte);
    Var
      Sv:real;
      i,j,k:byte;
    Begin
      for i:=n downto 2 do
       For j:=1 to i-1 do
    begin
     Sv:=u[j,i];
     For k:=1 to n do
      u[j,k]:=u[j,k]-Sv*u[i,k];
     r[j]:=r[j]-Sv*r[i];
    end;
    End;
 
  Begin
    Assign(f,'EX1.DAT'{ParamStr(1)});
    Reset(f);
    Readln(f,n);
    If n>Nmax then Exit;
    For i:=1 to n do
     Readln(f,x[i],y[i]);
    n:=n-1;
    For i:=1 to 2 do
     For j:=1 to 2 do
      begin
       u[i,j]:=0;
       r[i]:=0;
      end;
    For i:=1 to n do
     begin
      u[1,1]:=u[1,1]+1;
      u[1,2]:=u[1,2]+x[i];
      u[2,1]:=u[2,1]+x[i];
      u[2,2]:=u[2,2]+x[i]*x[i];
      r[1]:=r[1]+y[i];
      r[2]:=r[2]+x[i]*y[i];
     end;
    PStep(1,2);
    OStep(2);
    For i:=0 to 1 do
     a[i]:=r[i+1];
    For i:=1 to 3 do
     For j:=1 to 3 do
      begin
       u[i,j]:=0;
       r[i]:=0;
      end;
    For i:=1 to n do
     begin
      u[1,1]:=u[1,1]+1;
      u[1,2]:=u[1,2]+x[i];
      u[1,3]:=u[1,3]+x[i]*x[i];
      u[2,1]:=u[2,1]+x[i];
      u[2,2]:=u[2,2]+x[i]*x[i];
      u[2,3]:=u[2,3]+x[i]*x[i]*x[i];
      u[3,1]:=u[3,1]+x[i]*x[i];
      u[3,2]:=u[3,2]+x[i]*x[i]*x[i];
      u[3,3]:=u[3,3]+x[i]*x[i]*x[i]*x[i];
      r[1]:=r[1]+y[i];
      r[2]:=r[2]+x[i]*y[i];
      r[3]:=r[3]+x[i]*x[i]*y[i];
     end;
    PStep(1,3);
    OStep(3);
    For i:=0 to 2 do
     b[i]:=r[i+1];
    For i:=1 to 4 do
     For j:=1 to 4 do
      begin
       u[i,j]:=0;
       r[i]:=0;
      end;
    For i:=1 to n do
     begin
      u[1,1]:=u[1,1]+1;
      u[1,2]:=u[1,2]+x[i];
      u[1,3]:=u[1,3]+x[i]*x[i];
      u[1,4]:=u[1,4]+x[i]*x[i]*x[i];
      u[2,1]:=u[2,1]+x[i];
      u[2,2]:=u[2,2]+x[i]*x[i];
      u[2,3]:=u[2,3]+x[i]*x[i]*x[i];
      u[2,4]:=u[3,3]+x[i]*x[i]*x[i]*x[i];
      u[3,1]:=u[3,1]+x[i]*x[i];
      u[3,2]:=u[3,2]+x[i]*x[i]*x[i];
      u[3,3]:=u[3,3]+x[i]*x[i]*x[i]*x[i];
      u[3,4]:=u[3,4]+x[i]*x[i]*x[i]*x[i]*x[i];
      u[4,1]:=u[4,1]+x[i]*x[i]*x[i];
      u[4,2]:=u[4,2]+x[i]*x[i]*x[i]*x[i];
      u[4,3]:=u[4,3]+x[i]*x[i]*x[i]*x[i]*x[i];
      u[4,4]:=u[4,4]+x[i]*x[i]*x[i]*x[i]*x[i]*x[i];
      r[1]:=r[1]+y[i];
      r[2]:=r[2]+x[i]*y[i];
      r[3]:=r[3]+x[i]*x[i]*y[i];
      r[4]:=r[4]+x[i]*x[i]*x[i]*y[i];
     end;
    PStep(1,4);
    OStep(4);
    For i:=0 to 3 do
     c[i]:=r[i+1];
    n:=n+1;
    MaxX:=x[1];
    MinX:=x[1];
    For i:=2 to n do
     begin
      if MaxX<x[i] then MaxX:=x[i];
      if MinX>x[i] then MinX:=x[i];
     end;
    KX:=460/(MaxX-MinX);
    MaxY:=y[1];
    MinY:=y[1];
    For i:=2 to n do
     begin
      if MaxY<y[i] then MaxY:=y[i];
      if MinY>y[i] then MinY:=y[i];
     end;
    For i:=30 to 490 do
     begin
      if MaxY<f1((i-30)/KX+MinX) then MaxY:=f1((i-30)/KX+MinX);
      if MinY>f1((i-30)/KX+MinX) then MinY:=f1((i-30)/KX+MinX);
      if MaxY<f2((i-30)/KX+MinX) then MaxY:=f2((i-30)/KX+MinX);
      if MinY>f2((i-30)/KX+MinX) then MinY:=f2((i-30)/KX+MinX);
      if MaxY<f3((i-30)/KX+MinX) then MaxY:=f3((i-30)/KX+MinX);
      if MinY>f3((i-30)/KX+MinX) then MinY:=f3((i-30)/KX+MinX);
     end;
    KY:=420/(MaxY-MinY);
    OpenGraph;
    SetColor(White);
    Rectangle(20,20,500,460);
    Rectangle(500,20,620,460);
    Rectangle(500,20,620,70);
    Rectangle(500,140,620,230);
    if MinX*MaxX<=0 then
     begin
      Line(30-round(KX*MinX),30,30-round(KX*MinX),450);
      Line(30-round(KX*MinX),30,25-round(KX*MinX),50);
      Line(30-round(KX*MinX),30,35-round(KX*MinX),50);
      Line(25-round(KX*MinX),50,35-round(KX*MinX),50);
     end;
    if MinY*MaxY<=0 then
     begin
      Line(30,450+round(KY*MinY),490,450+round(KY*MinY));
      Line(490,450+round(KY*MinY),470,455+round(KY*MinY));
      Line(490,450+round(KY*MinY),470,445+round(KY*MinY));
      Line(470,455+round(KY*MinY),470,445+round(KY*MinY));
     end;
    SetColor(Red);
    Line(30,450-round(KY*(f1(MinX)-MinY)),
         490,450-round(KY*(f1(MaxX)-MinY)));
    SetColor(Green);
    MoveTo(30,450-round(KY*(f2(MinX)-MinY)));
    For i:=31 to 490 do
     begin
      LineTo(i,450-round(KY*(f2((i-30)/KX+MinX)-MinY)));
     end;
    SetColor(Blue);
    MoveTo(30,450-round(KY*(f3(MinX)-MinY)));
    For i:=31 to 490 do
     begin
      LineTo(i,450-round(KY*(f3((i-30)/KX+MinX)-MinY)));
     end;
    SetColor(Yellow);
    SetFillStyle(1,Yellow);
    For i:=1 to n do
     begin
      FillEllipse(30+round(KX*(x[i]-MinX)),450-round(KY*(y[i]-MinY)),2,2);
     end;
    SetColor(Red);
    For i:=0 to 1 do
     begin
      Str(i,St1);
      Str(a[i]:-10:5,St2);
      St3:='a'+St1+'='+St2;
      OutTextXY(510,i*20+30,St3);
     end;
    SetColor(Green);
    For i:=0 to 2 do
     begin
      Str(i,St1);
      Str(b[i]:-10:5,St2);
      St3:='b'+St1+'='+St2;
      OutTextXY(510,i*20+80,St3);
     end;
    SetColor(Blue);
    For i:=0 to 3 do
     begin
      Str(i,St1);
      Str(c[i]:-10:5,St2);
      St3:='c'+St1+'='+St2;
      OutTextXY(510,i*20+150,St3);
     end;
    Readln;
    CloseGraph;
  End.