svn/view/arraysc.pas

unit arraysc;
{$MODE DELPHI}{$H+}
//{$mode objfpc}{$H+}
//{$DEFINE MassArrays}
interface

uses
  Classes, SysUtils,
   base2, arrays, mathc,textfile;

type
    TComplexArray=array [0..1023] of Complex;
    PComplexArray=^TComplexArray;
type
   TArray2C=class;

   TArray1C=Class(TArray1)
      elements:PComplexArray;
      procedure subst(a:TArray1C);
      procedure add(a,b:TArray1C);
      procedure sub(a,b:TArray1C);
      procedure prod(a:TArray1C; b:TArray2C); overload;
      procedure prod(a:TArray2C; b:TArray1C); overload;
      procedure scalar(x:double; a:TArray1C); overload;
      procedure scalar(x:complex; a:TArray1C); overload;
      procedure con;overload;
      procedure con(x:double);overload;
      procedure con(x:complex);overload;
      procedure zer;overload;
      procedure zer(x:double);overload;
      procedure zer(x:complex);overload;
      procedure CROSS(a,b:TArray1C);
      destructor destroy;override;
      function InputDirective:string;
      function NewCopy:TArray1C;
      procedure MatPrint(ch:tTextDevice; direction:integer);override;
      procedure MatWrite(ch:tTextDevice);override;
      function kindlist:ansistring;
      procedure Read(list:TStringList;var p:integer);override; //list���p���������������������
      procedure AssignVarilen(list:TstringList);override;
      procedure LetWithTrace(ch:tTextDevice; name: ansistring;
                                             index1:complex; value: complex);
      function Array1N:TArray1N;
      function Array2N:TArray2N;   //������������������������
     private
      procedure init(lb1,ub1:integer);override;
      constructor createCopy(a:TArray1C);
      procedure CROSSsub(a,b:TArray1C);
    end;

   TArray2C=Class(TArray2)
      elements:PComplexArray;
      procedure subst(a:TArray2C);
      procedure add(a,b:TArray2C);
      procedure sub(a,b:TArray2C);
      procedure prod(a:TArray2C; b:TArray2C);
      procedure INV(a:TArray2C);
      procedure TRN(a:TArray2C);
      procedure scalar(x:double; a:TArray2C); overload;
      procedure scalar(x:complex;a:TArray2C); overload;
      procedure zer;overload;
      procedure zer(x:double);overload;
      procedure zer(x:complex);overload;
      procedure CON;overload;
      procedure CON(x:double);overload;
      procedure CON(x:complex);overload;
      procedure IDN;          overload;
      procedure IDN(x:double);overload;
      procedure IDN(x:complex);overload;
      destructor destroy;override;
      function InputDirective:string;
      function NewCopy:TArray2C;
      procedure MatPrint(ch:tTextDevice; direction:integer);override;
      procedure MatWrite(ch:tTextDevice);override;
      function kindlist:ansistring;
      procedure Read(list:TStringList;var p:integer);override; //list���p���������������������
      procedure LetWithTrace(ch:tTextDevice; name: ansistring;
                                             index1,index2:complex; value: complex);

      function Array2N:TArray2N;
     private
      procedure init(lb1,ub1,lb2,ub2:integer);override;
      constructor createCopy(a:TArray2C);
      procedure prodsub(a:TArray2C; b:TArray2C);
      procedure TRNSub(a:TArray2C);
      procedure determinant(var n:complex);
      function inverse:TArray2C;
    end;

   TArray3C=Class(TArray3)
      elements:PComplexArray;
      procedure subst(a:TArray3C);
      procedure add(a,b:TArray3C);
      procedure sub(a,b:TArray3C);
      procedure scalar(x:double; a:TArray3C);overload;
      procedure scalar(x:complex;a:TArray3C);overload;
      procedure zer;overload;
      procedure zer(x:double);overload;
      procedure zer(x:complex);overload;
      procedure CON;overload;
      procedure CON(x:double);overload;
      procedure CON(x:complex);overload;
      destructor destroy;override;
      function InputDirective:string;
      function NewCopy:TArray3C;
      procedure MatPrint(ch:tTextDevice; direction:integer);override;
      procedure MatWrite(ch:tTextDevice);override;
      function kindlist:ansistring;
      procedure Read(list:TStringList;var p:integer);override; //list���p���������������������
      procedure LetWithTrace(ch:tTextDevice; name: ansistring;
                                             index1,index2,index3:complex; value: complex);

     private
      procedure init(lb1,ub1,lb2,ub2,lb3, ub3:integer); override;
      constructor createCopy(a:TArray3C);

    end;

   TArray4C=Class(TArray4)
      elements:PComplexArray;
      procedure subst(a:TArray4C);
      procedure add(a,b:TArray4C);
      procedure sub(a,b:TArray4C);
      procedure scalar(x:double; a:TArray4C);overload;
      procedure scalar(x:complex;a:TArray4C);overload;
      procedure zer;overload;
      procedure zer(x:double);overload;
      procedure zer(x:complex);overload;
      procedure CON;overload;
      procedure CON(x:double);overload;
      procedure CON(x:complex);overload;
      destructor destroy;override;
      function InputDirective:string;
      function NewCopy:TArray4C;
      procedure MatPrint(ch:tTextDevice; direction:integer);override;
      procedure MatWrite(ch:tTextDevice);override;
      function kindlist:ansistring;
      procedure Read(list:TStringList;var p:integer);override; //list���p���������������������
      procedure LetWithTrace(ch:tTextDevice; name: ansistring;
                                   index1,index2,index3,index4:complex; value: complex);

     private
      procedure init(lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4:integer); override;
      constructor createCopy(a:TArray4C);

    end;


