summaryrefslogtreecommitdiff
path: root/vecmath.pas
diff options
context:
space:
mode:
Diffstat (limited to 'vecmath.pas')
-rw-r--r--vecmath.pas538
1 files changed, 538 insertions, 0 deletions
diff --git a/vecmath.pas b/vecmath.pas
new file mode 100644
index 0000000..5e374dd
--- /dev/null
+++ b/vecmath.pas
@@ -0,0 +1,538 @@
+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.
generated by cgit v1.2.3-54-g00ecf (git 2.44.0) at 2024-05-16 12:17:06 +0000