unit vecMath;

interface

uses math, Types, SysUtils, ComCtrls;

Type TDetMethode = (dmLaPlace,dmLaPlace0,dmGauss);
     TMatrix = class
               private
                _zanz,_sanz: integer;
                Komps: array of array of extended; // [zeile,spalte]
                procedure wsanz(sanz: integer);
                procedure wzanz(zanz: integer);
                function ra(z,s: integer): extended;
                procedure wa(z,s: integer; _a: extended);
                function rPart(l,t,r,b: integer): TMatrix;
               public
                Progressbar: TProgressbar;
                property xdim: integer read _sanz
                                       write wsanz;
                property ydim: integer read _zanz
                                       write wzanz;
                property a[z,s: integer]: extended read ra
                                                   write wa;
                property part[l,t,r,b: integer]: TMatrix read rPart;
                constructor create;
                destructor destroy; override;
                function asPoint: TPoint; overload;
                function asPoint(scale: extended): TPoint; overload;
                procedure assign(M: TMatrix);
                procedure assignId(dim: integer);
                procedure assign0;
                function duplicate: TMatrix;
                procedure plus(M: TMatrix);
                procedure minus(M: TMatrix);
                procedure times(x: extended);
                procedure times_(R: TMatrix);
                procedure _times(L: TMatrix);
                procedure cross(M: TMatrix);
                function det: extended; overload;
                function det(methode: TDetMethode): extended; overload;
                function det_laPlace(spare_optimize: boolean): extended;
                function det_gauss: extended;
                procedure transp;
                procedure orthogonalize; overload;
                procedure orthogonalize(normalize: boolean); overload;
                function len: extended; // only for vectors!
                function zlen(z: integer): extended;
                function slen(s: integer): extended;
                function toStr: string;
                function fromString(str: string): boolean;
               end;

implementation

function istZahl(str: string): boolean;
var hK: boolean;
begin
 result:=false;
 if length(str)=0 then exit;
 hK:=false;
 if str[1]='-' then delete(str,1,1);
 if length(str)=0 then exit;
 if not (str[1] in ['0'..'9']) then exit;
 while (length(str)>0) and ((str[1] in ['0'..'9']) or
                            ((str[1] = ',') and not hK)) do
  begin
   hK:=hK or (str[1]=',');
   delete(str,1,1);
  end;
 result:=length(str)=0;
end;

procedure TMatrix.wsanz(sanz: integer);
var i: integer;
begin
 if sanz>=0 then _sanz:=sanz;
 setlength(Komps,_zanz);
 for i:=0 to _zanz-1 do
  setlength(Komps[i],_sanz);
end;

procedure TMatrix.wzanz(zanz: integer);
var i: integer;
begin
 if zanz>=0 then _zanz:=zanz;
 setlength(Komps,_zanz);
 for i:=0 to _zanz-1 do
  setlength(Komps[i],_sanz);
end;

function TMatrix.ra(z,s: integer): extended;
begin
 if (z>=0) and (z<_zanz) and
    (s>=0) and (s<_sanz) then result:=Komps[z,s]
                         else result:=0;
end;

procedure TMatrix.wa(z,s: integer; _a: extended);
begin
 if (z>=0) and (z<_zanz) and
    (s>=0) and (s<_sanz) then Komps[z,s]:=_a;
end;

function TMatrix.rPart(l,t,r,b: integer): TMatrix;
var z,s: integer;
begin
 l:=max(l,0);
 t:=max(t,0);
 r:=max(l,r);
 b:=max(t,b);
 r:=min(r,xdim-1);
 b:=min(b,ydim-1);
 result:=TMatrix.create;
 result.ydim:=b-t+1;
 result.xdim:=r-l+1;
 for z:=t to b do
  for s:=l to r do
   result.a[z-t,s-l]:=a[z,s];
end;

constructor TMatrix.create;
begin
 inherited create;
 _zanz:=3;
 xdim:=1;
 Progressbar:=nil;