function dot(a,b:TArray1C):complex; overload;
function DET(a:TArray2C):complex;   overload;


implementation
uses base,vstack,float, baslibc;

procedure TArray1C.init(lb1,ub1:integer);
begin
    init0(lb1,ub1);
    if maxsize>0 then
      {$IFDEF MassArrays}
       Elements:=AllocMem(Maxsize*SizeOf(Complex));
      {$ELSE}
       Elements:=GetZeroMemory(Maxsize*SizeOf(Complex));
      {$ENDIF}
end;

destructor TArray1C.destroy;
begin
  {$IFDEF MassArrays}
    FreeMem(Elements,Maxsize*SizeOf(Complex));
  {$ELSE}
    FreeMemory(Maxsize*SizeOf(Complex));
  {$ENDIF}
    inherited destroy;
end;


constructor TArray1C.createCopy(a:TArray1C);
var
  i:integer;
begin
   TArray.create;
   Lbound1:=a.Lbound1;
   Size1:=a.Size1;
   Maxsize:=Size1;
   if maxsize>0 then
       begin
         Elements:=GetZeroMemory(Maxsize*SizeOf(Complex));
         for i:=0 to maxsize-1 do
             elements^[i]:=a.elements^[i];
       end;
end;

function TArray1C.NewCopy:TArray1C;
begin
   result:=TArray1C.createCopy(self)
end;


procedure TArray2C.init(lb1,ub1,lb2,ub2:integer);
begin
    init0(lb1,ub1,lb2,ub2);
    if maxsize>0 then
    {$IFDEF MassArrays}
       Elements:=AllocMem(Maxsize*SizeOf(Complex));
    {$ELSE}
       Elements:=GetZeroMemory(Maxsize*SizeOf(Complex));
    {$ENDIF}
end;

destructor TArray2C.destroy;
begin
   {$IFDEF MassArrays}
    FreeMem(Elements,Maxsize*SizeOf(Complex));
   {$ELSE}
    FreeMemory(Maxsize*SizeOf(Complex));
   {$ENDIF}
    inherited destroy;
end;

constructor TArray2C.createCopy(a:TArray2C);
var
i:integer;
begin
   TArray.create;
   Lbound1:=a.Lbound1;
   Size1:=a.Size1;
   Lbound2:=a.Lbound2;
   Size2:=a.Size2;
   Maxsize:=Size1*size2;
   if maxsize>0 then
       begin
         Elements:=GetZeroMemory(Maxsize*SizeOf(Complex));
         for i:=0 to maxsize-1 do
             elements^[i]:=a.elements^[i];
       end;
end;

function TArray2C.NewCopy:TArray2C;
begin
   result:=TArray2C.createCopy(self)
end;

procedure TArray3C.init(lb1,ub1,lb2,ub2,lb3, ub3:integer);
begin
    init0(lb1,ub1,lb2,ub2,lb3,ub3);
    if maxsize>0 then
    {$IFDEF MassArrays}
       Elements:=AllocMem(Maxsize*SizeOf(Complex));
    {$ELSE}
       Elements:=GetZeroMemory(Maxsize*SizeOf(Complex));
    {$ENDIF}
end;

destructor TArray3C.destroy;
begin
   {$IFDEF MassArrays}
     FreeMem(Elements,Maxsize*SizeOf(Complex));
   {$ELSE}
     FreeMemory(Maxsize*SizeOf(Complex));
   {$ENDIF}
    inherited destroy;
end;
 constructor TArray3C.createCopy(a:TArray3C);
var
  i:integer;
begin
   TArray.create;
   Lbound1:=a.Lbound1;
   Size1:=a.Size1;
   Lbound2:=a.Lbound2;
   Size2:=a.Size2;
   Lbound3:=a.Lbound3;
   Size3:=a.Size3;
   Maxsize:=Size1*size2*size3;
   if maxsize>0 then
       begin
         Elements:=GetZeroMemory(Maxsize*SizeOf(Complex));
         for i:=0 to maxsize-1 do
             elements^[i]:=a.elements^[i];
       end;
