{$A+,B-,D+,E-,F-,G+,I+,L+,N+,O-,P-,Q-,R-,S-,T-,V+,X+,Y+}
Program _Approc_from_file;{Аппроксимация функций методом наименьших квадратов}
{$E-,N+,G+,D+}
Uses Graph, CRT;
type
 Sto = array [ 1..70 ] of Extended;
const
  n = 1368 div 2 - 1;
  x0 : Extended = 1.0;
  xn : Extended = 4.0;
  Yb : Integer = 464;
var
  m, nn, k, p, s : Byte;
  ndot : Word;
  Dr, R, e : Integer;
  y : array [0..n] of Extended;
  xxx : array [0..n] of Extended;
  dx, ymax, ymin, x : Extended;
  A          : array [ 1..70 ] of Sto;{Массив }
  B          : Sto ;
  i, j : Longint;
  ff : Text;
Procedure Swap(p, q: Byte);
var
 j : Byte;
 X : Extended;
begin
 for j:=1 to nn do
   begin
    X := A[p,j];
    A[p,j] := A[q,j];
    A[q,j] := X;
   end;
   X := B[p];
   B[p] := B[q];
   B[q] := X;
end;
Procedure Divid(p : Byte; k : Extended);
var
 j : Byte;
begin
 for j := 1 to nn do
   A[p, j] := A[p, j] / k;
 B[p] := B[p] / k;
end;

Procedure Sub(p, q : Byte);
var
 j : Byte;
begin
 for j := 1 to nn do
   A[p, j] := A[p, j] - A[q, j];
 B[p] := B[p] - B[q];
end;

Function CheckNullStr(p : Byte): Boolean;
var
 i : Byte;
 k : Extended;
begin
 k := 0;
 for i := 1 to nn do
   k:=k + Abs (A[p, i]);
 CheckNullStr := ((k=0) and (B[p]=0));
end;

Function BadStr(p : Byte) : Boolean;
var
 i : Byte;
 k : Extended;
begin
 k := 0;
 for i := 1 to nn do
   k:=k + Abs (A[p, i]);
 BadStr := ((k=0) and (B[p]<>0));
end;

Procedure KillNullStr(r : Byte);
begin
 Swap(r, m);
 dec(m);
end;

Function IsZeroDiag : Boolean;
var
 i : Byte;
 q : Boolean;
begin
 q := false;
 for i := 1 to m + 1 do
   q := q or (A[i,i] = 0) ;
 IsZeroDiag := q;
end;

Procedure FindNZD ;
var
 i : Byte;
begin
 for i := k+1 to m + 1 do
   if A[i,k] <> 0.0 then
     begin
       Swap(k, i);
     end;
end;
Procedure SolveMattr;
var
  i : Byte;
Label
  Loop;