end;

destructor TMatrix.destroy;
begin
 inherited destroy;
end;

function TMatrix.asPoint: TPoint;
begin
 result:=asPoint(1);
end;

function TMatrix.asPoint(scale: extended): TPoint;
begin
 if xdim>ydim then result:=Point(round(a[0,0]*scale),
                                 round(a[0,1]*scale))
              else result:=Point(round(a[0,0]*scale),
                                 round(a[1,0]*scale));
end;

procedure TMatrix.assign(M: TMatrix);
var z,s: Integer;
begin
 ydim:=M.ydim;
 xdim:=M.xdim;
 for z:=0 to ydim-1 do
  for s:=0 to xdim-1 do
   a[z,s]:=M.a[z,s];
end;

function TMatrix.duplicate: TMatrix;
begin
 Result:=TMatrix.create;
 Result.assign(self);
end;

procedure TMatrix.plus(M: TMatrix);
var z,s: Integer;
begin
 if (M.xdim = xdim) and
    (M.ydim = ydim) then
  for z:=0 to ydim-1 do
   for s:=0 to xdim-1 do
    a[z,s]:=a[z,s]+M.a[z,s];
end;

procedure TMatrix.minus(M: TMatrix);
var z,s: Integer;
begin
 if (M.xdim = xdim) and
    (M.ydim = ydim) then
  for z:=0 to ydim-1 do
   for s:=0 to xdim-1 do
    a[z,s]:=a[z,s]-M.a[z,s];
end;

procedure TMatrix.times(x: extended);
var z,s: Integer;
begin
 for z:=0 to ydim-1 do
  for s:=0 to xdim-1 do
   a[z,s]:=x*a[z,s];
end;

function TMatrix.det: extended;
begin
  result:=det_laPlace(true);
end;

function TMatrix.det(methode: TDetMethode): extended;
begin
  case methode of
    dmLaPlace:  result:=det_laPlace(false);
    dmLaPlace0: result:=det_laPlace(true);
    dmGauss:    result:=det_gauss;
    else result:=-1;
  end{of case};
end;

function TMatrix.det_laPlace(spare_optimize: boolean): extended;
var p:      array of Integer;
    i,j,k,
    tiefe:  Integer;
    tmp:    extended;
    UT:     boolean;
begin
 result:=0;
 if xdim<>ydim then exit;
 setlength(p,xdim);
 for i:=0 to xdim-1 do
  p[i]:=i;
 if assigned(Progressbar) then begin
   Progressbar.Position:=0;
   Progressbar.Min:=0;
   Progressbar.Max:=1;
   tiefe:=-1;
   repeat
     inc(tiefe);
     i:=0;
     for j:=0 to xdim-1 do
       if a[tiefe,j]<>0 then inc(i);
     Progressbar.Max:=
       Progressbar.Max*i;
   until (Progressbar.Max>=100) or (tiefe>=xdim-1);
   Progressbar.Step:=1;
   Progressbar.Visible:=true;
 end
 else tiefe:=-2;
 repeat
  tmp:=1;
  for i:=0 to xdim-1 do
   begin
    for j:=i+1 to xdim-1 do
     if p[i]>p[j] then tmp:=-tmp;
    tmp:=tmp*a[i,p[i]];
   end;
  result:=result+tmp;
  if spare_optimize then begin
    i:=0;
    while (i<xdim-1) and (a[i,p[i]]<>0) do
      inc(i);
    for j:=i+1 to xdim-1 do
      p[j]:=-1;
  end
  else
    i:=xdim-1;
  UT:=true;
  while UT and (i>=0) do
   begin
    repeat
     inc(p[i]);
     UT:=spare_optimize and (p[i]<xdim) and (a[i,p[i]]=0);
     for j:=0 to i-1 do
      UT:=UT or (p[i]=p[j]);
    until not UT;
    if p[i]>=xdim then
     begin
      for j:=i to xdim-1 do
       p[i]:=-1;
      dec(i);
      UT:=true;
     end;
   end;
  if i<=tiefe then
    Progressbar.StepIt;
  if i>=0 then
   for j:=i+1 to xdim-1 do
    repeat
     inc(p[j]);
     UT:=false;
     for k:=0 to j-1 do
      UT:=UT or (p[j]=p[k]);
    until not UT;
 until UT and (i<0);
 if assigned(Progressbar) then
   Progressbar.Visible:=false;