end;

function TArray3C.NewCopy:TArray3C;
begin
   result:=TArray3C.createCopy(self)
end;


 procedure TArray4C.init(lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4:integer);
 begin
     init0(lb1,ub1,lb2,ub2,lb3,ub3,lb4,ub4);
     if maxsize>0 then
     {$IFDEF MassArrays}
        Elements:=AllocMem(Maxsize*SizeOf(Complex));
     {$ELSE}
        Elements:=GetZeroMemory(Maxsize*SizeOf(Complex));
     {$ENDIF}
 end;

 destructor TArray4C.destroy;
 begin
    {$IFDEF MassArrays}
      FreeMem(Elements,Maxsize*SizeOf(Complex));
    {$ELSE}
      FreeMemory(Maxsize*SizeOf(Complex));
    {$ENDIF}
     inherited destroy;
 end;

 constructor TArray4C.createCopy(a:TArray4C);
 var
   i:integer;
 begin
    TArray.create;
    Lbound1:=a.Lbound1;
    Size1:=a.Size1;
    Lbound2:=a.Lbound2;
    Size2:=a.Size2;
    Lbound3:=a.Lbound3;
    Size3:=a.Size3;
    Lbound4:=a.Lbound3;
    Size4:=a.Size4;
    Maxsize:=Size1*size2*size3*size4;
    if maxsize>0 then
        begin
          Elements:=GetZeroMemory(Maxsize*SizeOf(Complex));
          for i:=0 to maxsize-1 do
              elements^[i]:=a.elements^[i];
        end;
 end;

 function TArray4C.NewCopy:TArray4C;
 begin
    result:=TArray4C.createCopy(self)
 end;


{**************}
{MAT OPERATIONS}
{**************}

procedure TArray1C.subst(a:TArray1C);
var
   i:integer;
begin
   if a.size1>maxsize then setexception(5001);
   Size1:=a.Size1;
   for i:=0 to a.size1-1 do
       elements^[i]:=a.elements^[i]

end;

procedure TArray1C.add(a,b:TArray1C);
var
   c:TArray1C;
   i:integer;
begin
   if (a.Size1<>b.Size1) then
      setexception(6001);
   if a.Size>MaxSize then
      setexception(5001);
   c:=TArray1C.create(1,a.size1);
   try
     for i:=0 to c.size-1 do
         c.elements^[i]:=a.elements^[i]+b.elements^[i];
     subst(c)
   finally
     c.free;
   end;
end;

procedure TArray1C.sub(a,b:TArray1C);
var
   c:TArray1C;
   i:integer;
begin
   if (a.Size1<>b.Size1) then
      setexception(6001);
   if a.Size>MaxSize then
      setexception(5001);
   c:=TArray1C.create(1,a.size1);
   try
     for i:=0 to c.size-1 do
         c.elements^[i]:=a.elements^[i]-b.elements^[i];
     subst(c)
   finally
     c.free;
   end;
end;

procedure TArray1C.prod(a:TArray1C; b:TArray2C); overload;
var
   i,j:integer;
   c:TArray1C;
begin
   if a.size1<>b.size1 then setexception(6001);
   if maxsize<b.size2 then setexception(5001);
   c:=TArray1C.create(1, b.size2);
   //c.zer;
   try
     for j:=0 to b.size2-1 do
       for i:=0 to b.size1 -1 do
         c.elements^[j]:=c.elements^[j]+ a.elements^[i]*b.elements^[i*b.size2+j];
     subst(c);
   finally
     c.free;
   end;

end;

procedure TArray1C.prod(a:TArray2C; b:TArray1C); overload;
var
   i,j:integer;
   c:TArray1C;
begin
   if a.size2<>b.size1 then setexception(6001);
   if maxsize<a.size1 then setexception(5001);
   c:=TArray1C.create(1, a.size1);
   //c.zer;
   try
     for i:=0 to a.size1-1 do
       for j:=0 to a.size2 -1 do
         c.elements^[i]:=c.elements^[i]+ a.elements^[i*a.size2+j]*b.elements^[j];
     subst(c);
   finally
     c.free;
   end;
end;

procedure TArray1C.scalar(x:double; a:TArray1C);
var
  i:integer;
begin
  subst(a);
  for i:=0 to size-1 do
     elements^[i]:=x*elements^[i]
end;

procedure TArray1C.scalar(x:complex; a:TArray1C);
var
  i:integer;
begin
  subst(a);
  for i:=0 to size-1 do
     elements^[i]:=x*elements^[i]
end;

procedure TArray1C.CON(x:double);
var
  i:integer;
begin
  for i:=0 to size-1 do
     elements^[i]:=x