begin
k := 1;
Loop:
 While k <= m+1 do     {???}
   begin
    for i := 1 to m+1 do
      if not BadStr(i) then
         begin
           if CheckNullStr(i) then
             begin
               KillNullStr(i);
               goto loop;
             end;
         end
     else {ViewResult(0)};
  FindNZD;
      p:=0;
      for i := k to nn do
                    if A[k,i] = 0 then inc(p)
                                  else Break;
      for i := 1 to m + 1 do
        if A[i,k+p]<>0 then
          Divid(i, A[i,k+p]);
      for i :=1 to m + 1 do
        if (i<>(k+p)) and (A[i,k+p]<>0.0) then
          Sub(i,k+p);
    inc(k);
   end;
  for i:=1 to m + 1 do if A[i,i]<>0 then
    Divid(i,A[i,i]);
 i := 1;
 While i <= m + 1 do
   begin
{    if BadStr(i) then {ViewResult(0);}
    if CheckNullStr(i) then KillNullStr(i);
    inc(i);
   end;
end;
Function Xdg(x, dg : Extended):Extended;
begin
  Xdg := exp(dg*ln(x));
end;
Function Fi(x:Extended):Extended;
var
  tmp : Extended;
  k : Word;
begin
  tmp := 0;
  for k := 0 to m do tmp := tmp + B[k+1]*Xdg(x, 1.0*k);
  Fi := tmp;
end;
Procedure CalledKeyProc;
begin
  x := x0;
  dx := (xn-x0)/1000.0;
  MoveTo(20,Yb-Round((Fi(x0)-ymin)/(ymax-ymin)*450));
  SetColor(LightRed);
  While x <= xn do
  begin
    LineTo(20+Round((x-x0)/(xn-x0)*600),Yb-Round((Fi(x)-ymin)/(ymax-ymin)*450));
    x := x + dx;
  end;
end;
Function Max(Ar : array of Extended) : Extended;
var
  Mx : Word;
  j : Word;
begin
  Mx := 0;
  for j := 0 to ndot do
   if Ar[Mx]Ar[j] then Mn := j;
  Min := Ar[Mn];
end;
Procedure SetMattr;
var
  i, k : Word;
  tmp : Extended;
begin
  for s := 0 to m do {по строкам}
    begin
      for k := 0 to m do  {По элементам строки}
        begin
          tmp := 0.0;
          for i := 0 to ndot do
            tmp := tmp + Xdg(xxx[i], 1.0*(k+s));
          A[s+1, k+1] := tmp;
        end;
          tmp := 0.0;
          for i := 0 to ndot do
            tmp := tmp + y[i]*Xdg(xxx[i], 1.0*s);
      B[s+1] := tmp;
    end;
end;
Procedure ScreenMattr;
var
 i, j : Byte;
begin
Writeln;
 for i := 1 to m+1 do
   begin
{     for j := 1 to nn do
       Write(A[i, j]:1 :3,' ');}
    WriteLn('a',i-1,':| ',B[i]:1 :3);
   end;
ReadKey;
end;
Function Dispersia : Extended;
var
  i : Longint;
  tmp : Extended;
begin
  tmp := 0;  {!!!!!!!!!!!!!}
  for i := 0 to ndot do
    tmp := tmp + sqr(y[i]-Fi(xxx[i]));
  Dispersia := tmp/(ndot+1);
end;
begin
  if ParamCount=0 then
  begin
    WriteLn('Введите имя файла, содержащего данные...');
    WriteLn('Структура файла:');
    WriteLn('{Число точек}');
    WriteLn('{Степень полинома}');
    WriteLn('... {координаты точек} ...');
    Exit;
  end;
  Assign(ff, ParamStr(1));
  Reset(ff);
  ReadLn(ff, ndot);
  if ndot > n then
  begin
   Close(ff);
   WriteLn('Слишком много точек в файле!!!');
   Halt;
  end;
  Randomize;
  ClrScr;
  ReadLn(ff, m);
  ReadLn(ff, x0);
  ReadLn(ff, xn);
  nn := m+1;
  dec(ndot);
  for i := 0 to ndot do
  begin
    ReadLn(ff, xxx[i]);
    ReadLn(ff, y[i]);
  end;
  Close(ff);
  dx := (xn-x0)/ndot;
  WriteLn('Писло точек ', ndot+1);
  WriteLn('Степень полинома ', m);
  ymax := Max(y);
  ymin := Min(y);
  x0 := Min(xxx);
  xn := Max(xxx);
  if (xn <= 0) or (x0 <= 0) then
   begin
    for i := 0 to ndot do
      xxx[i] := xxx[i] + Abs(xn)+Abs(x0);
    x0 := Min(xxx);
    xn := Max(xxx);
   end;
  WriteLn('x0 ', x0);
  WriteLn('xn ', xn);
  if (x0 = 0) and (xn = 0) then
  begin
    WriteLn('Длина отрезка по оси абсцисс равна 0 !!!');
    Halt;
  end;
  SetMattr;
  SolveMattr;
  ScreenMattr;
  Dr := Detect;        {Инициализируем графику}
  InitGraph(Dr, R, '');
  e :=  GraphResult;
  if e <> grOk then {Если ошибка}
  begin
    WriteLn('Ошибка инициализации графики!');
    WriteLn(GraphErrorMsg(e));
    ReadKey;
    Halt($FFFF);
  end;
  MoveTo(20+Round((xxx[0]-x0)/(xn-x0)*600), Yb-Round((y[0]-ymin)/(ymax-ymin)*450));
  SetColor(LightGreen);
  Circle(20+Round((xxx[0]-x0)/(xn-x0)*600), Yb-Round((y[0]-ymin)/(ymax-ymin)*450), 1);
  for i := 1 to ndot do
  begin
    SetColor(Yellow);
    LineTo(20+Round((xxx[i]-x0)/(xn-x0)*600), Yb-Round((y[i]-ymin)/(ymax-ymin)*450));
    SetColor(LightGreen);
    Circle(20+Round((xxx[i]-x0)/(xn-x0)*600), Yb-Round((y[i]-ymin)/(ymax-ymin)*450), 1);
  end;
  ReadKey;
  CalledKeyProc;
  ReadKey;
  CloseGraph;
  dx := (xn-x0)/ndot;
  for i := 0 to ndot do
    WriteLn(i, ' : ', x0+dx*i:1:3, ' -> ', sqr(Fi(x0+dx*i)-y[i]));
  Write('Дисперсия равна: ');
  WriteLn(Dispersia:1);
  Write('Среднеквадратичное отклонение равно: ');
  WriteLn(sqrt(Dispersia):1);
  ReadKey;
end.program