home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / pascal / library / dos / math / matrix / matrixu.pas < prev   
Encoding:
Pascal/Delphi Source File  |  1988-06-13  |  11.9 KB  |  355 lines
  1. unit MATRIXU;
  2.  
  3.      {by       Josh Mitteldorf
  4.                6800 Scotforth Rd
  5.                Philadelphia, PA 19119
  6.                70441,1305 }
  7.  
  8. INTERFACE
  9. uses CRT;
  10.  
  11. type float     = double;
  12.     floatarray = array[0..99] of float;
  13.      vector    = ^floatarray;
  14.      matrix    = ^floatarray;
  15.  
  16. {You can change the definition of float above, to real or single or extended,
  17.  and all computations will be to the specified precision.  You must change to
  18.  float = real in order to run this unit without a coprocessor.
  19.  
  20.  In all these procedures and functions, the last argument d is the dimension.
  21.  All matrices are square, of dimension d x d.   Everything is accessed as a
  22.  pointer in order to facilitate the flexibility of being able to specify d when
  23.  the procedure is called.  Thus matrix and vector are just two different names
  24.  for the same thing, namely a pointer to array of type float.
  25.  
  26.  The calling unit may be structured as follows:
  27.  
  28.    const n=4;
  29.  
  30.    type float=real;
  31.         vector=array[1..n] of float;
  32.         matrix=array[1..n] of vector;
  33.  
  34.    var x,y    :vector;
  35.        z      :matrix;
  36.  
  37.     ...
  38.  
  39.     Prod(addr(z),addr(x),addr(y));
  40.     ...
  41.  
  42.    Then the element z[i,j] in the calling unit becomes a^[i*(d-1)+j] in the
  43.    present unit.
  44.  
  45.    CAUTION:  You can't declare the size of a matrix larger than what you
  46.              actually use as dimension in calling matrix procedures from this
  47.              unit.  In other words, the constant n in the above example
  48.              MUST be used as d in all calls to matrix functions in this unit.
  49.              (This is so that the elements of the matrix are all contiguous
  50.              in memory.)  For calls to functions that involve vectors only,
  51.              you may specify any d you want.
  52.  
  53. }
  54.  
  55. procedure ScalarMult(a :matrix; s :float;  b: matrix; d:integer);
  56.     {Multiplies matrix a by scalar s to give matrix b}
  57. procedure ScalarMultV( a :vector; s :float;  b :vector; d:integer);
  58.     {Multiplies vector a by scalar s to give vector b}
  59. procedure VecSum( a,b,c :vector;  d:integer);
  60.     {Adds vector a to vector b to give vector c}
  61. procedure VecDif( a,b,c :vector;  d:integer);
  62.     {Subtracts vector b from vector a to give vector c}
  63. procedure Diadic( a,b :vector;  c:matrix; d:integer);
  64.     {Multiplies every element of vector a by every element of b
  65.      to give the symmetric matrix c.}
  66. procedure VecPrint ( v:vector; d: integer);
  67.     {Prints a vector to the screen in one.}
  68. procedure MatPrint ( m: matrix; d:integer);
  69.     {Prints a matrix to the screen in one.}
  70. procedure Transpose( a :matrix;  b :matrix; d:integer);
  71.     {Transposes b[i,j]:=a[j,i].}
  72. procedure MatSum( a,b :matrix;  c :matrix; d:integer);
  73.     {Adds matrix a to matrix b to give matrix c.}
  74. procedure MatDif( a,b :matrix;  c :matrix; d:integer);
  75.     {Subtracts matrix b from matrix a to give matrix c.}
  76. procedure Cross( a,b,c :vector);
  77.     {Vector c is the cross product a x c. (Only makes sense for d=3.) }
  78. procedure Prod( a:matrix;  b,c :vector; d :integer);
  79.     {Matrix a is multiplied by vector b to give vector c.}
  80. procedure MatProd( a,b,c :matrix;  d:integer);
  81.     {Matrix a is multiplied by matrix b to give matrix c.}
  82. procedure Inverse( a,b :matrix;  d :integer);
  83.     {Matrix b is the inverse of matrix a.}
  84. procedure Solve( a :matrix;  b,c :vector;  d :integer);
  85.     {solves a.c=b for c, where a is a matrix and b and c are vectors.}
  86. procedure MatQuot( a,b,c :matrix;  d:integer);
  87.     {Divides matrix a by matrix b to give matrix c.}
  88. function Dot( a,b :vector; d:integer):float;
  89.     {Returns the scalar product of vectors a and b.}
  90. function Mag( a: vector; d:integer):float;
  91.     {Returns the sqrt of the sum of the squares of elements of a.
  92.      To obtain the magnitude of a square matrix, call Mag with last
  93.      argument = d².}
  94. function Det(aa:matrix; d :integer):float;
  95.     {Returns the determinant of matrix aa.}
  96.  
  97. IMPLEMENTATION
  98.  
  99. procedure ScalarMult;
  100.           var i,j   :integer;
  101.           begin
  102.           for i:=0 to pred(d) do
  103.               for j:=0 to pred(d) do
  104.                   b^[d*i+j]:=a^[d*i+j]*s;
  105.           end;
  106.  
  107. procedure ScalarMultV;
  108.           var i   :integer;
  109.           begin
  110.           for i:=0 to pred(d) do b^[i]:=a^[i]*s;
  111.           end;
  112.  
  113. procedure VecSum;
  114.           var i :integer;
  115.           begin
  116.           for i:=0 to pred(d) do c^[i]:=a^[i]+b^[i];
  117.           end;
  118.  
  119. procedure VecDif;
  120.           var i   :integer;
  121.           begin
  122.           for i:=0 to pred(d) do c^[i]:=a^[i]-b^[i];
  123.           end;
  124.  
  125. procedure Diadic;
  126.           var i,j   :integer;
  127.           begin
  128.           for i:=0 to pred(d) do
  129.               for j:=0 to pred(d) do
  130.                   c^[d*i+j]:=a^[i]*b^[j];
  131.           end;
  132.  
  133. procedure VecPrint;
  134.           var i: integer;
  135.           begin
  136.           for i:=0 to pred(d) do write(v^[i]:7:4,' ');
  137.           writeln;
  138.           end;
  139.  
  140. procedure MatPrint;
  141.           var i,j: integer;
  142.           begin
  143.           for i:=0 to pred(d) do begin
  144.               for j:=0 to pred(d) do write(m^[d*i+j]:7:4,' ');
  145.               writeln;
  146.               end;
  147.           end;
  148.  
  149. procedure Transpose;
  150.           var i,j: integer;
  151.           begin
  152.           for i:=0 to pred(d) do
  153.               for j:=0 to pred(d) do b^[d*i+j]:=a^[j*d+i];
  154.           end;
  155.  
  156. function Dot;
  157.          var i :integer; partial :float;
  158.          begin
  159.          partial:=0;
  160.          for i:=0 to pred(d) do partial:=partial+a^[i]*b^[i];
  161.          dot:= partial;
  162.          end;
  163.  
  164. function Mag;
  165.          begin
  166.          Mag:=sqrt(Dot(a,a,d));
  167.          end;
  168.  
  169. procedure MatSum;
  170.           var i,j  :integer;
  171.           begin
  172.           for i:=1 to d do
  173.               for j:=1 to d do
  174.                   c^[d*i+j]:=a^[d*i+j]+b^[d*i+j];
  175.           end;
  176.  
  177. procedure MatDif;
  178.           var i,j  :integer;
  179.           begin
  180.           for i:=0 to pred(d) do
  181.               for j:=0 to pred(d) do
  182.                   c^[d*i+j]:=a^[d*i+j]-b^[d*i+j];
  183.           end;
  184.  
  185. procedure Cross;
  186.          begin
  187.          c^[1]:=a^[2]*b^[3] - a^[3]*b^[2];
  188.          c^[2]:=a^[3]*b^[1] - a^[1]*b^[3];
  189.          c^[3]:=a^[1]*b^[2] - a^[2]*b^[1];
  190.          end;
  191.  
  192. procedure Prod;
  193.           var i,j,d1 :integer;  { Multiplies square matrix a by vector b }
  194.              sum     :float;    { to give vector c }
  195.              size    :word;
  196.  
  197.           begin
  198.           d1:=pred(d);
  199.           size:=d*sizeof(float);
  200.           for i:=0 to d1 do begin
  201.               sum:=0;
  202.               for j:=0 to d1 do sum:=sum+b^[j]*a^[d*i+j];
  203.               c^[i]:=sum;
  204.               end;
  205.           end;
  206.  
  207. function Det;
  208.          var i,j,k            :integer;
  209.              wr               :float;
  210.              a                :matrix;
  211.              size             :word;
  212.          begin
  213.          if (d>3) then begin
  214.             size:=sqr(d)*sizeof(float);
  215.             GetMem(a,size);
  216.             move(aa^,a^,size);
  217.             for i:=0 to d-2 do begin
  218.                 for j:= succ(i) to pred(d) do begin
  219.                     wr:=a^[i*d+i];
  220.                     if (wr=0) then begin det:=0; exit end;
  221.                     wr:=a^[j*d+i]/wr;
  222.                     for k:= succ(i) to pred(d) do a^[j*d+k]:= a^[j*d+k] - a^[i*d+k]*wr;
  223.                     end;
  224.                 end;
  225.             wr:=a^[0];
  226.             for i:=1 to pred(d) do wr:=wr*a^[i*d+i];
  227.             Det:=wr;
  228.             FreeMem(a,size);
  229.             end
  230.          else if (d=3) then Det:= aa^[0]*aa^[4]*aa^[8] +aa^[1]*aa^[5]*aa^[6] +aa^[2]*aa^[3]*aa^[7]
  231.                                  -aa^[0]*aa^[7]*aa^[5] -aa^[1]*aa^[3]*aa^[8] -aa^[2]*aa^[4]*aa^[6]
  232.          else if (d=2) then Det:=aa^[0]*aa^[3] - aa^[1]*aa^[2]
  233.          else if (d=1) then Det:=aa^[0]
  234.          else writeln('ERROR: DETERMINANT OF MATRIX OF ORDER ',d);
  235.          end;
  236.  
  237. {
  238.   THIS FUNCTION works fine but takes time proportional to factorial d,
  239.    so it's a useless mathematical curiosity.
  240.    The function above that replaces it takes time proportional to d^3.
  241. function Det;
  242.          var minor :matrix;
  243.              partial :float;
  244.              i,j,dmin,plusminus :integer;
  245.              size               :word;
  246.  
  247.          begin
  248.          size:=sqr(pred(d))*sizeof(float));
  249.          GetMem(minor,size);
  250.          dmin:=pred(d);
  251.          if (d>3) then begin
  252.             partial:=0;
  253.             plusminus:=1;
  254.             minor:=a^;
  255.             for i:=d downto 1 do begin
  256.                 partial:=partial+plusminus*a^[i*d+d]*det(minor,dmin);
  257.                 plusminus:=-plusminus;
  258.                 if i>1 then minor[pred(i)]:=a^[i];
  259.                 end;
  260.             Det:=partial;
  261.             end
  262.          else if (d=3) then
  263.          Det:= a^[1*d+1]*a^[2*d+2]*a^[3*d+3] +a^[1*d+2]*a^[2*d+3]*a^[3*d+1] +a^[1*d+3]*a^[2*d+1]*a^[3*d+2]
  264.              -a^[1*d+1]*a^[2*d+3]*a^[3*d+2] -a^[1*d+2]*a^[2*d+1]*a^[3*d+3] -a^[1*d+3]*a^[2*d+2]*a^[3*d+1]
  265.          else if (d=2) then Det:=a^[1*d+1]*a^[2*d+2] - a^[1*d+2]*a^[2*d+1]
  266.          else if (d=1) then Det:=a^[1*d+1]
  267.          else writeln('ERROR: DETERMINANT OF MATRIX OF ORDER ',d);
  268.          FreeMem(minor,size);
  269.          end;
  270. }
  271. procedure MatProd;
  272.           var i,j,k,pd :integer;
  273.               sum      :float;
  274.           begin
  275.           pd:=pred(d);
  276.           for i:=0 to pd do
  277.               for j:=0 to pd do begin
  278.                   sum:=0;
  279.                   for k:=0 to pd do sum:=sum + a^[i*d+k]*b^[k*d+j];
  280.                   c^[i*d+j]:=sum;
  281.                   end;
  282.           end;
  283.  
  284. procedure Inverse;
  285.           var determinant    :float;
  286.               minor          :matrix;
  287.               ai             :matrix;
  288.               plusminus,i,j,ii,jj,d1,d2          :integer;
  289.               flsize,sqsize,dsize,d1size,element :word;
  290.               ki                                 :char;
  291.  
  292.           begin
  293.           determinant:= det(a,d);
  294.           if (determinant=0) then writeln('Error: inverse of singular matrix.')
  295.           else begin
  296.             flsize:=sizeof(float);
  297.             dsize:=d*flsize;          {1 row of matrix}
  298.             d1size:=dsize - flsize;   {1 row of submatrix}
  299.             sqsize:=d1*d1size;        {Entire submatrix}
  300.             GetMem(minor,sqsize);
  301.             d1:=pred(d);
  302.             d2:=pred(d1);
  303.             for i:=0 to d1 do
  304.                 for j:=0 to d1 do begin
  305.                {First load minor^  with the minor of [i,j]}
  306.                    for ii:=0 to d2 do
  307.                        for jj:=0 to d2 do begin
  308.                            element:=ii*d+jj;
  309.                            if (ii>=i) then element:=element+d;
  310.                            if (jj>=j) then inc(element);
  311.                            minor^[ii*d1+jj]:=a^[element];
  312.                            end;
  313. {writeln('Minor of ',i,',',j,':'); MatPrint(minor,d1); writeln; ki:=readkey;}
  314.               {Then take the determinant of minor^.  Note that b^[i,j] is transpose
  315.                of position [j,i] of which determinant is taken.}
  316.                if ((i xor j) and 1 =0) then plusminus:=1 else plusminus:=-1;
  317.                b^[j*d+i]:=plusminus*Det(minor,d1) / determinant;
  318.                end;
  319.             FreeMem(minor,sqsize);
  320.             end;
  321.           end;
  322.  
  323. procedure Solve;
  324.           var  aa      :matrix;
  325.                i,j     :integer;
  326.                det0    :float;
  327.                size,d1 :word;
  328.           begin
  329.           det0:=Det(a,d);
  330.           if (det0=0) then writeln('ERROR: Singular argument in SOLVE.') else begin
  331.              size:=sqr(d)*sizeof(float);
  332.              GetMem(aa,size);
  333.              d1:=pred(d);
  334.              for i:=0 to d1 do begin
  335.                  move(a^,aa^,size);
  336.                  for j:=0 to d1 do aa^[j*d+i]:=b^[j];
  337.                  c^[i]:=Det(aa,d)/det0;
  338.                  end;
  339.              FreeMem(aa,size);
  340.              end;
  341.           end;
  342.  
  343. procedure MatQuot;
  344.           var bi  :matrix;
  345.              i,j  :integer;
  346.              size :word;
  347.           begin
  348.           size:=sqr(d)*sizeof(float);
  349.           GetMem(bi,size);
  350.           Inverse(b,bi,d);
  351.           MatProd(bi,a,c,d);
  352.           FreeMem(bi,size);
  353.           end;
  354.  
  355. end. {of MATRIX unit}