end;

procedure TArray1C.CON(x:complex);
var
  i:integer;
begin
  for i:=0 to size-1 do
     elements^[i]:=x
end;

procedure TArray1C.CON;
begin
  CON(1)
end;

procedure TArray1C.zer;overload;
var
   i:integer;
begin
   for i:=0 to size1-1 do
      elements^[i]:=0;
end;

procedure TArray1C.zer(x:double);overload;
begin
   zer
end;

procedure TArray1C.zer(x:complex);overload;
begin
   zer
end;

procedure TArray1C.CROSSsub(a,b:TArray1C);
var
   i:integer;
   x,y:complex;
begin
   for i:=0 to 2 do
     begin
       elements^[i mod 3]:=0;
       x:=a.elements^[(i+1) mod 3];
       x:=x*b.elements^[(i+2) mod 3];
       y:=b.elements^[(i+1) mod 3];
       y:=y*a.elements^[(i+2) mod 3];
       elements^[i mod 3]:=elements^[i mod 3]+x-y;
    end;
end;

procedure TArray1C.CROSS(a,b:TArray1C);
var
  c:TArray1C;
begin
   if (a.size<3) or (b.size<3) then setexception(6001);
   if MaxSize<3 then setexception(5001);
   c:=TArray1C.create(1,3);
   try
     c.CrossSub(a,b);
     subst(c);
   finally
     c.free;
   end;
end;


procedure TArray2C.subst(a:TArray2C);
var
   i:integer;
begin
   if a.Size>MaxSize then
      setexception(5001);
   resize(a.size1,a.size2);
   for i:=0 to size-1 do
          elements^[i]:=A.elements^[i]

end;

procedure TArray2C.add(a,b:TArray2C);
var
   c:TArray2C;
   i:integer;
begin
   if (a.Size1<>b.Size1) or (a.Size2<>b.Size2) then
      setexception(6001);
   if a.Size>MaxSize then
      setexception(5001);
   c:=TArray2C.create(1, a.size1, 1, a.size2);
   try
     for i:=0 to c.size-1 do
         c.elements^[i]:=a.elements^[i]+b.elements^[i];
     subst(c)
   finally
     c.free;
   end;
end;


procedure TArray2C.sub(a,b:TArray2C);
var
   c:TArray2C;
   i:integer;
begin
   if (a.Size1<>b.Size1) or (a.Size2<>b.Size2) then
      setexception(6001);
   if a.Size>MaxSize then
      setexception(5001);
   c:=TArray2C.create(1, a.size1, 1, a.size2);
   try
     for i:=0 to c.size-1 do
         c.elements^[i]:=a.elements^[i]-b.elements^[i];
     subst(c)
   finally
     c.free;
   end;
end;

procedure TArray2C.ProdSub(a:TArray2C; b:TArray2C);
var
  i,j,k:integer;
begin
   for i:=0 to size1-1 do
      for j:=0 to size2-1 do
         for k:=0 to a.size2-1 do
            elements^[size2*i+j]:=elements^[Size2*i+j]+a.elements^[a.Size2*i+k]*b.elements^[b.size2*k+j];
end;
procedure TArray2C.prod(a:TArray2C; b:TArray2C);
var
  c:TArray2C;
  //i,j,k:integer;
begin
   if (a.Size2<>b.Size1)  then
      setexception(6001);
   if a.Size1 * b.Size2 >MaxSize then
      setexception(5001);
   c:=TArray2C.create(1, a.size1, 1, b.size2);
   try
     c.ProdSub(a,b);
     subst(c);
   finally
     c.free;
   end;
end;

procedure TArray2C.INV(a:TArray2C);
var
   c:TArray2C;
begin
   c:=a.inverse;
   subst(c);
   c.free;
end;

procedure TArray2C.TRNSub(a:TArray2C);
var
   i,j:integer;
begin
   for i:=0 to size1 -1 do
     for j:=0 to size2 -1 do
        elements^[j+size2*i]:=a.elements^[i+a.Size2*j];
end;

procedure TArray2C.TRN(a:TArray2C);
var
   c:TArray2C;
begin
   c:=TArray2C.create(1, a.Size2, 1, a.size1);
   c.TRNSub(a);
   subst(c)
end;

procedure TArray2C.scalar(x:double; a:TArray2C);
var
  i:integer;
begin
  subst(a);
  for i:=0 to size-1 do
     elements^[i]:=x*elements^[i]
end;

procedure TArray2C.scalar(x:complex; a:TArray2C);
var
  i:integer;
begin
  subst(a);
  for i:=0 to size-1 do
     elements^[i]:=x*elements^[i]
end;

procedure TArray2C.con;
begin
  CON(1)
end;

procedure TArray2C.con(x:double);
var
  i:integer;
