home *** CD-ROM | disk | FTP | other *** search
/ Math Solutions 1995 October / Math_Solutions_CD-ROM_Walnut_Creek_October_1995.iso / pc / mac / discrete / lib / matrix.g < prev    next >
Encoding:
Text File  |  1993-05-05  |  39.3 KB  |  1,443 lines
  1. #############################################################################
  2. ##
  3. #A  matrix.g                    GAP library                  Martin Schoenert
  4. ##
  5. #A  @(#)$Id: matrix.g,v 3.22 1993/02/16 12:39:55 sam Rel $
  6. ##
  7. #Y  Copyright 1990-1992,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
  8. ##
  9. ##  This file contains  those  functions  that  mainly  deal  with  matrices.
  10. ##
  11. #H  $Log: matrix.g,v $
  12. #H  Revision 3.22  1993/02/16  12:39:55  sam
  13. #H  allowed '[]' as arguments in 'TriangulizeMat' etc.
  14. #H  (although '[]' is not a matrix in the sense of GAP)
  15. #H
  16. #H  Revision 3.21  1993/01/22  12:10:53  fceller
  17. #H  fixed 'NullMat' and 'IdentityMat' to allow ring argument
  18. #H
  19. #H  Revision 3.20  1993/01/20  18:25:25  sam
  20. #H  moved added functions to directory grp ...
  21. #H
  22. #H  Revision 3.19  1993/01/18  17:18:58  sam
  23. #H  added GL, SL, SP, GU, SU in 'MatricesOps'
  24. #H
  25. #H  Revision 3.18  1992/12/16  19:47:27  martin
  26. #H  replaced quoted record names with escaped ones
  27. #H
  28. #H  Revision 3.17  1992/11/25  13:41:51  fceller
  29. #H  improved 'MinimalPolynomial' (this is due to CRLG and EOB)
  30. #H
  31. #H  Revision 3.16  1992/10/19  14:08:50  fceller
  32. #H  renamed 'WeightVecFFE' to 'DepthVector'
  33. #H
  34. #H  Revision 3.15  1992/09/10  11:14:47  martin
  35. #H  changed 'DiagonalizeMat' to use a better pivoting strategy
  36. #H
  37. #H  Revision 3.14  1992/07/02  11:48:20  fceller
  38. #H  improved 'FiniteFieldMatricesOps.*Polynomial'
  39. #H
  40. #H  Revision 3.13  1992/07/01  11:39:07  fceller
  41. #H  added 'FiniteFieldMatricesOps.MinimalPolynomial' and
  42. #H  'FiniteFieldMatricesOps.CharacteristicPolynomial'
  43. #H
  44. #H  Revision 3.12  1992/06/05  09:53:08  fceller
  45. #H  added dispatchers for minimal and characteristic polynomial
  46. #H
  47. #H  Revision 3.11  1992/05/07  11:03:15  fceller
  48. #H  added 'OrderScalar'
  49. #H
  50. #H  Revision 3.10  1992/04/23  10:34:46  fceller
  51. #H  added functions for finite field matrices
  52. #H
  53. #H  Revision 3.9  1992/04/07  16:15:32  jmnich
  54. #H  adapted to changes in the finite field module
  55. #H
  56. #H  Revision 3.8  1992/04/03  16:25:56  martin
  57. #H  fixed 'MatricesOps.in'
  58. #H
  59. #H  Revision 3.7  1992/02/19  12:36:28  martin
  60. #H  fixed another bug in 'NullspaceMat'
  61. #H
  62. #H  Revision 3.6  1992/01/29  14:58:42  sam
  63. #H  changed 'TransposedMat' (does no longer abuse 'Add')
  64. #H
  65. #H  Revision 3.5  1992/01/16  12:57:00  martin
  66. #H  changed 'Matrices' to inherit from 'GroupElements'
  67. #H
  68. #H  Revision 3.4  1991/12/04  13:25:13  martin
  69. #H  added 'KroneckerProduct'
  70. #H
  71. #H  Revision 3.3  1991/11/08  16:13:41  martin
  72. #H  changed 'matrix.g' to read 'matgrp.g' and 'matring.g' automatically
  73. #H
  74. #H  Revision 3.2  1991/11/08  15:56:11  martin
  75. #H  changed everything for new domain concept
  76. #H
  77. #H  Revision 3.1  1991/06/03  10:47:23  fceller
  78. #H  added 'Solution'
  79. #H
  80. #H  Revision 3.0  1991/05/31  13:48:28  martin
  81. #H  initial revision under RCS
  82. ##
  83.  
  84.  
  85. #############################################################################
  86. ##
  87. #F  InfoMatrix1(...)  . . . . . . . . . . . information function for matrices
  88. #F  InfoMatrix2(...)  . . . . . . . . . . . information function for matrices
  89. ##
  90. if not IsBound(InfoMatrix1)  then InfoMatrix1 := Ignore;  fi;
  91. if not IsBound(InfoMatrix2)  then InfoMatrix2 := Ignore;  fi;
  92.  
  93.  
  94. #############################################################################
  95. ##
  96. #V  Matrices  . . . . . . . . . . . . . . . . . . . .  domain of all matrices
  97. #V  MatricesOps . . . . . . .  operation record of the domain of all matrices
  98. ##
  99. MatricesOps := Copy( GroupElementsOps );
  100.  
  101. Matrices := rec(
  102.     isDomain            := true,
  103.  
  104.     name                := "Matrices",
  105.  
  106.     isFinite            := false,
  107.     size                := "infinity",
  108.  
  109.     operations          := MatricesOps
  110. );
  111.  
  112. MatricesOps.\in := function ( obj, Matrices )
  113.     return IsMat( obj );
  114. end;
  115.  
  116.  
  117. #############################################################################
  118. ##
  119. #F  MatricesOps.SpecialLinearGroup( Matrices, <n>, <q> )
  120. #F  MatricesOps.GeneralLinearGroup( Matrices, <n>, <q> )
  121. #F  MatricesOps.SymplecticGroup( Matrices, <n>, <q> )
  122. #F  MatricesOps.GeneralUnitaryGroup( Matrices, <n>, <q> )
  123. #F  MatricesOps.SpecialUnitaryGroup( Matrices, <n>, <q> )
  124. ##
  125. MatricesOps.SpecialLinearGroup := function ( M, n, q )
  126.     return MatGroupLib( "SpecialLinearGroup", n, q );
  127. end;
  128.  
  129. MatricesOps.GeneralLinearGroup := function ( M, n, q )
  130.     return MatGroupLib( "GeneralLinearGroup", n, q );
  131. end;
  132.  
  133. MatricesOps.SymplecticGroup := function ( M, n, q )
  134.     return MatGroupLib( "SymplecticGroup", n, q );
  135. end;
  136.  
  137. MatricesOps.GeneralUnitaryGroup := function ( M, n, q )
  138.     return MatGroupLib( "GeneralUnitaryGroup", n, q );
  139. end;
  140.  
  141. MatricesOps.SpecialUnitaryGroup := function ( M, n, q )
  142.     return MatGroupLib( "SpecialUnitaryGroup", n, q );
  143. end;
  144.  
  145.  
  146. #############################################################################
  147. ##
  148. #F  DimensionsMat( <mat> )  . . . . . . . . . . . . .  dimensions of a matrix
  149. ##
  150. DimensionsMat := function ( mat )
  151.     return [ Length(mat), Length(mat[1]) ];
  152. end;
  153.  
  154.  
  155. #############################################################################
  156. ##
  157. #F  IdentityMat( <m> [, <F>] )  . . . . . . . identity matrix of a given size
  158. ##
  159. IdentityMat := function ( arg )
  160.     local   id, m, zero, one, row, i, k;
  161.  
  162.     # check the arguments and get dimension, zero and one
  163.     if Length(arg) = 1  then
  164.         m    := arg[1];
  165.         zero := 0;
  166.         one  := 1;
  167.     elif Length(arg) = 2  and IsField(arg[2])  then
  168.         m    := arg[1];
  169.         zero := arg[2].zero;
  170.         one  := arg[2].one;
  171.     elif Length(arg) = 2  and IsRing(arg[2])  then
  172.         if not IsBound(arg[2].one)  then
  173.             Error( "ring must be a ring-with-one" );
  174.         fi;
  175.         m    := arg[1];
  176.         zero := arg[2].zero;
  177.         one  := arg[2].one;
  178.     elif Length(arg) = 2  then
  179.         m    := arg[1];
  180.         zero := arg[2] - arg[2];
  181.         one  := arg[2] ^ 0;
  182.     else
  183.         Error("usage: IdentityMat( <m> [, <F>] )");
  184.     fi;
  185.  
  186.     # make an empty row
  187.     row := [];
  188.     for k  in [1..m]  do row[k] := zero;  od;
  189.  
  190.     # make the identity matrix
  191.     id := [];
  192.     for i  in [1..m]  do
  193.         id[i] := ShallowCopy( row );
  194.         id[i][i] := one;
  195.     od;
  196.  
  197.     # return the identity matrix
  198.     return id;
  199. end;
  200.  
  201.  
  202. #############################################################################
  203. ##
  204. #F  NullMat( <m>, <n> [, <F>] ) . . . . . . . . . null matrix of a given size
  205. ##
  206. NullMat := function ( arg )
  207.     local   null, m, n, zero, row, i, k;
  208.  
  209.     if Length(arg) = 2  then
  210.         m    := arg[1];
  211.         n    := arg[2];
  212.         zero := 0;
  213.     elif Length(arg) = 3  and IsField(arg[3])  then
  214.         m    := arg[1];
  215.         n    := arg[2];
  216.         zero := arg[3].zero;
  217.     elif Length(arg) = 3  and IsRing(arg[3])  then
  218.         m    := arg[1];
  219.         n    := arg[2];
  220.         zero := arg[3].zero;
  221.     elif Length(arg) = 3  then
  222.         m    := arg[1];
  223.         n    := arg[2];
  224.         zero := arg[3] - arg[3];
  225.     else
  226.         Error("usage: NullMat( <m>, <n> [, <F>] )");
  227.     fi;
  228.  
  229.     # make an empty row
  230.     row := [];
  231.     for k  in [1..n]  do row[k] := zero;  od;
  232.  
  233.     # make the null matrix
  234.     null := [];
  235.     for i  in [1..m]  do
  236.         null[i] := ShallowCopy( row );
  237.     od;
  238.  
  239.     # return the null matrix
  240.     return null;
  241. end;
  242.  
  243.  
  244. #############################################################################
  245. ##
  246. #F  RandomMat( <m>, <n> [, <R>] ) . . . . . . . . . . .  make a random matrix
  247. ##
  248. ##  'RandomMat' returns a random matrix with <m> rows and  <n>  columns  with
  249. ##  elements taken from the ring <R>, which defaults to 'Integers'.
  250. ##
  251. RandomMat := function ( arg )
  252.     local   mat, m, n, R, i, row, k;
  253.  
  254.     # check the arguments and get the list of elements
  255.     if Number(arg) = 2  then
  256.         m := arg[1];
  257.         n := arg[2];
  258.         R := Integers;
  259.     elif Number(arg) = 3  then
  260.         m := arg[1];
  261.         n := arg[2];
  262.         R := arg[3];
  263.     else
  264.         Error("usage: RandomMat( <m>, <n> [, <F>] )");
  265.     fi;
  266.  
  267.     # now construct the random matrix
  268.     mat := [];
  269.     for i  in [1..m]  do
  270.         row := [];
  271.         for k  in [1..n]  do
  272.             row[k] := Random( R );
  273.         od;
  274.         mat[i] := row;
  275.     od;
  276.  
  277.     return mat;
  278. end;
  279.  
  280.  
  281. #############################################################################
  282. ##
  283. #F  RandomInvertableMat( <m> [, <R>] )  . . . make a random invertable matrix
  284. ##
  285. ##  'RandomInvertableMat' returns a invertable   random matrix with  <m> rows
  286. ##  and columns  with elements  taken from  the  ring <R>, which defaults  to
  287. ##  'Integers'.
  288. ##
  289. RandomInvertableMat := function ( arg )
  290.     local   mat, m, R, i, row, k;
  291.  
  292.     # check the arguments and get the list of elements
  293.     if Number(arg) = 1  then
  294.         m := arg[1];
  295.         R := Integers;
  296.     elif Length(arg) = 2  then
  297.         m := arg[1];
  298.         R := arg[2];
  299.     else
  300.         Error("usage: RandomInvertableMat( <m> [, <R>] )");
  301.     fi;
  302.  
  303.     # now construct the random matrix
  304.     mat := [];
  305.     for i  in [1..m]  do
  306.         repeat
  307.             row := [];
  308.             for k  in [1..m]  do
  309.                 row[k] := Random( R );
  310.             od;
  311.             mat[i] := row;
  312.         until NullspaceMat( mat ) = [];
  313.     od;
  314.  
  315.     return mat;
  316. end;
  317.  
  318.  
  319. #############################################################################
  320. ##
  321. #F  RandomUnimodularMat( <m> )  . . . . . . . . . .  random unimodular matrix
  322. ##
  323. RandomUnimodularMat := function ( m )
  324.     local  mat, c, i, j, k, l, a, b, v, w, gcd;
  325.  
  326.     # start with the identity matrix
  327.     mat := IdentityMat( m );
  328.  
  329.     for c  in [1..m]  do
  330.  
  331.         # multiply two random rows with a random? unimodular 2x2 matrix
  332.         i := RandomList([1..m])+1;
  333.         repeat
  334.             j := RandomList([1..m])+1;
  335.         until j <> i;
  336.         repeat
  337.             a := Random( Integers );  b := Random( Integers );
  338.             gcd := Gcdex( a, b );
  339.         until gcd.gcd = 1;
  340.         v := mat[i];  w := mat[j];
  341.         mat[i] := gcd.coeff1 * v + gcd.coeff2 * w;
  342.         mat[j] := gcd.coeff3 * v + gcd.coeff4 * w;
  343.  
  344.         # multiply two random cols with a random? unimodular 2x2 matrix
  345.         k := RandomList([1..m])+1;
  346.         repeat
  347.             l := RandomList([1..m])+1;
  348.         until l <> k;
  349.         repeat
  350.             a := Random( Integers );  b := Random( Integers );
  351.             gcd := Gcdex( a, b );
  352.         until gcd.gcd = 1;
  353.         for i  in [1..m]  do
  354.             v := mat[i][k];  w := mat[i][l];
  355.             mat[i][k] := gcd.coeff1 * v + gcd.coeff2 * w;
  356.             mat[i][l] := gcd.coeff3 * v + gcd.coeff4 * w;
  357.         od;
  358.  
  359.     od;
  360.  
  361.     return mat;
  362. end;
  363.  
  364.  
  365. #############################################################################
  366. ##
  367. #F  TransposedMat( <mat> )  . . . . . . . . . . . . .  transposed of a matrix
  368. ##
  369. ##  'TransposedMat' returns the transposed of the matrix  <mat>,  i.e., a new
  370. ##  matrix <trans> such that '<trans>[<i>][<k>]' is '<mat>[<k>][<i>]'.
  371. ##
  372. TransposedMat := function ( mat )
  373.     local  trn, col, n, m, i, j;
  374.  
  375.     if mat = [] then return []; fi;
  376.  
  377.     # initialize the transposed
  378.     n := Length( mat[1] );
  379.     m := Length( mat );
  380.     trn := [];
  381.  
  382.     # copy the entries
  383.     for j  in [ 1 .. n ]  do
  384.         col := [];
  385.         for i  in [ 1 .. m ]  do
  386.             col[i] := mat[i][j];
  387.         od;
  388.         trn[j] := col;
  389.     od;
  390.  
  391.     # return the transposed
  392.     return trn;
  393. end;
  394.  
  395.  
  396. #############################################################################
  397. ##
  398. #F  KroneckerProduct(<mat1>,<mat2>) . . . . Kronecker product of two matrices
  399. ##
  400. KroneckerProduct := function ( mat1, mat2 )
  401.     local i, row1, row2, row, kroneckerproduct;
  402.     kroneckerproduct := [];
  403.     for row1  in mat1  do
  404.         for row2  in mat2  do
  405.             row := [];
  406.             for i  in row1  do
  407.                 Append( row, i * row2 );
  408.             od;
  409.             Add( kroneckerproduct, row );
  410.         od;
  411.     od;
  412.     return kroneckerproduct;
  413. end;
  414.  
  415.  
  416. #############################################################################
  417. ##
  418. #F  OrderMat( <mat> ) . . . . . . . . . . . . . . . . . . . order of a matrix
  419. #F  MatricesOps.Order(<Matrices>,<mat>) . . . . . . . . . . order of a matrix
  420. ##
  421. OrderMatLimit := 1000;
  422.  
  423. OrderMat := function ( mat )
  424.     local   ord, m, id, vec, v, o, i;
  425.  
  426.     # check that the argument is an invertable square matrix
  427.     m := Length(mat);
  428.     if m <> Length(mat[1])  then
  429.         Error("OrderMat: <mat> must be a square matrix");
  430.     fi;
  431.     if RankMat( mat ) <> m  then
  432.         Error("OrderMat: <mat> must be invertable");
  433.     fi;
  434.     id := mat ^ 0;
  435.  
  436.     # loop over the standard basis vectors
  437.     ord := 1;
  438.     for i  in [1..m]  do
  439.  
  440.         # compute the length of the orbit of the <i>th standard basis vector
  441.         vec := mat[i];
  442.         v   := vec * mat;
  443.         o   := 1;
  444.         while v <> vec  do
  445.             v := v * mat;
  446.             o := o + 1;
  447.             if OrderMatLimit = o  then
  448.                 Print("#W  OrderMat: warning, order be infinite\n");
  449.             fi;
  450.         od;
  451.  
  452.         # raise the matrix to this length (new mat will fix basis vector)
  453.         mat := mat ^ o;
  454.         ord := ord * o;
  455.     od;
  456.  
  457.     # return the order
  458.     return ord;
  459. end;
  460.  
  461. MatricesOps.Order := function ( Matrices, mat )
  462.     return OrderMat( mat );
  463. end;
  464.  
  465.  
  466. #############################################################################
  467. ##
  468. #F  TraceMat( <mat> ) . . . . . . . . . . . . . . . . . . . trace of a matrix
  469. ##
  470. TraceMat := function ( mat )
  471.     local   trc, m, i;
  472.  
  473.     # check that the element is a square matrix
  474.     m := Length(mat);
  475.     if m <> Length(mat[1])  then
  476.         Error("TraceMat: <mat> must be a square matrix");
  477.     fi;
  478.  
  479.     # sum all the diagonal entries
  480.     trc := mat[1][1];
  481.     for i  in [2..m]  do
  482.         trc := trc + mat[i][i];
  483.     od;
  484.  
  485.     # return the trace
  486.     return trc;
  487. end;
  488.  
  489.  
  490. #############################################################################
  491. ##
  492. #F  RankMat( <mat> )  . . . . . . . . . . . . . . . . . . .  rank of a matrix
  493. ##
  494. RankMat := function ( mat )
  495.     local   m, n, i, j, k, row, zero;
  496.  
  497.     # make a copy to avoid changing the original argument
  498.     m := Length(mat);
  499.     n := Length(mat[1]);
  500.     zero := 0*mat[1][1];
  501.     mat := ShallowCopy( mat );
  502.     InfoMatrix1("#I  RankMat called\n");
  503.     InfoMatrix2("#I    nonzero columns: \c");
  504.  
  505.     # run through all columns of the matrix
  506.     i := 0;
  507.     for k  in [1..n]  do
  508.  
  509.         # find a nonzero entry in this column
  510.         #N  26-Oct-91 martin if <mat> is an rational matrix look for a pivot
  511.         j := i + 1;
  512.         while j <= m and mat[j][k] = zero  do j := j+1;  od;
  513.  
  514.         # if there is a nonzero entry
  515.         if j <= m  then
  516.  
  517.             # increment the rank
  518.             InfoMatrix2(k," \c");
  519.             i := i + 1;
  520.  
  521.             # make its row the current row and normalize it
  522.             row := mat[j];  mat[j] := mat[i];  mat[i] := row[k]^-1 * row;
  523.  
  524.             # clear all entries in this column
  525.             for j  in [i+1..m] do
  526.                 if  mat[j][k] <> zero  then
  527.                     mat[j] := mat[j] - mat[j][k] * mat[i];
  528.                 fi;
  529.             od;
  530.  
  531.         fi;
  532.  
  533.     od;
  534.  
  535.     # return the rank (the number of linear independent rows)
  536.     InfoMatrix2("\n");
  537.     InfoMatrix1("#I  RankMat returns ",i,"\n");
  538.     return i;
  539. end;
  540.  
  541.  
  542. #############################################################################
  543. ##
  544. #F  DeterminantMat( <mat> ) . . . . . . . . . . . . . determinant of a matrix
  545. ##
  546. DeterminantMat := function ( mat )
  547.     local   det, sgn, row, zero, m, i, j, k, l;
  548.  
  549.     # check that the argument is a square matrix and get the size
  550.     m := Length(mat);
  551.     zero := 0*mat[1][1];
  552.     if m <> Length(mat[1])  then
  553.         Error("DeterminantMat: <mat> must be a square matrix");
  554.     fi;
  555.  
  556.     # make a copy to avoid changing the orginal argument
  557.     mat := ShallowCopy( mat );
  558.     InfoMatrix1("#I  DeterminantMat called\n");
  559.     InfoMatrix2("#I    nonzero columns: \c");
  560.  
  561.     # run through all columns of the matrix
  562.     i := 0;  det := 1;  sgn := 1;
  563.     for k  in [1..m]  do
  564.  
  565.         # find a nonzero entry in this column
  566.         #N  26-Oct-91 martin if <mat> is an rational matrix look for a pivot
  567.         j := i + 1;
  568.         while j <= m and mat[j][k] = zero  do j := j+1;  od;
  569.  
  570.         # if there is a nonzero entry
  571.         if j <= m  then
  572.  
  573.             # increment the rank
  574.             InfoMatrix2(k," \c");
  575.             i := i + 1;
  576.  
  577.             # make its row the current row
  578.             if i <> j  then
  579.                 row := mat[j];  mat[j] := mat[i];  mat[i] := row;
  580.                 sgn := -sgn;
  581.             fi;
  582.  
  583.             # clear all entries in this column
  584.             for j  in [i+1..m]  do
  585.                 mat[j] := (mat[i][k]*mat[j]-mat[j][k]*mat[i])/det;
  586.             od;
  587.  
  588.             det := mat[i][k];
  589.         fi;
  590.  
  591.     od;
  592.  
  593.     # return the determinant
  594.     if i < m  then det := zero;  fi;
  595.     InfoMatrix2("\n");
  596.     InfoMatrix1("#I  DeterminantMat returns ",sgn*det,"\n");
  597.     return sgn * det;
  598. end;
  599.  
  600.  
  601. #############################################################################
  602. ##
  603. #F  TriangulizeMat( <mat> ) . . . . . bring a matrix in upper triangular form
  604. ##
  605. TriangulizeMat := function ( mat )
  606.     local   m, n, i, j, k, row, zero;
  607.  
  608.     InfoMatrix1("#I  TriangulizeMat called\n");
  609.  
  610.     if mat <> [] then 
  611.  
  612.        # get the size of the matrix
  613.        m := Length(mat);
  614.        n := Length(mat[1]);
  615.        zero := 0*mat[1][1];
  616.        InfoMatrix2("#I    nonzero columns: \c");
  617.    
  618.        # run through all columns of the matrix
  619.        i := 0;
  620.        for k  in [1..n]  do
  621.    
  622.            # find a nonzero entry in this column
  623.            j := i + 1;
  624.            while j <= m and mat[j][k] = zero  do j := j + 1;  od;
  625.    
  626.            # if there is a nonzero entry
  627.            if j <= m  then
  628.    
  629.                # increment the rank
  630.                InfoMatrix2(k," \c");
  631.                i := i + 1;
  632.    
  633.                # make its row the current row and normalize it
  634.                row := mat[j];  mat[j] := mat[i];  mat[i] := row[k]^-1 * row;
  635.    
  636.                # clear all entries in this column
  637.                for j  in [1..m] do
  638.                    if  i <> j  and mat[j][k] <> zero  then
  639.                        mat[j] := mat[j] - mat[j][k] * mat[i];
  640.                    fi;
  641.                od;
  642.    
  643.            fi;
  644.    
  645.        od;
  646.  
  647.     fi;
  648.  
  649.     InfoMatrix2("\n");
  650.     InfoMatrix1("#I  TriangulizeMat returns\n");
  651. end;
  652.  
  653.  
  654. #############################################################################
  655. ##
  656. #F  BaseMat( <mat> )  . . . . . . . . . .  base for the row space of a matrix
  657. ##
  658. BaseMat := function ( mat )
  659.     local  base, m, n, i, k, zero;
  660.  
  661.     base := [];
  662.  
  663.     if mat <> [] then
  664.  
  665.        # make a copy to avoid changing the original argument
  666.        mat := ShallowCopy( mat );
  667.        m := Length(mat);
  668.        n := Length(mat[1]);
  669.        zero := 0*mat[1][1];
  670.    
  671.        # triangulize the matrix
  672.        TriangulizeMat( mat );
  673.    
  674.        # and keep only the nonzero rows of the triangular matrix
  675.        i := 1;
  676.        for k  in [1..n]  do
  677.            if i <= m  and mat[i][k] <> zero  then
  678.                Add( base, mat[i] );
  679.                i := i + 1;
  680.            fi;
  681.        od;
  682.  
  683.     fi;
  684.  
  685.     # return the base
  686.     return base;
  687. end;
  688.  
  689.  
  690. #############################################################################
  691. ##
  692. #F  NullspaceMat( <mat> ) . . . . . . basis of solutions of <vec> * <mat> = 0
  693. ##
  694. NullspaceMat := function ( mat )
  695.     local   nullspace, m, n, min, empty, i, k, row, tmp, zero, one;
  696.  
  697.     # triangulize the transposed of the matrix
  698.     mat := TransposedMat( Reversed( mat ) );
  699.     TriangulizeMat( mat );
  700.     m := Length(mat);
  701.     n := Length(mat[1]);
  702.     zero := 0*mat[1][1];
  703.     one  := mat[1][1]^0;
  704.     min := Minimum( m, n );
  705.  
  706.     # insert empty rows to bring the leading term of each row on the diagonal
  707.     empty := 0*mat[1];
  708.     i := 1;
  709.     while i <= Length(mat)  do
  710.         if i < n  and mat[i][i] = zero  then
  711.             for k in Reversed([i..Minimum(Length(mat),n-1)])  do
  712.                 mat[k+1] := mat[k];
  713.             od;
  714.             mat[i] := empty;
  715.         fi;
  716.         i := i+1;
  717.     od;
  718.     for i  in [ Length(mat)+1 .. n ]  do
  719.         mat[i] := empty;
  720.     od;
  721.  
  722.     # 'mat' now  looks  like  [ [1,2,0,2], [0,0,0,0], [0,0,1,3], [0,0,0,0] ],
  723.     # and the solutions can be read in those columns with a 0 on the diagonal
  724.     # by replacing this 0 by a -1, in  this  example  [2,-1,0,0], [2,0,3,-1].
  725.     nullspace := [];
  726.     for k  in Reversed([1..n]) do
  727.         if mat[k][k] = zero  then
  728.             row := [];
  729.             for i  in [1..k-1]  do row[n-i+1] := -mat[i][k];  od;
  730.             row[n-k+1] := one;
  731.             for i  in [k+1..n]  do row[n-i+1] := zero;  od;
  732.             Add( nullspace, row );
  733.         fi;
  734.     od;
  735.  
  736.     return nullspace;
  737. end;
  738.  
  739.  
  740. #############################################################################
  741. ##
  742. #F  DiagonalizeMat( <mat> ) . . . . . . . . . . diagonalize an integer matrix
  743. ##
  744. DiagonalizeMat := function ( mat )
  745.     local   i,          # current position
  746.             k, l,       # row and column loop variables
  747.             r, c,       # row and column index of minimal element
  748.             e, f,       # entries of the matrix
  749.             g,          # (extended) gcd of <e> and <f>
  750.             v, w,       # rows of the matrix or row and column norms
  751.             m,          # maximal entry, only for information
  752.             isClearCol; #
  753.  
  754.     InfoMatrix1("#I  DiagonalizeMat called\n");
  755.  
  756.     if mat <> [] then
  757.  
  758.        InfoMatrix2("#I    divisors: \c");
  759.        m := 0;
  760.    
  761.        # loop over all rows respectively columns of the matrix
  762.        for i  in [1..Minimum( Length(mat), Length(mat[1]) )]  do
  763.    
  764.            # compute the row and column norms
  765.            v := 0 * [1..Length(mat)];
  766.            w := 0 * [1..Length(mat[1])];
  767.            for k  in [i..Length(mat)]  do
  768.                for l  in [i..Length(mat[1])]  do
  769.                    e := mat[k][l];
  770.                    if   0 < e  then
  771.                        v[k] := v[k] + e;
  772.                        w[l] := w[l] + e;
  773.                    else
  774.                        v[k] := v[k] - e;
  775.                        w[l] := w[l] - e;
  776.                    fi;
  777.                od;
  778.            od;
  779.    
  780.            # find the element with the smallest absolute value in the matrix
  781.            f := 0;
  782.            for k  in [i..Length(mat)]  do
  783.                for l  in [i..Length(mat[1])]  do
  784.                    e := mat[k][l];
  785.                    if   0 < e and (f=0 or  e<f or  e=f and v[k]*w[l]<g) then
  786.                        f :=  e;  r := k;  c := l;  g := v[k]*w[l];
  787.                    elif e < 0 and (f=0 or -e<f or -e=f and v[k]*w[l]<g) then
  788.                        f := -e;  r := k;  c := l;  g := v[k]*w[l];
  789.                    fi;
  790.                    if   0 < e and m <  e  then
  791.                        m :=  e;
  792.                    elif e < 0 and m < -e  then
  793.                        m := -e;
  794.                    fi;
  795.                od;
  796.            od;
  797.    
  798.            # if there is no nonzero entry we are done
  799.            if f = 0  then
  800.                InfoMatrix2("\n");
  801.                InfoMatrix1("#I  DiagonalizeMat returns\n");
  802.                return;
  803.            fi;
  804.    
  805.            # move the minimal element to position 'mat[i][i]' make it positive
  806.            if i <> r  then
  807.                v := mat[i];  mat[i] := mat[r];  mat[r] := v;
  808.            fi;
  809.            if i <> c  then
  810.                for k  in [i..Length(mat)]  do
  811.                    e := mat[k][i];  mat[k][i] := mat[k][c];  mat[k][c] := e;
  812.                od;
  813.            fi;
  814.            if mat[i][i] < 0  then
  815.                mat[i] := - mat[i];
  816.            fi;
  817.    
  818.            # now clear the column i and the row i
  819.            isClearCol := false;
  820.            while not isClearCol  do
  821.    
  822.                # clear the column i using unimodular row operations such that
  823.                # mat[i][i] becomes gcd(mat[i][i],mat[r][i]) and mat[r][i] = 0
  824.                k := i + 1;
  825.                while k <= Length(mat)  do
  826.                    e := mat[i][i];  f := mat[k][i];
  827.                    if f mod e = 0  then
  828.                        mat[k] := mat[k] - f/e * mat[i];
  829.                    elif f <> 0  then
  830.                        g := Gcdex( e, f );
  831.                        v := mat[i];  w := mat[k];
  832.                        mat[i] := g.coeff1 * v + g.coeff2 * w;
  833.                        mat[k] := g.coeff3 * v + g.coeff4 * w;
  834.                    fi;
  835.                    k := k + 1;
  836.                od;
  837.    
  838.                isClearCol := true;
  839.    
  840.                # clear the row i using unimodular column operations such that
  841.                # mat[i][i] becomes gcd(mat[i][i],mat[i][c]) and mat[i][c] = 0
  842.                # after such an operation we may have to clear column i  again
  843.                l := i + 1;
  844.                while l <= Length(mat[1])  and isClearCol  do
  845.                    e := mat[i][i];  f := mat[i][l];
  846.                    if f mod e = 0  then
  847.                        mat[i][l] := 0;
  848.                    elif f <> 0  then
  849.                        g := Gcdex( e, f );
  850.                        for k  in [i..Length(mat)]  do
  851.                            v := mat[k][i];  w := mat[k][l];
  852.                            mat[k][i] := g.coeff1 * v + g.coeff2 * w;
  853.                            mat[k][l] := g.coeff3 * v + g.coeff4 * w;
  854.                        od;
  855.                        isClearCol := false;
  856.                    fi;
  857.                    l := l + 1;
  858.                od;
  859.    
  860.            od;
  861.    
  862.            InfoMatrix2(mat[i][i]," \c");
  863.        od;
  864.    
  865.        InfoMatrix2("\n");
  866.        InfoMatrix2("#I  maximal entry ",m,"\n");
  867.     fi;
  868.  
  869.     InfoMatrix1("#I  DiagonalizeMat returns\n");
  870. end;
  871.  
  872.  
  873. #############################################################################
  874. ##
  875. #F  ElementaryDivisorsMat( <mat> )   elementary divisors of an integer matrix
  876. ##
  877. ##  'ElementaryDivisors' returns a list of the elementary divisors, i.e., the
  878. ##  unique <d> with '<d>[<i>]' divides '<d>[<i>+1]' and <mat>  is  equivalent
  879. ##  to a diagonal matrix with the elements '<d>[<i>]' on the diagonal.
  880. ##
  881. ElementaryDivisorsMat := function ( mat )
  882.     local  divs, gcd, m, n, i, k;
  883.  
  884.     # make a copy to avoid changing the original argument
  885.     mat := Copy( mat );
  886.     m := Length(mat);  n := Length(mat[1]);
  887.  
  888.     # diagonalize the matrix
  889.     DiagonalizeMat( mat );
  890.  
  891.     # get the diagonal elements
  892.     divs := [];
  893.     for i  in [1..Minimum(m,n)]  do
  894.         divs[i] := mat[i][i];
  895.     od;
  896.  
  897.     # transform the divisors so that every divisor divides the next
  898.     for i  in [1..Length(divs)-1]  do
  899.         for k  in [i+1..Length(divs)]  do
  900.             if divs[i] <> 0  and divs[k] mod divs[i] <> 0  then
  901.                 gcd     := GcdInt( divs[i], divs[k] );
  902.                 divs[k] := divs[k] / gcd * divs[i];
  903.                 divs[i] := gcd;
  904.             fi;
  905.         od;
  906.     od;
  907.  
  908.     return divs;
  909. end;
  910.  
  911.  
  912. #############################################################################
  913. ##
  914. #F  SolutionMat(<mat>,<vec>)  . . . . . . . . . . .  one solution of equation
  915. ##
  916. ##  One solution <x> of <x> * <mat> = <vec> or 'false'.
  917. ##
  918. SolutionMat := function ( mat, vec )
  919.     local   h, v, tmp, i, l, r, s, c, zero;
  920.  
  921.     # Transpose <mat> and solve <mat> * x = <vec>.
  922.     mat  := TransposedMat( mat );
  923.     vec  := Copy( vec );
  924.     l    := Length( mat );
  925.     r    := 0;
  926.     zero := 0*mat[1][1];
  927.     InfoMatrix1("#I  SolutionMat called\n");
  928.     InfoMatrix2("#I    nonzero columns \c");
  929.  
  930.     # Run through all columns of the matrix.
  931.     c := 1;
  932.     while c <= Length( mat[ 1 ] ) and r < l  do
  933.  
  934.         # Find a nonzero entry in this column.
  935.         s := r + 1;
  936.         while s <= l and mat[ s ][ c ] = zero  do s := s + 1;  od;
  937.  
  938.         # If there is a nonzero entry,
  939.         if s <= l  then
  940.  
  941.             # increment the rank.
  942.             InfoMatrix2(c," \c");
  943.             r := r + 1;
  944.  
  945.             # Make its row the current row and normalize it.
  946.             tmp := mat[ s ][ c ] ^ -1;
  947.             v := mat[ s ];  mat[ s ] := mat[ r ];  mat[ r ] := tmp * v;
  948.             v := vec[ s ];  vec[ s ] := vec[ r ];  vec[ r ] := tmp * v;
  949.  
  950.             # Clear all entries in this column.
  951.             for s  in [ 1 .. Length( mat ) ]  do
  952.                 if s <> r and mat[ s ][ c ] <> zero  then
  953.                     tmp := mat[ s ][ c ];
  954.                     mat[ s ] := mat[ s ] - tmp * mat[ r ];
  955.                     vec[ s ] := vec[ s ] - tmp * vec[ r ];
  956.                 fi;
  957.             od;
  958.         fi;
  959.         c := c + 1;
  960.     od;
  961.  
  962.     # Find a solution.
  963.     for i  in [ r + 1 .. l ]  do
  964.         if vec[ i ] <> zero  then return false;  fi;
  965.     od;
  966.     h := [];
  967.     s := Length( mat[ 1 ] );
  968.     v := 0 * mat[ 1 ][ 1 ];
  969.     r := 1;
  970.     c := 1;
  971.     while c <= s and r <= l  do
  972.         while c <= s and mat[ r ][ c ] = zero  do
  973.             c := c + 1;
  974.             Add( h, v );
  975.         od;
  976.         if c <= s  then
  977.             Add( h, vec[ r ] );
  978.             r := r + 1;
  979.             c := c + 1;
  980.         fi;
  981.     od;
  982.     while c <= s  do Add( h, v );  c := c + 1;  od;
  983.  
  984.     InfoMatrix2("\n");
  985.     InfoMatrix1("#I  SolutionMat returns\n");
  986.     return h;
  987. end;
  988.  
  989.  
  990. #############################################################################
  991. ##
  992. #R  Read  . . . . . . . . . . . . .  read other function from the other files
  993. ##
  994. ReadLib( "matgrp" );
  995. ReadLib( "matring" );
  996.  
  997.  
  998. #############################################################################
  999. ##
  1000. #V  FieldMatrices . . . . . . . . . . . . . . . .  matrices over finite field
  1001. ##
  1002. FieldMatrices       := Copy( Matrices );
  1003. FieldMatrices.name  := "FieldMatrices";
  1004. FieldMatricesOps    := FieldMatrices.operations;
  1005.  
  1006.  
  1007. #############################################################################
  1008. ##
  1009. #F  MinimalPolynomial( <D>, <A> ) . . . . . . . . . . . .  minimal polynomial
  1010. ##
  1011. MinimalPolynomial := function( arg )
  1012.     local   D,  A;
  1013.  
  1014.     # get the domain <D> of <A>
  1015.     if Length(arg) = 2  then
  1016.         D := arg[1];
  1017.         A := arg[2];
  1018.     else
  1019.         A := arg[1];
  1020.         D := Domain([A]);
  1021.     fi;
  1022.  
  1023.     # return the minimal polynomial
  1024.     return D.operations.MinimalPolynomial(D, A);
  1025.  
  1026. end;
  1027.  
  1028.  
  1029. #############################################################################
  1030. ##
  1031. #F  FieldMatricesOps.MinimalPolynomial( <D>, <A> )  . . min polynomial of <A>
  1032. ##
  1033. FieldMatricesOps.MinimalPolynomial := function( D, A )
  1034.     local   I,  F,  L,  d,  p,  P,  M,  N,  one,  R,  h,  v,  w,  i,  j,  W;
  1035.  
  1036.     # try to find the field of <A>
  1037.     if not IsBound(D.baseRing)  then
  1038.     F := Field( Flat(A) );
  1039.     else
  1040.     F := D.baseRing;
  1041.     fi;
  1042.  
  1043.     # set up the other variables
  1044.     W := [];
  1045.     d := Length( A );
  1046.     P := [ Polynomial( F, [ F.one ] ) ];
  1047.     N := F.zero * A[1];
  1048.     M := ShallowCopy( N );
  1049.     Add( M, N[1] );
  1050.  
  1051.     # span the whole vector space
  1052.     for i  in [ 1 .. d ]  do
  1053.  
  1054.     # did we already see the <i>.th base vector
  1055.         if not IsBound( W[i] )  then
  1056.  
  1057.             # clear right and left sides <R>
  1058.             L := [];
  1059.             R := [];
  1060.  
  1061.             # start with the standard base vector
  1062.             v := ShallowCopy( N );
  1063.             v[i] := F.one;
  1064.  
  1065.             # <j>-1 gives the power of <A> we are looking at
  1066.             j := 1;
  1067.  
  1068.             # spin vector around and construct polynomial
  1069.             repeat
  1070.  
  1071.                 # compute the head of <v>
  1072.                 h := 1;
  1073.                 while h <= d and v[h] = F.zero  do
  1074.                     h := h + 1;
  1075.                 od;
  1076.  
  1077.               # start with appropriate polynomial x^(<j>-1)
  1078.                 p := ShallowCopy( M );
  1079.                 p[j] := F.one;
  1080.  
  1081.                 # divide by known left sides
  1082.                 w := v;
  1083.                 while h <= d and IsBound( L[h] )  do
  1084.                   p := p - w[h] * R[h];
  1085.                     w := w - w[h] * L[h];
  1086.                     while h <= d and w[h] = F.zero  do
  1087.                         h := h + 1;
  1088.                     od;
  1089.                 od;
  1090.  
  1091.                 # if <v> is not the zero vector try next power
  1092.                 if h <= d  then
  1093.                     R[h] := p * w[h]^-1;
  1094.                     L[h] := w * w[h]^-1;
  1095.                     W[h] := true;
  1096.                     j := j + 1;
  1097.                     v := v * A;
  1098.  
  1099.                 # otherwise, we know our polynomial
  1100.                 else
  1101.                     Add( P, Polynomial( F, p ) );
  1102.                 fi;
  1103.             until h > d;
  1104.         fi;
  1105.     od;
  1106.  
  1107.     # compute LCM of all polynomials
  1108.     return Lcm( P );
  1109.  
  1110. end;
  1111.  
  1112.  
  1113. #############################################################################
  1114. ##
  1115. #F  CharacteristicPolynomial( <D>, <A> )  . . . . . characteristic polynomial
  1116. ##
  1117. CharacteristicPolynomial := function( arg )
  1118.     local   D,  A;
  1119.  
  1120.     # get the domain <D> of <A>
  1121.     if Length(arg) = 2  then
  1122.         D := arg[1];
  1123.         A := arg[2];
  1124.     else
  1125.         A := arg[1];
  1126.         D := Domain([A]);
  1127.     fi;
  1128.  
  1129.     # return the characteristic polynomial
  1130.     return D.operations.CharacteristicPolynomial(D, A);
  1131.  
  1132. end;
  1133.  
  1134.  
  1135. #############################################################################
  1136. ##
  1137. #F  FieldMatricesOps.CharacteristicPolynomial( <D>, <A> ) . . char pol of <A>
  1138. ##
  1139. FieldMatricesOps.CharacteristicPolynomial := function( D, A )
  1140.     local   F,  L,  d,  p,  P,  M,  N,  one,  R,  h,  v,  w,  i,  j;
  1141.  
  1142.     # try to find the field of <A>
  1143.     if not IsBound(D.baseRing)  then
  1144.     F := Field( Flat(A) );
  1145.     else
  1146.     F := D.baseRing;
  1147.     fi;
  1148.  
  1149.     # set up the other variables
  1150.     L := [];
  1151.     d := Length( A );
  1152.     P := [ Polynomial( F, [ F.one ] ) ];
  1153.     N := F.zero * A[1];
  1154.     M := ShallowCopy( N );
  1155.     Add( M, N[1] );
  1156.  
  1157.     # span the whole vector space
  1158.     for i  in [ 1 .. d ]  do
  1159.  
  1160.     # did we already see the <i>.th base vector
  1161.         if not IsBound( L[i] )  then
  1162.  
  1163.             # clear right sides <R>
  1164.             R := [];
  1165.  
  1166.             # start with the standard base vector
  1167.             v := ShallowCopy( N );
  1168.             v[i] := F.one;
  1169.  
  1170.             # <j>-1 gives the power of <A> we are looking at
  1171.             j := 1;
  1172.  
  1173.             # spin vector around and construct polynomial
  1174.             repeat
  1175.  
  1176.                 # compute the head of <v>
  1177.                 h := 1;
  1178.                 while h <= d and v[h] = F.zero  do
  1179.                     h := h + 1;
  1180.                 od;
  1181.  
  1182.               # start with appropriate polynomial x^(<j>-1)
  1183.                 p := ShallowCopy( M );
  1184.                 p[j] := F.one;
  1185.  
  1186.                 # divide by known left sides
  1187.                 w := v;
  1188.                 while h <= d and IsBound( L[h] )  do
  1189.                     if IsBound( R[h] )  then
  1190.                         p := p - w[h] * R[h];
  1191.                     fi;
  1192.                     w := w - w[h] * L[h];
  1193.                     while h <= d and w[h] = F.zero  do
  1194.                         h := h + 1;
  1195.                     od;
  1196.                 od;
  1197.  
  1198.                 # if <v> is not the zero vector try next power
  1199.                 if h <= d  then
  1200.                     R[h] := p * w[h]^-1;
  1201.                     L[h] := w * w[h]^-1;
  1202.                     j := j + 1;
  1203.                     v := v * A;
  1204.  
  1205.                 # otherwise, we know our polynomial
  1206.                 else
  1207.                     Add( P, Polynomial( F, p ) );
  1208.                 fi;
  1209.             until h > d;
  1210.         fi;
  1211.     od;
  1212.  
  1213.     # compute the product of all polynomials
  1214.     return Product( P );
  1215.  
  1216. end;
  1217.  
  1218.  
  1219. #############################################################################
  1220. ##
  1221. #V  FiniteFieldMatrices    . . . . . . . . . . . . . . . . finite field matrices
  1222. ##
  1223. FiniteFieldMatrices       := Copy( FieldMatrices );
  1224. FiniteFieldMatrices.name  := "FiniteFieldMatrices";
  1225. FiniteFieldMatricesOps    := FiniteFieldMatrices.operations;
  1226.  
  1227.  
  1228. #############################################################################
  1229. ##
  1230. #F  FiniteFieldMatricesOps.OrderScalar( <R>, <A> )  . . . . order mod scalars
  1231. ##
  1232. ##  Return an integer n and a finite field element e such that <A>^n = eI.
  1233. ##
  1234. FiniteFieldMatricesOps.OrderScalar := function( R, A )
  1235.     local   p,  R;
  1236.  
  1237.     # construct the minimal polynomial of <A>
  1238.     p := R.operations.MinimalPolynomial( R, A );
  1239.  
  1240.     # check if <A> is invertible
  1241.     if p.coefficients[1] = 0 * p.coefficients[1]  then
  1242.     Error( "<A> must be invertible" );
  1243.     fi;
  1244.  
  1245.     # compute the order of <p>
  1246.     R := DefaultRing( p );
  1247.     return R.operations.OrderScalar( R, p );
  1248.  
  1249. end;
  1250.  
  1251.  
  1252. #############################################################################
  1253. ##
  1254. #F  FiniteFieldMatricesOps.Order( <R>, <A> )  . . . . . . . . .  order of <A>
  1255. ##
  1256. FiniteFieldMatricesOps.Order := function( R, A )
  1257.     local   o;
  1258.  
  1259.     o := R.operations.OrderScalar( R, A );
  1260.     return o[1] * OrderFFE(o[2]);
  1261.  
  1262. end;
  1263.  
  1264.  
  1265. #############################################################################
  1266. ##
  1267. #F  FiniteFieldMatricesOps.MinimalPolynomial( <D>, <A> )  . .  min polynomial
  1268. ##
  1269. FiniteFieldMatricesOps.MinimalPolynomial := function( D, A )
  1270.     local   I,  F,  L,  d,  p,  P,  M,  N,  one,  R,  h,  v,  w,  i,  W;
  1271.  
  1272.     # try to find the field of <A>
  1273.     if not IsBound(D.baseRing)  then
  1274.     F := Field( List( A, x -> Field(x).root ) );
  1275.     else
  1276.     F := D.baseRing;
  1277.     fi;
  1278.  
  1279.     # set up the other variables
  1280.     W := [];
  1281.     d := Length(A);
  1282.     P := [ Polynomial(F, [F.one]) ];
  1283.     N := F.zero * A[1];
  1284.     IsVector(N);
  1285.  
  1286.     # span the whole vector space
  1287.     for i  in [ 1 .. d ]  do
  1288.  
  1289.     # did we already see the <i>.th base vector
  1290.         if not IsBound(W[i])  then
  1291.  
  1292.             # clear right and left sides <R>
  1293.             L := [];
  1294.             R := [];
  1295.  
  1296.             # start with the standard base vector
  1297.             v := ShallowCopy(N);
  1298.             v[i] := F.one;
  1299.  
  1300.             # <j>-1 gives the power of <A> we are looking at
  1301.             M := [];
  1302.  
  1303.             # spin vector around and construct polynomial
  1304.             repeat
  1305.  
  1306.                 # compute the head of <v>
  1307.                 h := DepthVector(v);
  1308.  
  1309.               # start with appropriate polynomial x^(<j>-1)
  1310.                 p := Copy(M);
  1311.                 Add(p, F.one);
  1312.                 Add(M, F.zero);
  1313.  
  1314.                 # divide by known left sides
  1315.                 w := Copy(v);
  1316.                 while h <= Length(w) and IsBound(L[h])  do
  1317.                 AddCoeffs(p, R[h], -w[h]);
  1318.                     AddCoeffs(w, L[h], -w[h]);
  1319.                     h := DepthVector(w);
  1320.                 od;
  1321.  
  1322.                 # if <v> is not the zero vector try next power
  1323.                 if h <= Length(w)  then
  1324.                     R[h] := p * w[h]^-1;
  1325.                     L[h] := w * w[h]^-1;
  1326.  
  1327.             # enter vectors seen in <W>
  1328.                     while h <= Length(w) and IsBound(W[h])  do
  1329.                 AddCoeffs(w, W[h], -w[h]);
  1330.                    h := DepthVector(w);
  1331.                     od;
  1332.                     if h <= Length(w) then 
  1333.                            W[h] := w * w[h]^-1;
  1334.                     fi;
  1335.  
  1336.             # next power of <A>
  1337.                     v := v * A;
  1338.             h := 1;
  1339.  
  1340.                 # otherwise, we know our polynomial
  1341.                 else
  1342.                     Add(P, Polynomial(F, p));
  1343.             h := 0;
  1344.                 fi;
  1345.             until h = 0;
  1346.         fi;
  1347.     od;
  1348.  
  1349.     # compute LCM of all polynomials
  1350.     return Lcm(P);
  1351.  
  1352. end;
  1353.  
  1354.  
  1355. #############################################################################
  1356. ##
  1357. #F  FiniteFieldMatricesOps.CharacteristicPolynomial( <D>, <A> ) . .  char pol
  1358. ##
  1359. FiniteFieldMatricesOps.CharacteristicPolynomial := function( D, A )
  1360.     local   F,  L,  d,  p,  P,  M,  N,  one,  R,  h,  v,  w,  i;
  1361.  
  1362.     # try to find the field of <A>
  1363.     if not IsBound(D.baseRing)  then
  1364.     F := Field( List( A, x -> Field(x).root ) );
  1365.     else
  1366.     F := D.baseRing;
  1367.     fi;
  1368.  
  1369.     # set up the other variables
  1370.     L := [];
  1371.     d := Length(A);
  1372.     P := [ Polynomial( F, [F.one] ) ];
  1373.     N := F.zero * A[1];
  1374.     IsVector(N);
  1375.  
  1376.     # span the whole vector space
  1377.     for i  in [ 1 .. d ]  do
  1378.  
  1379.     # did we already see the <i>.th base vector
  1380.         if not IsBound(L[i])  then
  1381.  
  1382.             # clear right sides <R>
  1383.             R := [];
  1384.  
  1385.             # start with the standard base vector
  1386.             v := ShallowCopy(N);
  1387.             v[i] := F.one;
  1388.  
  1389.             # <j>-1 gives the power of <A> we are looking at
  1390.             M := [];
  1391.  
  1392.             # spin vector around and construct polynomial
  1393.             repeat
  1394.  
  1395.                 # compute the head of <v>
  1396.                 h := DepthVector(v);
  1397.  
  1398.               # start with appropriate polynomial x^(<j>-1)
  1399.                 p := Copy(M);
  1400.                 Add(p, F.one);
  1401.                 Add(M, F.zero);
  1402.  
  1403.                 # divide by known left sides
  1404.                 w := Copy(v);
  1405.                 while h <= Length(w) and IsBound(L[h])  do
  1406.                     if IsBound(R[h])  then
  1407.                     AddCoeffs(p, R[h], -w[h]);
  1408.                     fi;
  1409.                     AddCoeffs(w, L[h], -w[h]);
  1410.                     h := DepthVector(w);
  1411.                 od;
  1412.  
  1413.                 # if <v> is not the zero vector try next power
  1414.                 if h <= Length(w)  then
  1415.                     R[h] := p * w[h]^-1;
  1416.                     L[h] := w * w[h]^-1;
  1417.                     v := v * A;
  1418.  
  1419.                 # otherwise, we know our polynomial
  1420.                 else
  1421.                     Add(P, Polynomial(F, p));
  1422.                 fi;
  1423.             until h > Length(w);
  1424.         fi;
  1425.     od;
  1426.  
  1427.     # compute the product of all polynomials
  1428.     return Product(P);
  1429.  
  1430. end;
  1431.  
  1432.  
  1433. #############################################################################
  1434. ##
  1435. #E  Emacs . . . . . . . . . . . . . . . . . . . . . . . local emacs variables
  1436. ##
  1437. ## Local Variables:
  1438. ## mode:           outline
  1439. ## outline-regexp: "#F\\|#V\\|#E"
  1440. ## eval:           (hide-body)
  1441. ## End:
  1442. ##
  1443.