program Matrix_Process;
(******************************************************
* Демонстрация работы с двумерными массивами          *
* (матрицами)                                         *
******************************************************)
{ uses power1, gamma1, integral;}

 const
  NMAX = 20;
  MMAX = 20;

type
  Matrix = array [1..NMAX, 1..MMAX] of extended;

var
  Funk, Matr: Matrix;
  n,m: integer;
  Str:string;
  file1,file2:string;
  s1,s2:char;

(***************************************************
            Процедура ввода матрицы из файла
****************************************************)
procedure EnterMatrixFromFile;

var
  i,j: integer;
  f:text;

begin
 assign(f,'matrixdat.txt');
 reset(f);
 i := 1;

 while not eof(f) do
   begin
   j := 1 ;
     while not eoln(f) do
         begin
             read(f,Matr[i, j]);
             {writeln('Matr[',i,',',j,']',Matr[i,j]);}
             inc(j);
         end;

       readln(f); inc(i);
   end;
    close(f);
    n:=i-1;  {Є®«ЁзҐбўв® бва®Є}
    m:=j-1;  {Є®«ЁзҐбвў® бв®«Ўж®ў}
    writeln('ђ §¬Ґа Ёб室­®© ¬ ваЁжл: ‘ва®Є n=',n,',','‘в®«Ўж®ў m=',m);
    end;

(******************************************************
Функция возведения действительного числа а в степень b
******************************************************)

Function Power (a, b : extended): extended; { a^b }

begin

if (a > 0) and (b<>0) then Power := exp(b * ln (a));
if (a < 0) and (b<>0) then Power := exp(b * ln (-a));
if b = 0 then Power := 1;
if a = 0 then Power := 0;

end;

(******************************************************
Процедура печати матрицы реальной размерности (n,m) на
экране, располагающая одну строку матрицы на одной строке
экрана
******************************************************)
procedure PrintMatrix (Var Matrix_: Matrix; nn, mm: integer);
var
  i, j: integer;
begin
  for i:=1 to nn do begin
    for j:=1 to mm do
     writeln('Matrix_[',i,',',j,']',Matrix_[i,j]);
    writeln;
  end;
end;

(******************************************************
Процедура печати матрицы реальной размерности (n,m) в файл,
располагающая одну строку матрицы на одной строке
******************************************************)
procedure PrintMatrixToFile (Var Matrix_: Matrix; Var String_: string; nn: integer);
var
  i, j: integer;
  g:text;
begin
 assign(g,String_);
 rewrite(g);

  for i:=1 to nn do
  begin
    for j:=1 to nn do
     begin
      write (g,matrix_[i,j]);
      write (g,' ');
     end;
    writeln(g);
  end;
  close(g);
  end;

(*******************************************************
                     Гамма-Функция Эйлера
*******************************************************)
(******************************************************)
function Gamma(x: extended): extended;
const
pi	= 3.1415926;

var
i,j	: integer;
y,gam	: extended;

begin		{ gamma function }
If x=0 then gamma:=1;
  if x>0.0 then
    begin
      y:=x+2.0;
      gam:=sqrt(2*pi/y)*exp(y*ln(y)+(1-1/(30*y*y))/(12*y)-y);
      Gamma:=gam/(x*(x+1))
    end
  else		{ x<0 }
    begin
      j:=0;
      y:=x;
      repeat
	j:=j+1;
	y:=y+1.0
      until y>0.0;
      gam:=gamma(y);		{ recursive call }
      for i:=0 to j-1 do
	gam:=gam/(x+i);
      Gamma:=gam
    end	{ x<0 }
end;		{ gamma function }
(******************************************************)

(**********************************************************************
Программа вычисления определенного интеграла
                tk

             I= S power(t,-1+m/2)*exp(-t/2) dt

                tn
по методу Симпсона.

Значение определенного интеграла по
формуле Симпсона вычисляется по формуле:

 ISimps=2*h/3*(0.5*Q(A)+2*Q(A+h)+Q(A+2*h)+2*Q(A+3*h)+...
                                                +2*Q(B-h)+0.5*Q(B))

где A=0 и B - нижняя и верхняя границы интервала интегрирования,
    N - число разбиений интервала интегрирования,
    h=(B-A)/N, причем N должно быть четным.
***********************************************************************)