begin
  for i:=0 to size-1 do
     elements^[i]:=x
end;

procedure TArray2C.con(x:complex);
var
  i:integer;
begin
  for i:=0 to size-1 do
     elements^[i]:=x
end;

procedure TArray2C.zer;overload;
var
   i:integer;
begin
   for i:=0 to size-1 do
      elements^[i]:=0;
end;

procedure TArray2C.zer(x:double);overload;
begin
   zer
end;

procedure TArray2C.zer(x:complex);overload;
begin
   zer
end;


procedure TArray2C.IDN;          overload;
begin
   IDN(1)
end;

procedure TArray2C.IDN(x:double);overload;
var
   i:integer;
begin
   if size1<>size2 then setexception(6004);
   zer;
   for i:=0 to size1-1 do
       elements^[i+size2*i]:=x
end;

procedure TArray2C.IDN(x:complex);overload;
var
   i:integer;
begin
   if size1<>size2 then setexception(6004);
   zer;
   for i:=0 to size1-1 do
       elements^[i+size2*i]:=x
end;


procedure TArray3C.subst(a:TArray3C);
var
   i:integer;
begin
   if a.Size>MaxSize then
      setexception(5001);
   resize(a.size1,a.size2,a.size3);
   for i:=0 to size-1 do
          elements^[i]:=a.elements^[i]

end;

procedure TArray3C.add(a,b:TArray3C);
var
   c:TArray3C;
   i:integer;
begin
   if (a.Size1<>b.Size1) or (a.Size2<>b.Size2) or (a.Size3<>b.Size3)then
      setexception(6001);
   if a.Size>MaxSize then
      setexception(5001);
   c:=TArray3C.create(1, a.size1, 1, a.size2, 1, a.size3);
   try
     for i:=0 to c.size-1 do
         c.elements^[i]:=a.elements^[i]+b.elements^[i];
     subst(c)
   finally
     c.free;
   end;
end;


procedure TArray3C.sub(a,b:TArray3C);
var
   c:TArray3C;
   i:integer;
begin
   if (a.Size1<>b.Size1) or (a.Size2<>b.Size2) or (a.Size3<>b.Size3)then
      setexception(6001);
   if a.Size>MaxSize then
      setexception(5001);
   c:=TArray3C.create(1, a.size1, 1, a.size2, 1, a.size3);
   try
     for i:=0 to c.size-1 do
         c.elements^[i]:=a.elements^[i]-b.elements^[i];
     subst(c)
   finally
     c.free;
   end;
end;


procedure TArray3C.scalar(x:double; a:TArray3C);
var
  i:integer;
begin
  subst(a);
  for i:=0 to size-1 do
     elements^[i]:=x*elements^[i]
end;

procedure TArray3C.scalar(x:complex; a:TArray3C);
var
  i:integer;
begin
  subst(a);
  for i:=0 to size-1 do
     elements^[i]:=x*elements^[i]
end;

procedure TArray3C.con;
begin
  CON(1)
end;

procedure TArray3C.con(x:double);
var
  i:integer;
begin
  for i:=0 to size-1 do
     elements^[i]:=x
end;

procedure TArray3C.con(x:complex);
var
  i:integer;
begin
  for i:=0 to size-1 do
     elements^[i]:=x
end;

procedure TArray3C.zer;overload;
var
   i:integer;
begin
   for i:=0 to size-1 do
      elements^[i]:=0;
end;

procedure TArray3C.zer(x:double);overload;
begin
   zer
end;

procedure TArray3C.zer(x:complex);overload;
begin
   zer
end;


procedure TArray4C.subst(a:TArray4C);
var
   i:integer;
begin
   if a.Size>MaxSize then
      setexception(5001);
   resize(a.size1,a.size2,a.size3,a.size4);
   for i:=0 to size-1 do
          elements^[i]:=a.elements^[i]

end;

procedure TArray4C.add(a,b:TArray4C);
var
   c:TArray4C;
   i:integer;
begin
   if (a.Size1<>b.Size1) or (a.Size2<>b.Size2) or (a.Size3<>b.Size3)
                                               or (a.Size4<>b.Size4)then
      setexception(6001);
   if a.Size>MaxSize then
      setexception(5001);
   c:=TArray4C.create(1, a.size1, 1, a.size2, 1, a.size3, 1, a.size4);
   try
     for i:=0 to c.size-1 do
         c.elements^[i]:=a.elements^[i]+b.elements^[i];
     subst(c)
   finally
     c.free;
   end;
end;


procedure TArray4C.sub(a,b:TArray4C);
var
   c:TArray4C;
   i:integer;