end;

function TMatrix.det_gauss: extended;
var M:       TMatrix;
    I,J,K:   Longint;
    tmp,fak: extended;
    vz:      Longint;
begin
  result:=0;
  if xdim<>ydim then exit;
  M:=TMatrix.create;
  M.assign(self);
  vz:=1;
  for I:=xdim-1 downto 0 do begin
    J:=I;
    while (J>=0) and (M.a[I,J]=0) do
      dec(J);
    if J<0 then begin
      result:=0;
      M.free;
      exit;
    end;
    if J<>I then
      for K:=0 to xdim-1 do begin
        tmp:=M.a[K,I];
        M.a[K,I]:=M.a[K,J];
        M.a[K,J]:=tmp;
        vz:=-vz;
      end;
    for J:=I-1 downto 0 do begin
      fak:=M.a[J,I]/M.a[I,I];
      for K:=0 to xdim-1 do
        M.a[J,K]:=M.a[J,K] - fak * M.a[I,K];
    end;
  end;
  result:=vz;
  for I:=0 to xdim-1 do
    result:=result*M.a[I,I];
  M.free;
end;

procedure TMatrix.transp;
var xd,yd,z,s: integer;
    tmp:       extended;
begin
 xd:=ydim;
 yd:=xdim;
 xdim:=max(xd,yd);
 ydim:=max(xd,yd);
 for z:=0 to ydim-1 do
  for s:=z+1 to xdim-1 do
   begin
    tmp:=a[z,s];
    a[z,s]:=a[s,z];
    a[s,z]:=tmp;
   end;
 xdim:=xd;
 ydim:=yd;
end;

function TMatrix.toStr: string;
var z,s: integer;
begin
 result:='';
 for z:=0 to ydim-1 do
  begin
   result:=result+#$0D#$0A;
   for s:=0 to xdim-1 do
    result:=result+floattostr(a[z,s])+' ';
   delete(result,length(result),1);
  end;
 delete(result,1,2);
end;

function TMatrix.fromString(str: string): boolean;
var str2:  string;
    z,s,i: integer;