Function Simps(b:extended;m:integer):extended;

   var
      sum, h: extended;
      N,j:Integer;
   begin
     N:=1000;
     {if Odd(N) then N:=N+1;}
     h:=(b-a)/N;
     sum:=0.5*(power(b,-1+m/2)*exp(-b/2));
     for j:=1 to N-1 do
       sum:=sum+(j mod 2+1)*power(j*h,-1+m/2)*exp(-j*h/2);
       Simps:=2*h*sum/3
   end;


(******************************************************
Процедура попарной обработки строк прямоугольной матрицы
******************************************************)

procedure Funky(Var Matrix_: Matrix; nn, mm: integer);
 var
  i,j,k,p,StSv: integer;
  Sum1,Sum2,Sum3: extended;

  dd:text;

 begin
 assign (dd,'matres.txt');
 rewrite(dd);
    for k:=1 to nn do
     begin
      for p:=1 to nn do
        begin
            Sum1:=0;
            Sum2:=0;
            Sum3:=0;
               for j:=1 to mm do
               begin
               StSv:=0;
               Sum1:=Sum1+Matrix_[k,j];
               Sum2:=Sum2+Matrix_[p,j];

               end;
               {writeln('Sum1=',Sum1);
               writeln('Sum2=',Sum2);
               readln;}

               for j:=1 to mm do
               begin
                 If ((Matrix_[k,j]+Matrix_[p,j])<>0) then
                   begin
                     If (Sum1=0) and (Sum2<>0) then

                     Sum3:=Sum3+((1/(Matrix_[k,j]+Matrix_[p,j]))*
                           (-Matrix_[p,j]/Sum2)*(-Matrix_[p,j]/Sum2));
                     If (Sum2=0) and (Sum1<>0) then
                     Sum3:=Sum3+( (1/(Matrix_[k,j]+Matrix_[p,j]))*
                           (Matrix_[k,j]/Sum1)*(Matrix_[k,j]/Sum1));

                     If Sum1*Sum2<>0 then
                        Sum3:=Sum3+((1/(Matrix_[k,j]+Matrix_[p,j]))*
                           (Matrix_[k,j]/Sum1-Matrix_[p,j]/Sum2)*
                           (Matrix_[k,j]/Sum1-Matrix_[p,j]/Sum2))

                   end
                 else
                     inc(StSv);
                     Sum3:=Sum3;
               end;
               {Funk[k,p]:=Sum1*Sum2*Sum3;}

               StSv:=mm-StSv;
               Funk[k,p]:= Simps(Sum1*Sum2*Sum3,StSv)/
                          (Power(2,StSv/2)*Gamma(StSv/2));
               write(dd,Funk[k,p]); write(dd,' ');

        end;
        writeln(dd);
     end;
close(dd);
 end;

(*******************************************************
Процедура замены в файле символа точки на символ запятой
*******************************************************)
 procedure ExchSymbolInFile (var String1,String2:string; var symbol_1,symbol_2:char);
 var
 char_:char;
 g,f:text;

begin
  assign (g,String1);
  reset(g);
  assign (f,String2);
  rewrite(f);

 while not eof(g) do
  begin

   while not eoln(g) do
         begin
             read(g,char_);
             if char_=symbol_1 then
                  write(f,symbol_2)
             else
              write (f,char_);
         end;
     readln(g);
     writeln(f);
  end;
  close(f);
  close(g);
end;


 begin

    file1:='matrix.txt';
    file2:='matrixdat.txt';
    s1:='.'; s2:=',';
    ExchSymbolInFile(file1,file2,s2,s1);

    EnterMatrixFromFile;

    Funky (Matr,n,m);

    file1:='mat_res.txt';
    file2:='mat_res_v.txt';
    ExchSymbolInFile(file1,file2,s1,s2);


 end.