begin
   if (a.Size1<>b.Size1) or (a.Size2<>b.Size2) or (a.Size3<>b.Size3)
                                               or (a.Size4<>b.Size4)then
      setexception(6001);
   if a.Size>MaxSize then
      setexception(5001);
   c:=TArray4C.create(1, a.size1, 1, a.size2, 1, a.size3, 1, a.size4);
   try
     for i:=0 to c.size-1 do
         c.elements^[i]:=a.elements^[i]-b.elements^[i];
     subst(c)
   finally
     c.free;
   end;
end;


procedure TArray4C.scalar(x:double; a:TArray4C);
var
  i:integer;
begin
  subst(a);
  for i:=0 to size-1 do
     elements^[i]:=x*elements^[i]
end;

procedure TArray4C.scalar(x:complex; a:TArray4C);
var
  i:integer;
begin
  subst(a);
  for i:=0 to size-1 do
     elements^[i]:=x*elements^[i]
end;

procedure TArray4C.con;
begin
  CON(1)
end;

procedure TArray4C.con(x:double);
var
  i:integer;
begin
  for i:=0 to size-1 do
     elements^[i]:=x
end;

procedure TArray4C.con(x:complex);
var
  i:integer;
begin
  for i:=0 to size-1 do
     elements^[i]:=x
end;

procedure TArray4C.zer;overload;
var
   i:integer;
begin
   for i:=0 to size-1 do
      elements^[i]:=0;
end;

procedure TArray4C.zer(x:double);overload;
begin
   zer
end;

procedure TArray4C.zer(x:complex);overload;
begin
   zer
end;

{*******}
{MAT I/O}
{*******}

procedure TArray1C.MatPrint(ch:tTextDevice; direction:integer);
var
   i:integer;
begin
   ch.newlineifneed;
   for i:=0 to size1 -1 do
     begin
        if direction<>0 then ch.NewZone;
        ch.AppendStrV2(CStr(elements^[i])+' ')
     end;
   ch.newline;
   ch.newline;
end;

procedure TArray2C.MatPrint(ch:tTextDevice; direction:integer);
var
   i,j:integer;
begin
   ch.newlineifneed;
   for i:=0 to size1 -1 do
     begin
       for j:=0 to size2-1 do
         begin
          if direction<>0 then ch.NewZone;
          ch.AppendStrV2(CStr(elements^[i*size2 + j])+' ')
         end;
       ch.newline;
     end;
   ch.newline
end;

procedure TArray3C.MatPrint(ch:tTextDevice; direction:integer);
var
   i,j,k:integer;
begin
   ch.newlineifneed;
   for i:=0 to size1 -1 do
     begin
       for j:=0 to size2-1 do
         begin
           for k:=0 to size3-1 do
              begin
                if direction<>0 then ch.NewZone;
                ch.AppendStrV2(CSTR(elements^[(i*size2 + j)*size3 + k ])+' ')
              end;
            ch.newline;
         end;
       ch.newline;
     end;
end;

procedure TArray4C.MatPrint(ch:tTextDevice; direction:integer);
var
   i,j,k,l:integer;
begin
   ch.newlineifneed;
   for i:=0 to size1 -1 do
     begin
       for j:=0 to size2-1 do
         begin
           for k:=0 to size3-1 do
              begin
                for l:=0 to size4-1 do
                  begin
                    if direction<>0 then ch.NewZone;
                    ch.AppendStrV2(CSTR(elements^[((i*size2 + j)*size3 + k)*size4 +l ])+' ')
                  end;
                ch.newline;
              end;
            ch.newline;
         end;
       ch.newline;
     end;
end;

procedure TArray1C.MatWrite(ch:tTextDevice);
var
   i:integer;
begin
   for i:=0 to size1 -1 do
     begin
        if i>0 then
           ch.WriteSeparator(false);
        ch.AppendStrV2(elements^[i])
     end;
end;

procedure TArray2C.MatWrite(ch:tTextDevice);
var
   i:integer;
begin
   for i:=0 to size -1 do
     begin
        if i>0 then
           ch.WriteSeparator((ch.RecType=rcCSV) and (i mod size1=0));
        ch.AppendStrV2(elements^[i])
     end;
end;

procedure TArray3C.MatWrite(ch:tTextDevice);
var
   i:integer;
begin
   for i:=0 to size -1 do
     begin
        if i>0 then
           ch.WriteSeparator((ch.RecType=rcCSV) and (i mod size1=0));
        ch.AppendStrV2(elements^[i])
     end;
end;

procedure TArray4C.MatWrite(ch:tTextDevice);
var
   i:integer;
begin
   for i:=0 to size -1 do
     begin
        if i>0 then
           ch.WriteSeparator((ch.RecType=rcCSV) and (i mod size1=0));
        ch.AppendStrV2(elements^[i])
     end;
end;

procedure TArray1C.Read(list:TStringList;var p:integer); //list���p���������������������
var
  i:integer;
