// Unit MATRICE

// OBJECT : matrix and vector aritmetic (inversion, addition, multiplication ...)
// first written in Turbo-Pascal (02-05-1991)
// adapted to Delphi 4 (28-05-1999)

// NOTES

// Note 1 : these routines operate on dynamical square matrix and dynamical vectors
// (RCOMat and RCOVec),whose dimension should have been be declared with
// the InitMat(Dim) procedure (error if if not initialize).

// Note 2 : the index (i,j..) are "zero based", this is a constraint of
// Delphi4 dynamical matrix structure.

// Note 3 : the "if true ..." instructions are remainings of the TurboPascal 5 version,
// I should have transform them into "try .. except .." instructions in order to
// manage memory overflow errors. If somebody do that, I can mail me the upgraded version.

// Note 4 : I checked this unit with the small program at the end of that file.
// The time required to invert a 200 x 200 extented matrix is of about 2 seconds,
// (AMD K6-2 overclocked at 400 MHz), and the accuracy compatible with extended
// accuracy (10E-17 .. 10E-18)

// if you find that unit usefull, please let me know :
// remi CORMIER, rco@club-internet.fr

UNIT Matrice;

//determinant为行列式的值

INTERFACE
uses SysUtils,Dialogs,PerlRegEx,classes;
// I would'nt worry if you change "RCOMat" to "YourNameMat" ...
TYPE      RCOMat = array of array of extended;
          RCOVec = array of extended;

