eigen.pas

Eigen mathematical function for delphi 5 and higher

function sgn(x:double):double;
begin
 if x<0 then sgn:=-1 else sgn:=1;
 if x=0 then sgn:=0;
end;

function Eigen(Vstup:Variant;Values:Boolean):Variant;
var M,i,j,k,l,IU,IW:Integer;
    E,U,V,r36,r37,r38,r39,r40,Pom:Double;
    Z,S,R:Variant;
begin
 if not(VarIsArray(Vstup)) then begin
  if Values then begin
   Eigen:=VarArrayCreate([1,1],varDouble);
   Eigen[1]:=Vstup;
  end
  else begin
   Eigen:=VarArrayCreate([1,1,1,1],varDouble);
   Eigen[1,1]:=1;
  end;
  Exit;
 end;
 M:=VarArrayHighBound(Vstup,1);
 if M<>VarArrayHighBound(Vstup,2) then begin
  Eigen:='N/A';
  Exit;
 end;
 if M<1 then begin
  Eigen:='N/A';
  Exit;
 end;
 if M=1 then begin
  if Values then begin
   Eigen:=VarArrayCreate([1,1],varDouble);
   Eigen[1]:=Vstup[1,1];
  end
  else begin
   Eigen:=VarArrayCreate([1,1,1,1],varDouble);
   Eigen[1,1]:=1;
  end;
  Exit;
 end;
 Z:=VarArrayCreate([1,M,1,M],varDouble);
 S:=VarArrayCreate([1,M],varDouble);
 R:=VarArrayCreate([1,M],varDouble);
 for i:=1 to M-1 do
  for j:=i+1 to M do
   Vstup[i,j]:=Vstup[j,i];
 for i:=1 to M do begin
  S[i]:=0;
  R[i]:=0;
  for j:=1 to M do
   Z[i,j]:=Vstup[i,j];
 end;
 for i:=M downto 2 do begin
  IU:=i-1;
  V:=0;
  E:=0;
  if IU>=2 then begin
   for k:=1 to IU do
    V:=V+Abs(Z[i,k]);
    if V<>0 then begin
     for k:=1 To IU do begin
      Z[i,k]:=Z[i,k]/V;
      E:=E+Z[i,k]*Z[i,k];
     end;
     r36:=Z[i,IU];
     Pom:=r36+1E-36;
     r37:=-Sqrt(Abs(E))*Sgn(Pom);
     S[i]:=V*r37;
     E:=E-r36*r37;
     Z[i,IU]:=r36-r37;
     r36:=0;
     for j:=1 to IU do begin
      Z[j,i]:=Z[i,j]/(V*E);
      r37:=0;
      for k:=1 to j do r37:=r37+Z[j,k]*Z[i,k];
      if IU>=(j+1) then
       for k:=j+1 to IU do r37:=r37+Z[k,j]*Z[i,k];
      S[j]:=r37/E;
      r36:=r36+S[j]*Z[i,j];
     end;
     r38:=r36/(E+E);
     for j:=1 to IU do begin
      r36:=Z[i,j];
      r37:=S[j]-r38*r36;
      S[j]:=r37;
      for k:=1 to j do Z[j,k]:=Z[j,k]-r36*S[k]-r37*Z[i,k];
     end;
     for k:=1 to IU do Z[i,k]:=V*Z[i,k];
    end
    else begin
     //else v<>0
     S[i]:=Z[i,IU];
    end;
  end
  else begin
   //else IU>=2
   S[i]:=Z[i,IU];
  end;
  R[i]:=E;
 end;
 S[1]:=0;
 R[1]:=0;
 for i:=1 To M do begin
  IU:=i-1;
  if R[i]<>0 then
   for j:=1 to IU do begin
    r37:=0;
    for k:=1 to IU do r37:=Z[i,k]*Z[k,j]+r37;
    for k:=1 to IU do Z[k,j]:=Z[k,j]-r37*Z[k,i];
   end;
  R[i]:=Z[i,i];
  Z[i,i]:=1;
  if IU>=1 then
   for j:=1 to IU do begin
    Z[i,j]:=0;
    Z[j,i]:=0;
   end;
 end;
 for i:=2 To M do S[i-1]:=S[i];
 S[M]:=0;
 U:=0;
 V:=0;
 for l:=1 to M do begin
  J:=0;
  E:=(0.000000001)*(Abs(R[l])+Abs(S[l]));
  if V<E then V:= E;
  IW:=l;
  while ((Abs(S[IW])>V) and (IW<M)) do Inc(IW);
  if IW<>l then
  Repeat
   if j=100 then begin
    Eigen:='N/A';
    Exit;
   end;
   Inc(j);
   r36:=(R[l+1]-R[l])/(2*S[l]);
   r37:=Sqrt(Abs(r36*r36+1));
   Pom:=r36+1E-35;
   E:=R[l]-S[l]/(r36+r37*Sgn(Pom));
   for i:=l to M do R[i]:=R[i]-E;
   U:=U+E;
   r36:=R[IW];
   r38:=1;
   r39:=0;
   for i:=(IW-1) downto l do begin
    r40:=r38*S[i];
    E:=r38*r36;
    if Abs(r36)>=Abs(S[i]) then begin
     r38:=S[i]/r36;
     r37:=Sqrt(Abs(r38*r38+1));
     S[i+1]:=r39*r36*r37;
     r39:=r38/r37;
     r38:=1/r37;
    end
    else begin
     r38:=r36/S[i];
     r37:=Sqrt(Abs(r38*r38+1));
     S[i+1]:=r39*S[i]*r37;
     r39:=1/r37;
     r38:=r38*r39;
    end;
    r36:=r38*R[i]-r39*r40;
    R[i+1]:=E+r39*(r38*r40+r39*R[i]);
    for k:=1 to M do begin
     E:=Z[k,i+1];
     Z[k,i+1]:=r39*Z[k,i]+r38*E;
     Z[k,i]:=r38*Z[k,i]-r39*E;
    end;
   end;
   S[l]:=r39*r36;
   R[L]:=r38*r36;
   if Abs(S[l])=V then Break;
  Until Abs(S[l])<=V;
  R[l]:=U+R[l];
 end;
 For i:=1 to M-1 do begin
  k:=i;
  r36:=R[i];
  for j:=1+1 to M do
   if R[j]<r36 then begin
    k:=j;
    r36:=R[j];
   end;
  if i<>k then begin
   R[k]:=R[i];
   R[i]:=r36;
   for j:=1 to M do begin
    r36:=Z[j,i];
    Z[j,i]:=Z[j,k];
    Z[j,k]:=r36;
   end;
  end;
 end;
 if Values then
  Eigen:=R
 else
  Eigen:=Z;
end;