begin
  for i:=0 to size-1 do
    begin
      elements^[i]:=StrToFloat(list.Strings[p]);
      inc(p)
    end;
end;

procedure TArray2C.Read(list:TStringList;var p:integer); //list���p���������������������
var
  i:integer;
begin
  for i:=0 to size-1 do
    begin
      elements^[i]:=StrToFloat(list.Strings[p]);
      inc(p)
    end;
end;

procedure TArray3C.Read(list:TStringList;var p:integer); //list���p���������������������
var
  i:integer;
begin
  for i:=0 to size-1 do
    begin
      elements^[i]:=StrToFloat(list.Strings[p]);
      inc(p)
    end;
end;

procedure TArray4C.Read(list:TStringList;var p:integer); //list���p���������������������
var
  i:integer;
begin
  for i:=0 to size-1 do
    begin
      elements^[i]:=StrToFloat(list.Strings[p]);
      inc(p)
    end;
end;

function TArray1C.kindlist:ansistring;
begin
  result:=StringOfChar('n',size);
end;

function TArray2C.kindlist:ansistring;
begin
  result:=StringOfChar('n',size);
end;

function TArray3C.kindlist:ansistring;
begin
  result:=StringOfChar('n',size);
end;

function TArray4C.kindlist:ansistring;
begin
  result:=StringOfChar('n',size);
end;

function TArray1C.InputDirective:string;
begin
    result:=StringOfChar('n',size1);
end;

function TArray2C.InputDirective:string;
begin
    result:=StringOfChar('n',size1*Size2);
end;

function TArray3C.InputDirective:string;
begin
    result:=StringOfChar('n',size1*Size2*Size3);
end;

function TArray4C.InputDirective:string;
begin
    result:=StringOfChar('n',size1*Size2*Size3*size4);
end;

procedure TArray1C.AssignVarilen(list:TstringList);
var
  i:integer;
begin
  ReSize(list.count);
  with list do
    for i:=0 to Count-1 do
      elements^[i]:=StrToFloat(Strings[i]);
end;

{**********}
{ DOT & DET}
{**********}

function dot(a,b:TArray1C):complex;
var
   i:integer;
begin
   if a.size1<>b.size1 then setexception(6001);
   result:=0;
   try
     for i:=0 to a.size-1 do
       result:=result+a.elements^[i]*b.elements^[i]
   except
      on EOverFlow do
      begin
        {$IFNDEF Windows} asm finit end; set8087cw(controlword); {$ENDIF}
         raise EExtype.create(1009)
      end;
   end;
end;

function DET(a:TArray2C):complex;
begin
   a.Determinant(result)
end;

{**************}
{INVERSE MATRIX}
{**************}


procedure matinv(size:integer; p:PcomplexArray;  pv:PIntArray; var det:complex);
{$MAXFPUREGISTERS 0}
  function a(i,j:integer):Pcomplex;
  begin
     result:=@p^[i+j*size]
  end;
var
  i,j,k,tmp:integer;
  t,u,temp:complex;
  eps:double;
label
  EXIT;
begin
  eps:=1; FEPS(eps); eps:=eps/2;
  for k:=0 to size-1 do pv^[k]:=k;
  det:=1;
  for k:=0 to size-1 do
     begin
        i:=k;
        while  (i<size) and (a(i,k)^=0.0) do inc(i);
        if i=size then
           begin det:=0.0; goto EXIT end
        else if i<>k then
           begin
              tmp:=pv^[i]; pv^[i]:=pv^[k]; pv^[k]:=tmp;
              for j:=0 to size-1 do
                  begin  temp:=a(i,j)^; a(i,j)^:=a(k,j)^; a(k,j)^:=temp end;
              det:=-det;
           end;

        t:=a(k,k)^;
        det:=det*t;
        for i:=0 to size-1 do
            a(k,i)^:=a(k,i)^/t;
        a(k,k)^:=1.0/t;
        for j:=0 to size-1 do
          if j<>k then
           begin
            u:=a(j,k)^;
            for i:=0 to k-1 do
                 begin
                   a(j,i)^:=a(j,i)^-a(k,i)^*u;
                 end;
            a(j,k)^:=-u/t;
            for i:=k+1 to size-1 do
                 begin
                   temp:=a(j,i)^-a(k,i)^*u;
                   if mathc.abs(temp)<mathc.abs(a(j,i)^)*EPS then temp:=0;
                   a(j,i)^:=temp;
                 end
           end;
        idle;
     end;
  EXIT:
end;

procedure TArray2C.determinant(var n:complex);
{$MAXFPUREGISTERS 0}
var
  i,j:integer;
  det:complex;
  a:PcomplexArray;
  pv:PIntArray;