begin
 result:=false;
 if length(str)=0 then exit;
 i:=length(str);
 if str[length(str)]<>#13 then str:=str+#13;
 while pos(#$0A,str)>0 do
  delete(str,pos(#$0A,str),1);
 while (i>0) do
  begin
   if str[i]=#13 then str:=copy(str,1,i-1)+' '+copy(str,i,length(str)+1-i);
   dec(i);
  end;
 str2:=str;
 for z:=0 to ydim-1 do
  begin
   for s:=0 to xdim-1 do
    begin
     while (length(str)>0) and (str[1]=' ') do
      delete(str,1,1);
     if pos(' ',str)=0 then exit;
     if not istZahl(copy(str,1,pos(' ',str)-1)) then exit;
     delete(str,1,pos(' ',str));
     while (length(str)>0) and (str[1]=' ') do
      delete(str,1,1);
    end;
   if (length(str)=0) or not (str[1]=#$0D) then exit;
   delete(str,1,1);
  end;
 result:=true;
 for z:=0 to ydim-1 do
  begin
   for s:=0 to xdim-1 do
    begin
     while (length(str2)>0) and (str2[1]=' ') do
      delete(str2,1,1);
     if pos(' ',str2)=0 then exit;
     a[z,s]:=strtofloat(copy(str2,1,pos(' ',str2)-1));
     delete(str2,1,pos(' ',str2));
     while (length(str2)>0) and (str2[1]=' ') do
      delete(str2,1,1);
    end;
   if (length(str2)=0) or not (str2[1]=#$0D) then exit;
   delete(str2,1,1);
  end;
end;

procedure TMatrix.times_(R: TMatrix);
var L:     TMatrix;
    z,s,i: integer;
begin
 if R.ydim<>xdim then exit;
 L:=duplicate;
 xdim:=R.xdim;
 for z:=0 to ydim-1 do
  for s:=0 to xdim-1 do
   begin
    a[z,s]:=0;
    for i:=0 to R.ydim-1 do
     a[z,s]:=a[z,s] + L.a[z,i]*R.a[i,s];
   end;
end;

procedure TMatrix._times(L: TMatrix);
var R:     TMatrix;
    z,s,i: integer;
begin
 if ydim<>L.xdim then exit;
 R:=duplicate;
 ydim:=L.ydim;
 for z:=0 to ydim-1 do
  for s:=0 to xdim-1 do
   begin
    a[z,s]:=0;
    for i:=0 to R.ydim-1 do
     a[z,s]:=a[z,s] + L.a[z,i]*R.a[i,s];
   end;
end;

procedure TMatrix.cross(M: TMatrix);
var x,y,z: extended;
begin
 if (xdim<>1) or (M.xdim<>1) or (ydim<>3) or (M.ydim<>3) then exit; // ziemlich speziell ...
 x:=a[1,0]*M.a[2,0]-a[2,0]*M.a[1,0];
 y:=a[2,0]*M.a[0,0]-a[0,0]*M.a[2,0];
 z:=a[0,0]*M.a[1,0]-a[1,0]*M.a[0,0];
 a[0,0]:=x;
 a[1,0]:=y;
 a[2,0]:=z;
end;

procedure TMatrix.orthogonalize;
begin
 orthogonalize(false);
end;

procedure TMatrix.orthogonalize(normalize: boolean);
var zT,z,s:      integer;
    diff,me,opp: TMatrix;
    l:           extended;
begin
 if det=0 then exit;

 diff:=TMatrix.create;
 diff.xdim:=xdim;
 diff.ydim:=1;
 for zT:=0 to ydim-1 do
  begin
   diff.assign0;
   for z:=0 to zT-1 do
    begin
     me:=part[0,zT,xdim-1,zT];
     l:=me.len;
     me.transp;
     opp:=part[0,z,xdim-1,z];
     me._times(opp);
     l:=me.a[0,0]/l;
     me.destroy;
     opp.times(l);
     diff.plus(opp);
     opp.destroy;
    end;
   for s:=0 to xdim-1 do
    a[zT,s]:=a[zT,s]-diff.a[0,s];
   if normalize then
    begin
     me:=part[0,zT,xdim-1,zT];
     l:=me.len;
     me.destroy;
     for s:=0 to xdim-1 do
      a[zT,s]:=a[zT,s]/l;
    end;
  end;
end;

procedure TMatrix.assignId(dim: integer);
var z,s: integer;
begin
 if dim<1 then exit;
 xdim:=dim;
 ydim:=dim;
 for z:=0 to dim-1 do
  for s:=0 to dim-1 do
   a[z,s]:=byte(z=s);
end;

procedure TMatrix.assign0;
var z,s: integer;
begin
 for z:=0 to ydim-1 do
  for s:=0 to xdim-1 do
   a[z,s]:=0;
end;

function TMatrix.len: extended;
begin
 if xdim=1 then result:=slen(0)
           else result:=zlen(0);
end;

function TMatrix.zlen(z: integer): extended;
var s: integer;
begin
 result:=0;
 for s:=0 to xdim-1 do
  result:=result + sqr(a[z,s]);
 result:=sqrt(result);
end;

function TMatrix.slen(s: integer): extended;
var z: integer;
begin
 result:=0;
 for z:=0 to ydim-1 do
  result:=result + sqr(a[z,s]);
 result:=sqrt(result);
end;

end.