VAR       ErreurMat : integer; { code d'erreur, Cf ci apr鑣 :   }

CONST     MatCorrect   = 0;
          MatNonInit   = 1;
          MatDebordTas = 2;
          MatDetNul    = 3;
//两个行列相同的矩阵相加
function MatAddMat(A,B:RCOMat;i,j:integer):RCOMat;

//将矩阵转成相反数，求A的负
function FuMat(A:RCOMat;i,j:integer):RCOMat;

//小数转分数
function Fraction(decimal: double): string;

//从文本文件中读取矩阵到A ，i,j为返回的矩阵长度
procedure MatFromText(var A:RCOMat;path:string;var i,j:integer);

//用MessageBox显示出矩阵，isT表示是否用分数显示，默认为用小数显示
procedure MatToMessageBox(A:RCOMat;i,j:integer;isT:boolean=false;isToText:boolean=true);

// Returns (W) = (U)-(V) ; U,V,W : Vectors
PROCEDURE DifVec(U,V:RCOVec;var W:RCOVec);

// Returns the current dimension of matrix and vectors
FUNCTION DimMatrices:integer;

// Initialize the matrix and vectors dimension with Dim
PROCEDURE InitMat(Dim:integer);

// Returns [B], the invert matrix of [A] ([A]x[B]=[I]), Det is the determinant of [A] (and [B])
//求A的逆阵B   只有方阵才有逆 返回det为A的行列式的值
PROCEDURE InvMat(A:RCOMat;var B:RCOMat;var Det:extended{MatElem});

// Returns [C] = [A]x[B]
PROCEDURE MatMat(A,B:RCOMat;var C:RCOMat);

// Returns [C] = [A]x[B] ; same as MatMat, except that [A] and [B] can be
// rectangular. [A] have Dim rows and M columns, [B] have M rows and P colums :
// [A]N,M x [B]M,P -> [C]N,P (with N = DimMat)
PROCEDURE MatMatR(A,B:RCOMat;M,P:integer;var C:RCOMat);

// Returns the unit matrix (one on diagonal, zero everywhere else)返回单元矩阵
PROCEDURE MatUnt(var MU:RCOMat);

// Returns the (W)=[M]x(V) ; [M] RCOMat, (V) and (W) RCOVec
PROCEDURE MatVec(M:RCOMat;V:RCOVec;var W:RCOVec);

// Returns the norm of vector (U)
FUNCTION Norme(U:RCOVec)    : extended;

// Returns the norm of [A], [A] can be rectangular, with n row and m columns
FUNCTION NorMat(A:RCOMat;n,m:integer):extended;

// Returns the scalar product of the vector (U) and (V)
FUNCTION PrdScal(U,V:RCOVec) : extended;

// Multiply the matrix [A] with x, the result in matrix [B]
PROCEDURE ScalMat(x:extended;A:RCOMat;var B:RCOMat);

// Multiply the vector (U) with x, the result in vector (V)
PROCEDURE ScalVec(x:extended;U:RCOVec;var V:RCOVec);

// Addition of two square matrix [A] and [B], result in [C]
PROCEDURE SomMat(A,B:RCOMat;var C:RCOMat);

// Addition of two vectors, (U) and (V), result in (W)
PROCEDURE SomVec(U,V:RCOVec;var W:RCOVec);

// Solution of the linear system [A]x(X)=(B)
// vector (X) is the solution
// check Det for valid solution (should be <> 0)
PROCEDURE SysLin(A:RCOMat;var X:RCOVec;B:RCOVec;var Det:extended{MatElem});

// Returns the transposed matrix of [A] into [B] (Bij = Aji)   转置行列式
PROCEDURE TransMat(A:RCOMat;var B:RCOMat;n,m:integer); //n,m为A的行、列数

// Returns the text (in french) associated with the error # no
FUNCTION MsgErreurMat(no:integer):string;

IMPLEMENTATION

VAR   Ini       : boolean;
      DimMat    : integer;

TYPE RCoVecE   = array of integer;

function MatAddMat(A,B:RCOMat;i,j:integer):RCOMat;
var C:RCOMat;
    m,n:integer;
begin
    InitMat(i);
    SetLength(C,i,j);
    for m := 0 to i - 1 do
      for n := 0 to j - 1 do
        C[m,n] := A[m,n]+B[m,n];
    result := C;
end;

function FuMat(A:RCOMat;i,j:integer):RCOMat;
var B:RCOMat;
    temp,temp1:integer;
begin
   setlength(B,i,j);
   for temp := 0 to i - 1 do
     for temp1 := 0 to j - 1 do
       B[temp,temp1] := -A[temp,temp1];
   result := B;
end;

function Fraction(decimal: double): string;
var
intNumerator, intDenominator, intNegative: integer; // 声明整数变量为长整数
dblFraction, dblDecimal, dblAccuracy, dblinteger: Double; // 声明浮点数为双精度
temp:string;
begin
dblDecimal := decimal; //取得目标小数
if trunc(decimal) = decimal then // 如果是整数，则直转
    result := floattostr(decimal)
else
begin
    if abs(decimal) > 1 then //如果小数大于1 如 10.24 ，进行拆解
    begin
      dblinteger := trunc(decimal); //取出整数部分
      dblDecimal := abs(frac(decimal)); //取出小数部分
    end
    else dblDecimal := decimal;

    dblAccuracy := 0.01; //设置精度
    intNumerator := 0; //初始分子为0
    intDenominator := 1; //初始分母为1
    intNegative := 1; //符号标记为正
    if dblDecimal < 0 then intNegative := -1; //如果目标为负，设置负标志位
    dblFraction := 0; //设置分数值为 0/1
    while Abs(dblFraction - dblDecimal) > dblAccuracy do //如果当前没有达到精度要求就继续循环
    begin
      if Abs(dblFraction) > Abs(dblDecimal) then //如果我们的分数大于目标
        intDenominator := intDenominator + 1 //增加分母
      else //否则
        intNumerator := intNumerator + intNegative; //增加分子
      dblFraction := intNumerator / intDenominator; //计算新的分数
    end;
   // edit2.Text := inttostr(intNumerator) + '/' + inttostr(intDenominator);
    if abs(decimal) > 1 then //如果小数大于1 如 10.24 ，进行拆解
     //结果：dblinteger为整数部分， intNumerator为分子，intDenominator为分母
     begin
     if dblinteger<0 then
       begin dblinteger := -dblinteger;
       temp := floattostr(dblinteger*intDenominator+intNumerator) + '/' + inttostr(intDenominator);
       //result := floattostr(dblinteger) + '又' + inttostr(intNumerator) + '/' + inttostr(intDenominator)
       result := '-'+temp;
       end else
       begin
       temp := floattostr(dblinteger*intDenominator+intNumerator) + '/' + inttostr(intDenominator);
       result :=temp;
       end;
     end
    else
      result := inttostr(intNumerator) + '/' + inttostr(intDenominator);
end;
end;

procedure MatFromText(var A:RCOMat;path:string;var i,j:integer);
var f:TextFile;
    n,m:integer;
    Reg:TPerlRegEx;
    List:TStrings;
    s:string;
begin
i:=0; j:=0;
Reg:=TPerlRegEx.Create(nil);
List:=TStringList.Create;
n:=0;
assignFile(f,path);
reset(f);
while not eof(f) do
begin
    readln(f,s);
    if not (trim(s)='') then      //空行情况的处理
    i:=i+1;
end;
closefile(f);
reset(f);
s:='';
while not eof(f) do
begin
    Readln(f,s); s:=trim(s);
    if not (s='') then          //空行情况的处理
    begin
      reg.Subject := s;
      reg.RegEx := ' *';                  //以空格为分隔符（空格个数不限）
      reg.Split(list,MaxInt);
      j:=list.Count;
      SetLength(A,i,j);
      for m := 0 to j - 1 do
      begin
         A[n,m]:=strtofloat(list[m]);
      end;
       n:=n+1;
    end;
   list.Clear;
end;
list.Free;
closefile(f);
end;

procedure MatToMessageBox(A:RCOMat;i,j:integer;isT:boolean=false;isToText:boolean=true);
var m,n:integer;
    str:string;
    f:textfile;
begin
   str:='';
   for m := 0 to i - 1 do
   begin
     str:=str+#13;
     for n := 0 to j - 1 do
       begin
          if isT then
             str := str+Fraction(A[m,n])+'     '
          else
            str:=str+format('%.4f',[A[m,n]])+'     ';
       end;
   end;
   showmessage(str);
    str :='';

   if isToText then                              //将结果保存到文本文件
   begin
   assignfile(f,'.\data\result.txt');
   rewrite(f);                         //清空文本
   closefile(f);
   assignfile(f,'.\data\result.txt');
   Append(f);
   for m := 0 to i - 1 do
     for n := 0 to j - 1 do
     begin
      str:=str+format('%.4f',[A[m,n]])+'     ';
      if n=j-1 then
      begin
      writeln(f,str); str := '';
      end;
     end;
   closefile(f);
   end;
end;

FUNCTION MsgErreurMat(no:integer):string;
BEGIN
MsgErreurMat:='';
case no of
    MatCorrect   : Result:='succ鑣';
    MatNonInit   : Result:='unit?non initialis閑';
    MatDebordTas : Result:='pas assez de m閙oire';
    MatDetNul    : Result:='d閠erminant nul';
    end;
END;

FUNCTION NorMat(A:RCOMat;n,m:integer):extended;
VAR Nor : extended;
    i,j : integer;
BEGIN
Nor:=1;
for i:=0 to pred(n) do
    for j:=0 to pred(m) do
      Nor:=Nor*sqr(A[i,j]);
Nor:=sqrt(Nor/n/m);
Result:=Nor;
END;

FUNCTION DimMatrices:integer;
BEGIN
DimMatrices:=DimMat;
END;

PROCEDURE Init;
begin
if not Ini then
    begin
      ErreurMat:=MatDebordTas;
      write('Unit?MATRICES non initialis閑');// I think that write generate the write
      halt(1);                                // instruction would generate an exception
    end;
ErreurMat:=0;
end;

PROCEDURE InitMat(Dim:integer);
begin
DimMat:=Dim;
Ini:=true;
Init;
end;

PROCEDURE MatMat(A,B:RCOMat;var C:RCOMat);
var i, j , k : integer;
    Cij      : extended;

begin
Init;
if true
    then
      begin
        for i:=0 to pred(DimMat) do
          for j:=0 to pred(DimMat) do
            begin
              Cij:=0;
              for k:=0 to pred(DimMat) do
                Cij:=Cij+A[i,k]*B[k,j];
              C[i,j]:=Cij;
            end;
        ErreurMat:=0;
      end
    else ErreurMat:=MatDebordTas;
end;

PROCEDURE SomMat(A,B:RCOMat;var C:RCOMat);
var i,j : integer;

BEGIN
if true
    then
      begin
        for i:=0 to pred(DimMat) do
          for j:=0 to pred(DimMat) do
            C[i,j]:=A[i,j]+B[i,j];
        ErreurMat:=0;
      end
    else ErreurMat:=MatDebordTas;
END;

PROCEDURE MatMatR(A,B:RCOMat;M,P:integer;var C:RCOMat);
var i,j,k       : integer;
    Cij         : extended;
BEGIN
if true
    then
      begin
        for i:=0 to pred(DimMat) do
          for j:=0 to pred(P) do
            begin
              Cij:=0;
              for k:=0 to pred(M) do
                Cij:=Cij+A[i,k]*B[k,j];
              C[i,j]:=Cij;
            end;
        ErreurMat:=0;
      end
    else ErreurMat:=MatDebordTas;
END;

PROCEDURE TransMat(A:RCOMat;var B:RCOMat;n,m:integer); //n,m为A的行、列数
var i, j     : integer;
begin
Init;
if true
    then
      begin
        for i:=0 to pred(m) do
          for j:=0 to pred(n) do
            B[i,j]:=A[j,i];
        ErreurMat:=0;
      end
    else ErreurMat:=MatDebordTas;
end;

PROCEDURE MatUnt(var MU:RCOMat);
var i,j     : integer;
begin
Init;
for i:=0 to pred(DimMat) do
    for j:=0 to pred(DimMat) do
      if i=j then MU[i,j]:=1
             else MU[i,j]:=0;
end;

PROCEDURE MatVec(M:RCOMat;V:RCOVec;var W:RCOVec);
var i, j     : integer;
    Wi       : extended;
begin
Init;
if true
    then
      begin
        for i:=0 to pred(DimMat) do
          begin
            Wi:=0;
            for j:=0 to pred(DimMat) do
              Wi:=Wi+M[i,j]*V[j];
            W[i]:=Wi;
          end;
        ErreurMat:=0;
      end
    else ErreurMat:=MatDebordTas;
end;

FUNCTION PrdScal(U,V:RCOVec) : extended;
var i : integer;
    s : extended;
begin
Init;
s:=0;
for i:=0 to pred(DimMat) do
    s:=s+U[i]*V[i];
Result:=s;
end;

PROCEDURE ScalVec(x:extended;U:RCOVec;var V:RCOVec);
var i        : integer;
begin
Init;
for i:=0 to pred(DimMat) do
    V[i]:=x*U[i];
end;

PROCEDURE ScalMat(x:extended;A:RCOMat;var B:RCOMat);
var i,j : integer;
BEGIN
for i:=0 to pred(DimMat) do
    for j:=0 to pred(DimMat) do
      B[i,j]:=x*A[i,j];
END;

PROCEDURE AddVec (Signe:integer;U,V:RCOVec;var W:RCOVec);
var i        : integer;
begin
Init;
for i:=0 to pred(DimMat) do
    W[i]:=U[i]+Signe*V[i];
end;

PROCEDURE SomVec(U,V:RCOVec;var W:RCOVec);
begin
AddVec(1,U,V,W);
end;

PROCEDURE DifVec(U,V:RCOVec;var W:RCOVec);
begin
AddVec(-1,U,V,W);
end;

FUNCTION Norme(U:RCOVec)    : extended;
var i : integer;
    n : extended;
begin
Init;
n:=0;
for i:=0 to pred(DimMat) do
    n:=n+sqr(U[i]);
Result:=sqrt(n);
end;

PROCEDURE TrianguleO(A:RCOMat;var AT:RCOMat;var Jx,Iy:RCOVecE;var Det:extended);
var i, k     : integer;
    im , jm : integer;
    ATkk     : extended;
    C        : extended;
    Max      : extended;
    SignDet : integer;

PROCEDURE InitJxIy;
var i:integer;
begin
    for i:=0 to pred(DimMat) do
    begin
      Jx[i]:=i;
      Iy[i]:=i;
    end;
end;

PROCEDURE CopiAAT;
var i,j:integer;
begin
    for i:=0 to pred(DimMat) do
    for j:=0 to pred(DimMat) do
      AT[i,j]:=A[i,j];
end;

PROCEDURE Reduit;
var i,j:integer;
begin
    ATkk:=AT[k,k];
      for i:=k+1 to pred(DimMat) do
        begin
          C   :=AT[i,k]/ATkk;
          for j:=k+1 to pred(DimMat) do
            AT[i,j]:=AT[i,j]-C*AT[k,j];
          for j:=0 to k do
            if j<>k
              then AT[i,j]:=AT[i,j]-C*AT[k,j]
              else AT[i,j]:=-C;
        end;
end;

PROCEDURE ChercheMax(var Max:extended;var im,jm:integer);
var i,j : integer;
      M   : extended;
begin
    Max:=abs(AT[k,k]);
    im:=k;jm:=k;
    for i:=k to pred(DimMat) do
      for j:=k to pred(DimMat) do
        begin
          M:=abs(AT[i,j]);
          if M>Max then
            begin
              Max:=M;
              im:=i;jm:=j;
            end;
        end;
end;

PROCEDURE PermuteCol;
var Ind : integer;
      i   : integer;
      x   : extended;
begin
    for i:=0 to pred(DimMat) do
      begin
        x:=AT[i,k];
        AT[i,k]:=AT[i,jm];
        AT[i,jm]:=x;
      end;
    Ind   :=Jx[k];
    Jx[k] :=Jx[jm];
    Jx[jm]:=Ind;
    SignDet:=-SignDet;
end;

PROCEDURE PermuteLig;
var Ind : integer;
      j   : integer;
      x   : extended;
begin
    for j:=0 to pred(DimMat) do
      begin
        x:=AT[k,j];
        AT[k,j]:=AT[im,j];
        AT[im,j]:=x;
      end;
    Ind:=Iy[k];
    Iy[k]:=Iy[im];
    Iy[im]:=ind;
    SignDet:=-SignDet;
end;

begin
InitJxIy;
CopiAAT;
k:=0;
SignDet:=1;
while (k<=pred(DimMat)-1) and (ErreurMat=0) do
    begin
      ChercheMax(Max,im,jm);
      if Max<>0
        then
          begin
            if k<>jm then PermuteCol;
            if k<>im then PermuteLig;
            Reduit;
          end
        else ErreurMat:=MatDetNul;
      k:=k+1;
    end;
Det:=SignDet;
for i:=0 to pred(DimMat) do
    Det:=Det*AT[i,i];
if Det=0 then ErreurMat:=MatDetNul;
end;

PROCEDURE ResouTrO(A:RCOMat;var X,Y:RCOVec;Jx,Iy:RCOVecE);
var i, j     : integer;
    PVX      : RCOVec;
    Xi       : extended;

begin
if true
    then
      begin
        SetLength(PVX,DimMat);
        for i:=pred(DimMat) downto 0 do
          begin
            Xi:=Y[Iy[i]];
            for j:=0 to i-1 do
              Xi:=Xi+A[i,j]*Y[Iy[j]];
            for j:=i+1 to pred(DimMat) do
              Xi:=Xi-A[i,j]*PVX[j];
            PVX[i]:=Xi/A[i,i];
          end;
        for i:=0 to pred(DimMat) do
          X[Jx[i]]:=PVX[i];
        Setlength(PVX,0);
        ErreurMat:=0;
      end
    else ErreurMat:=MatDebordTas;
end;

PROCEDURE SysLin(A:RCOMat;var X:RCOVec;B:RCOVec;var Det:extended);
var AT       : RCOMat;
    Jx,Iy    : RCOVecE;

PROCEDURE Alloc;
begin
    SetLength(AT,DimMat,DimMat);
    SetLength(Jx,DimMat);
    SetLength(Iy,DimMat);
end;

PROCEDURE DesAlloc;
begin
    SetLength(Iy,0);
    SetLength(Jx,0);
    SetLength(At,0,0);
end;

begin
Init;
if true
    then
      begin
        Alloc;
        TrianguleO(A,AT,Jx,Iy,Det);
        if Det<>0
          then ResouTrO(AT,X,B,Jx,Iy);
        DesAlloc;
      end
    else ErreurMat:=MatDebordTas;
end;

PROCEDURE InvMat(A:RCOMat;var B:RCOMat;var Det:extended);
var j        : integer;
    AT       : RCOMAt;
    PVX      : RCOVec;
    Jx,Iy    : RCOVecE;

PROCEDURE Alloc;
begin
    Setlength(AT,DimMat,DimMat);
    SetLength(PVX,DimMat);
    SetLength(Jx,DimMat);
    SetLength(Iy,DimMat);
end;

PROCEDURE DesAlloc;
begin
    SetLength(Iy,0);
    SetLength(Jx,0);
    SetLength(PVX,0);
    SetLength(AT,0,0);
end;

PROCEDURE initVX;
var i : integer;
begin
    for i:=0 to pred(DimMat) do
      PVX[i]:=0;
    PVX[j]:=1;
end;

PROCEDURE TransfertCol;
var i : integer;
begin
    for i:=0 to pred(DimMat) do
      B[i,j]:=PVX[i];
end;

begin
Init;
if true
    then
      begin
        Alloc;
        TrianguleO(A,AT,Jx,Iy,Det);
        if ErreurMat=0
          then
            begin
              for j:=0 to pred(DimMat) do
              begin
                InitVX;
                ResouTrO(AT,PVX,PVX,Jx,Iy);
                TransfertCol;
              end;
            end;
        DesAlloc;
      end
    else ErreurMat:=MatDebordTas;
end;

initialization
Ini:=false;
DimMat:=0;
END.

// Test program :
{
procedure TForm1.FormCreate(Sender: TObject);
Const Dim = 200;
var   A,B,C : RCOMat;
      V,W,U : RCOVec;
      i,j   : integer;
      Det   : extended;
begin
InitMat(Dim);
SetLength(A,Dim,Dim);
SetLength(B,Dim,Dim);
SetLength(C,Dim,Dim);
Randomize;
for i:=0 to pred(Dim) do
    begin
      V[i]:=random*10;
      for j:=0 to pred(Dim) do
      A[i,j]:=random*10;
    end;
InvMat(A,B,Det);
MatMat(A,B,C);
// [C] should be the unit matrix
end;
}