Метод наименьших квадратов (Pascal)

program MNK;   {Метод наименьших квадратов}

uses
   Crt;
const
   E=0.0001;
   N=3;
   Nmax=11;
type
   matrix = array[1..N,1..N+1] of real;
   vector = array[1..Nmax] of real;
   vector1 = array[1..Nmax-1] of real;
   vectorAlpha = array[1..N] of real;
const
   Y: vector
     =(3.240,3.984,4.473,4.724,4.772,4.662,4.445,4.166,3.869,3.585,3.339);
   X: vector
     =(  0.0,  1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,  9.0, 10.0);
   Ypr: vector1
     =(3.644,4.259,4.626,4.770,4.734,4.564,4.310,4.018,3.724,3.456);
   Xpr: vector1
     =(  0.5,  1.5,  2.5,  3.5,  4.5,  5.5,  6.5,  7.5,  8.5,  9.5);
var
   A: matrix;
   pog,max: real;
   KEY:BOOLEAN;
   S,M:REAL;

procedure Gauss(N:integer;A:matrix);  {Метод Гаусса}
var
   X: vectorAlpha;
   B: matrix;
   i,j: integer;

  Procedure PRINT;
     VAR
        I,J:INTEGER;
     BEGIN
        FOR I:=1 TO N DO
          BEGIN
            FOR J:=1 TO N+1 DO
              WRITE (A[I,J]:9:3, '    ');
            WRITELN
          END;
        WRITELN
     END;

   Procedure PR_XOD;                   {Прямой ход метода Гаусса}
     VAR
       K,I,J:INTEGER;
     BEGIN
       KEY:=TRUE;
       FOR K:=1 TO N-1 DO
         BEGIN
           IF A[K,K]<>0
             THEN
               FOR J:=K+1 TO N DO
                 BEGIN
                   M:=A[J,K]/A[K,K];   {Коэффициент М}
                   FOR I:=K TO N+1 DO
                     A[J,I]:=A[J,I]-A[K,I]*M;
                   A[J,K]:=0
                 END
             ELSE
               BEGIN
                 KEY:=FALSE;
                 WRITELN ('Алгоритм для решения не применим!');
               END
         END;
     END;

   Procedure PROVERKA;                {Проверка на наличие решения}
     BEGIN
       IF ABS(A[N,N])<E
         THEN
           IF ABS(A[N,N+1])<E
             THEN
               BEGIN
                 KEY:=FALSE;
                 WRITE ('Система имеет бесконечное число решений!')
               END
             ELSE
               BEGIN
                 KEY:=FALSE;
                 WRITE ('‘Система не имеет решений!');
               END
     END;

   Procedure OB_XOD;                   {Обратный ход метода Гаусса}
     VAR
       K,I,J:INTEGER;
     BEGIN
       IF KEY=TRUE
         THEN
           BEGIN
             FOR K:=N DOWNTO 1 DO
               BEGIN
                 S:=A[K,N+1];
                 FOR J:=N DOWNTO K+1 DO
                   S:=S-A[K,J]*X[J];
                 X[K]:=S/A[K,K];
               END;
             WRITELN;
             WRITELN ('Корни системы уравнений:');
             FOR I:=1 TO N DO
               WRITELN('X[',I,']=',X:5:3)
           END
         ELSE
           EXIT
     END;

procedure Diag;
var
   Determinant: real;
   i: integer;
begin
   Determinant:=1;
   for i:=1 to N do
      Determinant:=Determinant*A[i,i];
   writeln;
   writeln('„ҐвҐа¬Ё***в=',Determinant:9:4);
end;

   BEGIN
     WRITELN ('                  Исходная матрица:');
     PRINT;
     WRITELN;
     PR_XOD;
     PROVERKA;
     WRITELN;
     WRITELN;
     WRITELN ('                 Треугольная матрица:');
     PRINT;
     Diag;
     OB_XOD;
   END;

function power(x:real;n:word):real;
var
   result: real;
begin
   result:=1.0;
   while (n>0) do
      begin
         result:=result*x;
         dec(n);
      end;
   power:=result;
end;

var
   i,j,l: integer;

BEGIN                                 {Основная программа}
   clrscr;
   for i:=1 to N do
      for j:=1 to N+1 do
         a[i,j] := 0;
   for i:=1 to N do
      for j:=1 to N do
         for l:=1 to Nmax do
            a[i,j] := a[i,j] + power(X[l],i+j-2);
   for i:=1 to N do
      for l:=1 to Nmax do
         a[i,N+1] := a[i,N+1] + power(X[l],i-1)*(Y[l]);
   Gauss(N,A);
   Readkey
end.