begin
  if size1=size2 then
    begin
      getmem(pv,size1*sizeof(integer));
      getmem(a,size1*size2*sizeof(complex));
         try
            for i:=0 to size1-1 do
              for j:=0 to size2-1 do
                a^[i+j*size1]:=elements^[i*size2+j];
            matinv(size1,a,pv,det);
            n:=det;
         finally
            freemem(a,size1*size2*sizeof(complex));
            freemem(pv,size1*sizeof(integer));
         end
    end
  else
     setexception(6002);
end;

function TArray2C.inverse:TArray2C;
var
  i,j:integer;
  det:complex;
  p:PcomplexArray;
  pv:PIntArray;
begin
   result:=nil;
   if size1=size2 then
   begin
     getmem(pv,size1*sizeof(integer)+size1*size2*sizeof(complex));
     try
        try
            p:=@pv^[size1];
            for i:=0 to size1-1 do
              for j:=0 to size2-1 do
                p^[i+j*size1]:=elements^[i*size2+j];
            matinv(size1,p,pv,det);
            if det=0 then
               setexception(3009)
            else
              begin
                result:=NewCopy;
                if result=nil then
                   setexception(ArraySizeOverflow)
                else
                begin
                  try
                    for i:=0 to size1-1 do
                      for j:=0 to size2-1 do
                        result.elements^[i*size2+pv^[j]]:=p^[i+j*size1];
                  except
                    on EMathError do
                       begin
                      {$IFNDEF Windows} asm finit end; set8087cw(controlword); {$ENDIF}
                       end;
                    on EDivByZero do
                       begin
                      {$IFNDEF Windows} asm finit end; set8087cw(controlword); {$ENDIF}
                       end;
                  end;
                end;
              end;
        finally
            freemem(pv,size1*sizeof(integer)+size1*size2*sizeof(complex));
        end
     except
           result.free;
           result:=nil;
           raise;
     end;
   end
  else
     setexception(6003);
end;

function TArray1C.Array1N:TArray1N;
var
   i:integer;
begin
   result:=TArray1N.create(Lbound,Ubound);
   try
     for i:=0 to size1 -1 do
       result.elements^[i]:=testreal(elements^[i]);
   except
     result.free;
     result:=nil;
     raise;
   end;
end;

function TArray1C.Array2N:TArray2N;
var
   i:integer;
begin
   result:=TArray2N.create(Lbound,Ubound, 0,1);
   try
     for i:=0 to size1 -1 do
       begin
          result.elements^[i*2]:=(elements^[i]).x;
          result.elements^[i*2+1]:=(elements^[i]).y;
       end;
   except
     result.free;
     result:=nil;
     raise;
   end;
end;


function TArray2C.Array2N:TArray2N;
var
   i:integer;
begin
   result:=TArray2N.create(Lbound(1),Ubound(1), Lbound(2), Ubound(2));
   try
     for i:=0 to size - 1 do
       result.elements^[i]:=testreal(elements^[i]);
   except
     result.free;
     result:=nil;
     raise;
   end;
end;

procedure TArray1C.LetWithTrace(ch:tTextDevice; name: ansistring; index1:complex; value: complex);
begin
   elements^[index(index1)]:=value;
   ch.PRINT([],rsNone, false ,[' '+name+'('+Format17(LongintRound(testreal(index1)))+')=',
                                                 TComplex.create(value), TNewLine.create]);
end;

procedure TArray2C.LetWithTrace(ch:tTextDevice; name: ansistring; index1,index2:complex; value: complex);
begin
   elements^[index(index1,index2)]:=value;
   ch.PRINT([],rsNone, false ,[' '+name+'('+Format17(LongintRound(testreal(index1)))+ ','+
                                            Format17(LongintRound(testreal(index2)))+ ')=',
                                                TComplex.create(value), TNewLine.create]);
end;

procedure TArray3C.LetWithTrace(ch:tTextDevice; name: ansistring; index1,index2,index3:complex; value: complex);
begin
   elements^[index(index1,index2,index3)]:=value;
   ch.PRINT([],rsNone, false ,[' '+name+'('+Format17(LongintRound(testreal(index1)))+ ','+
                                            Format17(LongintRound(testreal(index2)))+ ','+
                                            Format17(LongintRound(testreal(index3)))+ ')=',
                                                 TComplex.create(value), TNewLine.create]);
end;

procedure TArray4C.LetWithTrace(ch:tTextDevice; name: ansistring; index1,index2,index3,index4:complex; value: complex);
begin
   elements^[index(index1,index2,index3,index4)]:=value;
   ch.PRINT([],rsNone, false ,[' '+name+'('+Format17(LongintRound(testreal(index1)))+ ','+
                                            Format17(LongintRound(testreal(index2)))+ ','+
                                            Format17(LongintRound(testreal(index3)))+ ','+
                                            Format17(LongintRound(testreal(index4)))+ ')=',
                                                 TComplex.create(value), TNewLine.create]);
end;

end.
