axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / REP2.spad

)abbrev package REP2 RepresentationPackage2
++ Authors: Holger Gollan, Johannes Grabmeier
++ Date Created: 10 September 1987
++ Date Last Updated: 20 August 1990
++ Reference:
++   R. A. Parker: The Computer Calculation of Modular Characters
++   (The Meat-Axe), in M. D. Atkinson (Ed.), Computational Group Theory
++   Academic Press, Inc., London 1984
++   H. Gollan, J. Grabmeier: Algorithms in Representation Theory and
++    their Realization in the Computer Algebra System Scratchpad,
++    Bayreuther Mathematische Schriften, Heft 33, 1990, 1-23.
++ Description: 
++ RepresentationPackage2 provides functions for working with
++ modular representations of finite groups and algebra.
++ The routines in this package are created, using ideas of R. Parker,
++ (the meat-Axe) to get smaller representations from bigger ones, 
++ finding sub- and factormodules, or to show, that such the
++ representations are irreducible.
++ Note that most functions are randomized functions of Las Vegas type
++ every answer is correct, but with small probability
++ the algorithm fails to get an answer.

RepresentationPackage2(R) : SIG == CODE where 
  R : Ring

  OF    ==> OutputForm
  I     ==> Integer
  L     ==> List
  SM    ==> SquareMatrix
  M     ==> Matrix
  NNI   ==> NonNegativeInteger
  V     ==> Vector
  PI    ==> PositiveInteger
  B     ==> Boolean
  RADIX ==> RadixExpansion
 
  SIG ==> with
 
    completeEchelonBasis : V V R  ->  M R
      ++ completeEchelonBasis(lv) completes the basis lv assumed
      ++ to be in echelon form of a subspace of R**n (n the length
      ++ of all the vectors in lv with unit vectors to a basis of
      ++ R**n. It is assumed that the argument is not an empty
      ++ vector and that it is not the basis of the 0-subspace.
      ++ Note that the rows of the result correspond to the vectors 
      ++ of the basis.

    createRandomElement : (L M R, M R) -> M R
      ++ createRandomElement(aG,x) creates a random element of the group
      ++ algebra generated by aG. 

    if R has EuclideanDomain then   -- using rowEchelon

      cyclicSubmodule  : (L M R, V R) ->   V V R
        ++ cyclicSubmodule(lm,v) generates a basis as follows.
        ++ It is assumed that the size n of the vector equals the number
        ++ of rows and columns of the matrices. Then the matrices generate
        ++ a subalgebra, say \spad{A}, of the algebra of all square matrices of
        ++ dimension n. V R is an \spad{A}-module in the natural way.
        ++ cyclicSubmodule(lm,v) generates the R-Basis of Av as
        ++ described in section 6 of R. A. Parker's "The Meat-Axe".
        ++ Note that in contrast to the description in "The Meat-Axe" and to
        ++ standardBasisOfCyclicSubmodule the result is in echelon form.

      standardBasisOfCyclicSubmodule  : (L M R, V R) ->   M R
        ++ standardBasisOfCyclicSubmodule(lm,v) returns a matrix as follows.
        ++ It is assumed that the size n of the vector equals the number
        ++ of rows and columns of the matrices.  Then the matrices generate
        ++ a subalgebra, say \spad{A}, 
        ++ of the algebra of all square matrices of
        ++ dimension n. V R is an \spad{A}-module in the natural way.
        ++ standardBasisOfCyclicSubmodule(lm,v) calculates a matrix whose
        ++ non-zero column vectors are the R-Basis of Av achieved
        ++ in the way as described in section 6 
        ++ of R. A. Parker's "The Meat-Axe".
        ++ Note that in contrast to cyclicSubmodule, the result is not
        ++ in echelon form.

    if R has Field then  -- only because of inverse in SM

      areEquivalent? : (L M R, L M R, B, I) -> M R
        ++ areEquivalent?(aG0,aG1,randomelements,numberOfTries) tests
        ++ whether the two lists of matrices, all assumed of same
        ++ square shape, can be simultaneously conjugated by a non-singular
        ++ matrix. If these matrices represent the same group generators,
        ++ the representations are equivalent. 
        ++ The algorithm tries
        ++ numberOfTries times to create elements in the
        ++ generated algebras in the same fashion. If their ranks differ,
        ++ they are not equivalent. If an
        ++ isomorphism is assumed, then 
        ++ the kernel of an element of the first algebra
        ++ is mapped to the kernel of the corresponding element in the
        ++ second algebra. Now consider the one-dimensional ones.
        ++ If they generate the whole space (for example, irreducibility !)
        ++ we use standardBasisOfCyclicSubmodule to create the
        ++ only possible transition matrix. The method checks whether the 
        ++ matrix conjugates all corresponding matrices from aGi.
        ++ The way to choose the singular matrices is as in meatAxe.
        ++ If the two representations are equivalent, this routine
        ++ returns the transformation matrix TM with
        ++ aG0.i * TM = TM * aG1.i for all i. If the representations
        ++ are not equivalent, a small 0-matrix is returned.
        ++ Note that the case
        ++ with different sets of group generators cannot be handled.

      areEquivalent? : (L M R, L M R) -> M R
        ++ areEquivalent?(aG0,aG1) calls areEquivalent?(aG0,aG1,true,25).
        ++ Note that the choice of 25 was rather arbitrary.

      areEquivalent? : (L M R, L M R, I) -> M R
        ++ areEquivalent?(aG0,aG1,numberOfTries) calls
        ++ areEquivalent?(aG0,aG1,true,25).
        ++ Note that the choice of 25 was rather arbitrary.

      isAbsolutelyIrreducible? : (L M R, I) -> B
        ++ isAbsolutelyIrreducible?(aG, numberOfTries) uses
        ++ Norton's irreducibility test to check for absolute
        ++ irreduciblity, assuming if a one-dimensional kernel is found.
        ++ As no field extension changes create "new" elements
        ++ in a one-dimensional space, the criterium stays true
        ++ for every extension. The method looks for one-dimensionals only
        ++ by creating random elements (no fingerprints) since
        ++ a run of meatAxe would have proved absolute irreducibility
        ++ anyway.

      isAbsolutelyIrreducible? : L M R -> B
        ++ isAbsolutelyIrreducible?(aG) calls
        ++ isAbsolutelyIrreducible?(aG,25).
        ++ Note that the choice of 25 was rather arbitrary.

      split : (L M R, V R)  ->  L L M R
        ++ split(aG, vector) returns a subalgebra \spad{A} of all 
        ++ square matrix of dimension n as a list of list of matrices,
        ++ generated by the list of matrices aG, where n denotes both
        ++ the size of vector as well as the dimension of each of the
        ++ square matrices.
        ++ V R is an A-module in the natural way.
        ++ split(aG, vector) then checks whether the cyclic submodule
        ++ generated by vector is a proper submodule of V R.
        ++ If successful, it returns a two-element list, which contains
        ++ first the list of the representations of the submodule,
        ++ then the list of the representations of the factor module.
        ++ If the vector generates the whole module, a one-element list
        ++ of the old representation is given.
        ++ Note that a later version this should call the other split.

      split:  (L M R, V V R) -> L L M R
        ++ split(aG,submodule) uses a proper submodule of R**n
        ++ to create the representations of the submodule and of
        ++ the factor module.

    if (R has Finite) and (R has Field) then

      meatAxe  : (L M R, B, I, I)  ->  L L M R
        ++ meatAxe(aG,randomElements,numberOfTries, maxTests) returns
        ++ a 2-list of representations as follows.
        ++ All matrices of argument aG are assumed to be square
        ++ and of equal size.
        ++ Then \spad{aG} generates a subalgebra, say \spad{A}, of the algebra
        ++ of all square matrices of dimension n. V R is an A-module
        ++ in the usual way.  
        ++ meatAxe(aG,numberOfTries, maxTests) creates at most
        ++ numberOfTries random elements of the algebra, tests
        ++ them for singularity. If singular, it tries at most maxTests
        ++ elements of its kernel to generate a proper submodule.
        ++ If successful, a 2-list is returned: first, a list 
        ++ containing first the list of the
        ++ representations of the submodule, then a list of the
        ++ representations of the factor module.
        ++ Otherwise, if we know that all the kernel is already
        ++ scanned, Norton's irreducibility test can be used either
        ++ to prove irreducibility or to find the splitting.
        ++ If randomElements is false, the first 6 tries
        ++ use Parker's fingerprints.

      meatAxe : L M R -> L L M R
        ++ meatAxe(aG) calls meatAxe(aG,false,25,7) returns
        ++ a 2-list of representations as follows.
        ++ All matrices of argument \spad{aG} are assumed to be square
        ++ and of
        ++ equal size. Then \spad{aG} generates a subalgebra, 
        ++ say \spad{A}, of the algebra
        ++ of all square matrices of dimension n. V R is an A-module
        ++ in the usual way.  
        ++ meatAxe(aG) creates at most 25 random elements 
        ++ of the algebra, tests
        ++ them for singularity. If singular, it tries at most 7
        ++ elements of its kernel to generate a proper submodule.
        ++ If successful a list which contains first the list of the
        ++ representations of the submodule, then a list of the
        ++ representations of the factor module is returned.
        ++ Otherwise, if we know that all the kernel is already
        ++ scanned, Norton's irreducibility test can be used either
        ++ to prove irreducibility or to find the splitting.
        ++ Notes: the first 6 tries use Parker's fingerprints.
        ++ Also, 7 covers the case of three-dimensional kernels over 
        ++ the field with 2 elements.

      meatAxe: (L M R, B) ->  L L M R
        ++ meatAxe(aG, randomElements) calls meatAxe(aG,false,6,7),
        ++ only using Parker's fingerprints, if randomElemnts is false.
        ++ If it is true, it calls meatAxe(aG,true,25,7),
        ++ only using random elements.
        ++ Note that the choice of 25 was rather arbitrary.
        ++ Also, 7 covers the case of three-dimensional kernels over the field
        ++ with 2 elements.

      meatAxe : (L M R, PI)  ->  L L M R
        ++ meatAxe(aG, numberOfTries) calls 
        ++ meatAxe(aG,true,numberOfTries,7).
        ++ Notes: 7 covers the case of three-dimensional 
        ++ kernels over the field with 2 elements.

      scanOneDimSubspaces: (L V R, I)  ->  V R
        ++ scanOneDimSubspaces(basis,n) gives a canonical representative
        ++ of the n-th one-dimensional subspace of the vector space
        ++ generated by the elements of basis, all from R**n.
        ++ The coefficients of the representative are of shape
        ++ (0,...,0,1,*,...,*), * in R. If the size of R
        ++ is q, then there are (q**n-1)/(q-1) of them.
        ++ We first reduce n modulo this number, then find the
        ++ largest i such that +/[q**i for i in 0..i-1] <= n.
        ++ Subtracting this sum of powers from n results in an
        ++ i-digit number to basis q. This fills the positions of the
        ++ stars.
        -- would prefer to have (V V R,....   but nullSpace  results
        -- in L V R
 
  CODE ==> add
 
    -- import of domain and packages
    import OutputForm
 
    -- declarations and definitions of local variables and
    -- local function
 
    -- blockMultiply(a,b,li,n) assumes that a has n columns
    -- and b has n rows, li is a sublist of the rows of a and
    -- a sublist of the columns of b. The result is the
    -- multiplication of the (li x n) part of a with the
    -- (n x li) part of b. We need this, because just matrix
    -- multiplying the parts would require extra storage.
    blockMultiply: (M R, M R, L I, I) -> M R
    blockMultiply(a, b, li, n) ==
      matrix([[ +/[a(i,s) * b(s,j) for s in 1..n ] _
        for j in li  ]  for i in li])
 
    -- is local, because one should know all the results for smaller i
    fingerPrint: (NNI, M R, M R, M R) -> M R
    fingerPrint (i : NNI, a : M R, b : M R, x :M R) ==
      -- i > 2 only gives the correct result if the value of x from
      -- the parameter list equals the result of fingerprint(i-1,...)
      (i::PI) = 1 => x := a + b + a*b
      (i::PI) = 2 => x := (x + a*b)*b
      (i::PI) = 3 => x := a + b*x
      (i::PI) = 4 => x := x + b
      (i::PI) = 5 => x := x + a*b
      (i::PI) = 6 => x := x - a + b*a
      error "Sorry, but there are only 6 fingerprints!"
      x
 
    completeEchelonBasis(basis) ==
      dimensionOfSubmodule : NNI := #basis
      n : NNI := # basis.1
      indexOfVectorToBeScanned : NNI := 1
      row : NNI := dimensionOfSubmodule
      completedBasis : M R := zero(n, n)
      for i in 1..dimensionOfSubmodule repeat
        completedBasis := setRow_!(completedBasis, i, basis.i)
      if #basis <= n then
        newStart : NNI := 1
        for j in 1..n
          while indexOfVectorToBeScanned <= dimensionOfSubmodule repeat
            if basis.indexOfVectorToBeScanned.j = 0 then
              completedBasis(1+row,j) := 1  --put unit vector into basis
              row := row + 1
            else
              indexOfVectorToBeScanned := indexOfVectorToBeScanned + 1
            newStart : NNI := j + 1
        for j in newStart..n repeat
          completedBasis(j,j) := 1  --put unit vector into basis
      completedBasis
 
    createRandomElement(aG,algElt) ==
      numberOfGenerators : NNI := #aG
      -- randomIndex := randnum numberOfGenerators
      randomIndex := 1+(random()$Integer rem numberOfGenerators)
      algElt := algElt * aG.randomIndex
      -- randomIndxElement := randnum numberOfGenerators
      randomIndex := 1+(random()$Integer rem numberOfGenerators)
      algElt + aG.randomIndex
 
    if R has EuclideanDomain then

      cyclicSubmodule (lm : L M R, v : V R)  ==
        basis : M R := rowEchelon matrix list entries v
        -- normalizing the vector
        -- all these elements lie in the submodule generated by v
        furtherElts : L V R := [(lm.i*v)::V R for i in 1..maxIndex lm]
        --furtherElts has elements of the generated submodule. It will
        --will be checked whether they are in the span of the vectors
        --computed so far. Of course we stop if we have got the whole
        --space.
        while (^null furtherElts) and (nrows basis < #v)  repeat
          w : V R := first furtherElts
          nextVector : M R := matrix list entries w -- normalizing the vector
          -- will the rank change if we add this nextVector
          -- to the basis so far computed?
          addedToBasis : M R := vertConcat(basis, nextVector)
          if rank addedToBasis ^= nrows basis then
             basis := rowEchelon addedToBasis  -- add vector w to basis
             updateFurtherElts : L V R := _
               [(lm.i*w)::V R for i in 1..maxIndex lm]
             furtherElts := append (rest furtherElts, updateFurtherElts)
          else
             -- the vector w lies in the span of matrix, no updating
             -- of the basis
             furtherElts := rest furtherElts
        vector [row(basis, i) for i in 1..maxRowIndex basis]
 
      standardBasisOfCyclicSubmodule (lm : L M R, v : V R)  ==
        dim   : NNI := #v
        standardBasis : L L R := list(entries v)
        basis : M R := rowEchelon matrix list entries v
        -- normalizing the vector
        -- all these elements lie in the submodule generated by v
        furtherElts : L V R := [(lm.i*v)::V R for i in 1..maxIndex lm]
        --furtherElts has elements of the generated submodule. It will
        --will be checked whether they are in the span of the vectors
        --computed so far. Of course we stop if we  have got the whole
        --space.
        while (^null furtherElts) and (nrows basis < #v)  repeat
          w : V R := first furtherElts
          nextVector : M R := matrix list entries w  -- normalizing the vector
          -- will the rank change if we add this nextVector
          -- to the basis so far computed?
          addedToBasis : M R := vertConcat(basis, nextVector)
          if rank addedToBasis ^= nrows basis then
             standardBasis := cons(entries w, standardBasis)
             basis := rowEchelon addedToBasis  -- add vector w to basis
             updateFurtherElts : L V R := _
               [lm.i*w for i in 1..maxIndex lm]
             furtherElts := append (rest furtherElts, updateFurtherElts)
          else
             -- the vector w lies in the span of matrix, therefore
             -- no updating of matrix
             furtherElts := rest furtherElts
        transpose matrix standardBasis
 
    if R has Field then  -- only because of inverse in Matrix
 
      -- as conditional local functions, *internal have to be here
 
      splitInternal: (L M R, V R, B) -> L L M R
      splitInternal(algebraGenerators : L M R, vector: V R,doSplitting? : B) ==
        n : I := # vector    -- R-rank of representation module =
                               -- degree of representation
        submodule : V V R := cyclicSubmodule (algebraGenerators,vector)
        rankOfSubmodule : I := # submodule  -- R-Rank of submodule
        submoduleRepresentation    : L M R := nil()
        factormoduleRepresentation : L M R := nil()
        if n ^= rankOfSubmodule then
          messagePrint "  A proper cyclic submodule is found."
          if doSplitting? then   -- no else !!
            submoduleIndices : L I := [i for i in 1..rankOfSubmodule]
            factormoduleIndices : L I := [i for i in (1+rankOfSubmodule)..n]
            transitionMatrix : M R := _
              transpose completeEchelonBasis submodule
            messagePrint "  Transition matrix computed"
            inverseTransitionMatrix : M R :=  _
             autoCoerce(inverse transitionMatrix)$Union(M R,"failed")
            messagePrint "  The inverse of the transition matrix computed"
            messagePrint "  Now transform the matrices"
            for i in 1..maxIndex algebraGenerators repeat
              helpMatrix : M R := inverseTransitionMatrix * algebraGenerators.i
              -- in order to not create extra space and regarding the fact
              -- that we only want the two blocks in the main diagonal we
              -- multiply with the aid of the local function blockMultiply
              submoduleRepresentation := cons( blockMultiply( _
                helpMatrix,transitionMatrix,submoduleIndices,n), _
                submoduleRepresentation)
              factormoduleRepresentation := cons( blockMultiply( _
                helpMatrix,transitionMatrix,factormoduleIndices,n), _
                factormoduleRepresentation)
          [reverse submoduleRepresentation, reverse _
            factormoduleRepresentation]
        else -- represesentation is irreducible
          messagePrint "  The generated cyclic submodule was not proper"
          [algebraGenerators]
 
      irreducibilityTestInternal: (L M R, M R, B) -> L L M R
      irreducibilityTestInternal(algebraGenerators,_
          singularMatrix,split?) ==
        algebraGeneratorsTranspose : L M R := [transpose _
           algebraGenerators.j for j in 1..maxIndex algebraGenerators]
        xt : M R := transpose singularMatrix
        messagePrint _
         "  We know that all the cyclic submodules generated by all"
        messagePrint _
         "    non-trivial element of the singular matrix under view are"
        messagePrint _
         "    not proper, hence Norton's irreducibility test can be done:"
        -- actually we only would need one (!) non-trivial element from
        -- the kernel of xt, such an element must exist as the transpose
        -- of a singular matrix is of course singular. Question: Can
        -- we get it more easily from the kernel of x = singularMatrix?
        kernel : L V R := nullSpace xt
        result : L L M R :=  _
          splitInternal(algebraGeneratorsTranspose,first kernel,split?)
        if null rest result then  -- this means first kernel generates
          -- the whole module
          if 1 = #kernel then
             messagePrint "  Representation is absolutely irreducible"
          else
             messagePrint "  Representation is irreducible, but we don't know "
             messagePrint "    whether it is absolutely irreducible"
        else
          if split? then
            messagePrint _
             "  Representation is not irreducible and it will be split:"
            -- these are the dual representations, so calculate the
            -- dual to get the desired result, "transpose inverse"
            -- improvements??
            for i in 1..maxIndex result repeat
              for j in 1..maxIndex (result.i) repeat
                mat : M R := result.i.j
                result.i.j := _
                  transpose autoCoerce(inverse mat)$Union(M R,"failed")
          else
            messagePrint _
             "  Representation is not irreducible, use meatAxe to split"
        -- if "split?" then dual representation interchange factor
        -- and submodules, hence reverse
        reverse result
 
      -- exported functions for FiniteField-s.
 
      areEquivalent? (aG0, aG1) ==
        areEquivalent? (aG0, aG1, true, 25)
 
      areEquivalent? (aG0, aG1, numberOfTries) ==
        areEquivalent? (aG0, aG1, true, numberOfTries)
 
      areEquivalent? (aG0, aG1, randomelements, numberOfTries) ==
          result : B := false
          transitionM : M R := zero(1, 1)
          numberOfGenerators  : NNI := #aG0
          -- need a start value for creating random matrices:
          -- if we switch to randomelements later, we take the last
          -- fingerprint.
          if randomelements then   -- random should not be from I
             --randomIndex  : I   := randnum numberOfGenerators
             randomIndex := 1+(random()$Integer rem numberOfGenerators)
             x0 : M R := aG0.randomIndex
             x1 : M R := aG1.randomIndex
          n : NNI := #row(x0,1)   -- degree  of representation
          foundResult : B := false
          for i in 1..numberOfTries until foundResult repeat
            -- try to create a non-singular element of the algebra
            -- generated by "aG". If only two generators,
            -- i < 7 and not "randomelements" use Parker's  fingerprints
            -- i >= 7 create random elements recursively:
            -- x_i+1 :=x_i * mr1 + mr2, where mr1 and mr2 are randomly
            -- chosen elements form "aG".
            if i = 7 then randomelements := true
            if randomelements then
               --randomIndex := randnum numberOfGenerators
               randomIndex := 1+(random()$Integer rem numberOfGenerators)
               x0 := x0 * aG0.randomIndex
               x1 := x1 * aG1.randomIndex
               --randomIndex := randnum numberOfGenerators
               randomIndex := 1+(random()$Integer rem numberOfGenerators)
               x0 := x0 + aG0.randomIndex
               x1 := x1 + aG1.randomIndex
            else
               x0 := fingerPrint (i, aG0.0, aG0.1 ,x0)
               x1 := fingerPrint (i, aG1.0, aG1.1 ,x1)
            -- test singularity of x0 and x1
            rk0 : NNI := rank x0
            rk1 : NNI := rank x1
            rk0 ^= rk1 =>
              messagePrint  "Dimensions of kernels differ"
              foundResult := true
              result := false
            -- can assume dimensions are equal
            rk0 ^= n - 1 =>
              -- not of any use here if kernel not one-dimensional
              if randomelements then
                messagePrint  "Random element in generated algebra does"
                messagePrint  "  not have a one-dimensional kernel"
              else
                messagePrint  "Fingerprint element in generated algebra does"
                messagePrint  "  not have a one-dimensional kernel"
            -- can assume dimensions are equal and equal to n-1
            if randomelements then
              messagePrint  "Random element in generated algebra has"
              messagePrint  "  one-dimensional kernel"
            else
              messagePrint  "Fingerprint element in generated algebra has"
              messagePrint  "  one-dimensional kernel"
            kernel0 : L V R := nullSpace x0
            kernel1 : L V R := nullSpace x1
            baseChange0 : M R := standardBasisOfCyclicSubmodule(_
              aG0,kernel0.1)
            baseChange1 : M R := standardBasisOfCyclicSubmodule(_
              aG1,kernel1.1)
            (ncols baseChange0) ^= (ncols baseChange1) =>
              messagePrint "  Dimensions of generated cyclic submodules differ"
              foundResult := true
              result := false
            -- can assume that dimensions of cyclic submodules are equal
            (ncols baseChange0) = n =>   -- full dimension
              transitionM := baseChange0 * _
                autoCoerce(inverse baseChange1)$Union(M R,"failed")
              foundResult := true
              result := true
              for j in 1..numberOfGenerators while result repeat
                if (aG0.j*transitionM) ^= (transitionM*aG1.j) then
                  result := false
                  transitionM := zero(1 ,1)
                  messagePrint _
                   "  There is no isomorphism, as the only possible one"
                  messagePrint "    fails to do the necessary base change"
            -- can assume that dimensions of cyclic submodules are not "n"
            messagePrint _
              "  Generated cyclic submodules have equal, but not full"
            messagePrint  "    dimension, hence we can not draw any conclusion"
          -- here ends the for-loop
          if not foundResult then
            messagePrint  " "
            messagePrint  "Can neither prove equivalence nor inequivalence."
            messagePrint  "  Try again."
          else
            if result then
              messagePrint  " "
              messagePrint  "Representations are equivalent."
            else
              messagePrint  " "
              messagePrint  "Representations are not equivalent."
          transitionM
 
      isAbsolutelyIrreducible?(aG) == isAbsolutelyIrreducible?(aG,25)
 
      isAbsolutelyIrreducible?(aG, numberOfTries) ==
        result : B := false
        numberOfGenerators  : NNI := #aG
        -- need a start value for creating random matrices:
        -- randomIndex  : I   := randnum numberOfGenerators
        randomIndex := 1+(random()$Integer rem numberOfGenerators)
        x : M R := aG.randomIndex
        n : NNI := #row(x,1)   -- degree  of representation
        foundResult : B := false
        for i in 1..numberOfTries until foundResult repeat
          -- try to create a non-singular element of the algebra
          -- generated by "aG", dimension of its kernel being 1.
          -- create random elements recursively:
          -- x_i+1 :=x_i * mr1 + mr2, where mr1 and mr2 are randomly
          -- chosen elements form "aG".
          -- randomIndex := randnum numberOfGenerators
          randomIndex := 1+(random()$Integer rem numberOfGenerators)
          x := x * aG.randomIndex
          --randomIndex := randnum numberOfGenerators
          randomIndex := 1+(random()$Integer rem numberOfGenerators)
          x := x + aG.randomIndex
          -- test whether rank of x is n-1
          rk : NNI := rank x
          if rk = n - 1 then
            foundResult := true
            messagePrint "Random element in generated algebra has"
            messagePrint "  one-dimensional kernel"
            kernel : L V R := nullSpace x
            if n=#cyclicSubmodule(aG, first kernel) then
              result := (irreducibilityTestInternal(aG,x,false)).1 ^= nil()_
                 $(L M R)
            else -- we found a proper submodule
              result := false
              --split(aG,kernel.1) -- to get the splitting
          else -- not of any use here if kernel not one-dimensional
            messagePrint "Random element in generated algebra does"
            messagePrint "  not have a one-dimensional kernel"
        -- here ends the for-loop
        if not foundResult then
          messagePrint "We have not found a one-dimensional kernel so far,"
          messagePrint "  as we do a random search you could try again"
        --else
        --  if not result then
        --    messagePrint "Representation is not irreducible."
        --  else
        --    messagePrint "Representation is irreducible."
        result
 
      split(algebraGenerators: L M R, vector: V R) ==
        splitInternal(algebraGenerators, vector, true)
 
      split(algebraGenerators : L M R, submodule: V V R)== --not zero submodule
        n : NNI := #submodule.1 -- R-rank of representation module =
                                -- degree of representation
        rankOfSubmodule : I := (#submodule) :: I --R-Rank of submodule
        submoduleRepresentation    : L M R := nil()
        factormoduleRepresentation : L M R := nil()
        submoduleIndices : L I := [i for i in 1..rankOfSubmodule]
        factormoduleIndices : L I := [i for i in (1+rankOfSubmodule)..(n::I)]
        transitionMatrix : M R := _
          transpose completeEchelonBasis submodule
        messagePrint "  Transition matrix computed"
        inverseTransitionMatrix : M R :=
          autoCoerce(inverse transitionMatrix)$Union(M R,"failed")
        messagePrint "  The inverse of the transition matrix computed"
        messagePrint "  Now transform the matrices"
        for i in 1..maxIndex algebraGenerators repeat
          helpMatrix : M R := inverseTransitionMatrix * algebraGenerators.i
          -- in order to not create extra space and regarding the fact
          -- that we only want the two blocks in the main diagonal we
          -- multiply with the aid of the local function blockMultiply
          submoduleRepresentation := cons( blockMultiply( _
            helpMatrix,transitionMatrix,submoduleIndices,n), _
            submoduleRepresentation)
          factormoduleRepresentation := cons( blockMultiply( _
            helpMatrix,transitionMatrix,factormoduleIndices,n), _
            factormoduleRepresentation)
        cons(reverse submoduleRepresentation, list( reverse _
          factormoduleRepresentation)::(L L M R))
 
 
    -- the following is "under"  "if R has Field", as there are compiler
    -- problems with conditinally defined local functions, it
    -- doesn't know, that "FiniteField" has "Field".
 
 
      -- we are scanning through the vectorspaces
      if (R has Finite) and (R has Field) then
 
        meatAxe(algebraGenerators, randomelements, numberOfTries, _
           maxTests) ==
          numberOfGenerators  : NNI := #algebraGenerators
          result : L L M R := nil()$(L L M R)
          q   : PI  := size()$R:PI
          -- need a start value for creating random matrices:
          -- if we switch to randomelements later, we take the last
          -- fingerprint.
          if randomelements then   -- random should not be from I
             --randomIndex  : I   := randnum numberOfGenerators
             randomIndex := 1+(random()$Integer rem numberOfGenerators)
             x : M R := algebraGenerators.randomIndex
          foundResult : B := false
          for i in 1..numberOfTries until foundResult repeat
            -- try to create a non-singular element of the algebra
            -- generated by "algebraGenerators". If only two generators,
            -- i < 7 and not "randomelements" use Parker's  fingerprints
            -- i >= 7 create random elements recursively:
            -- x_i+1 :=x_i * mr1 + mr2, where mr1 and mr2 are randomly
            -- chosen elements form "algebraGenerators".
            if i = 7 then randomelements := true
            if randomelements then
               --randomIndex := randnum numberOfGenerators
               randomIndex := 1+(random()$Integer rem numberOfGenerators)
               x := x * algebraGenerators.randomIndex
               --randomIndex := randnum numberOfGenerators
               randomIndex := 1+(random()$Integer rem numberOfGenerators)
               x := x + algebraGenerators.randomIndex
            else
               x := fingerPrint (i, algebraGenerators.1,_
                 algebraGenerators.2 , x)
            -- test singularity of x
            n : NNI := #row(x, 1)  -- degree  of representation
            if (rank x) ^= n then  -- x singular
              if randomelements then
                 messagePrint "Random element in generated algebra is singular"
              else
                 messagePrint _
                  "Fingerprint element in generated algebra is singular"
              kernel : L V R := nullSpace x
              -- the first number is the maximal number of one dimensional
              -- subspaces of the kernel, the second is a user given
              -- constant
              numberOfOneDimSubspacesInKernel : I := (q**(#kernel)-1)quo(q-1)
              numberOfTests : I :=  _
                min(numberOfOneDimSubspacesInKernel, maxTests)
              for j in 1..numberOfTests repeat
                 --we create an element in the kernel, there is a good
                 --probability for it to generate a proper submodule, the
                 --called "split" does the further work:
                 result := _
                   split(algebraGenerators,scanOneDimSubspaces(kernel,j))
                 -- we had "not null rest result" directly in the following
                 -- if .. then, but the statment there foundResult := true
                 -- didn't work properly
                 foundResult :=  not null rest result
                 if foundResult then
                   leave -- inner for-loop
                   -- finish here with result
                 else  -- no proper submodule
                   -- we were not successfull, i.e gen. submodule was
                   -- not proper, if the whole kernel is already scanned,
                   -- Norton's irreducibility test is used now.
                   if (j+1)>numberOfOneDimSubspacesInKernel then
                     -- we know that all the cyclic submodules generated
                     -- by all non-trivial elements of the kernel are proper.
                     foundResult := true
                     result : L L M R := irreducibilityTestInternal (_
                       algebraGenerators,x,true)
                     leave  -- inner for-loop
              -- here ends the inner for-loop
            else  -- x non-singular
              if randomelements then
                messagePrint _
                 "Random element in generated algebra is non-singular"
              else
                messagePrint _
                 "Fingerprint element in generated algebra is non-singular"
          -- here ends the outer for-loop
          if not foundResult then
             result : L L M R := [nil()$(L M R), nil()$(L M R)]
             messagePrint " "
             messagePrint "Sorry, no result, try meatAxe(...,true)"
             messagePrint "  or consider using an extension field."
          result
 
        meatAxe (algebraGenerators) ==
          meatAxe(algebraGenerators, false, 25, 7)
 
        meatAxe (algebraGenerators, randomElements?) ==
          randomElements? => meatAxe (algebraGenerators, true, 25, 7)
          meatAxe(algebraGenerators, false, 6, 7)
 
        meatAxe (algebraGenerators:L M R, numberOfTries:PI) ==
          meatAxe (algebraGenerators, true, numberOfTries, 7)
 
        scanOneDimSubspaces(basis,n) ==
          -- "dimension" of subspace generated by "basis"
          dim : NNI := #basis
          -- "dimension of the whole space:
          nn : NNI := #(basis.1)
          q : NNI := size()$R
          -- number of all one-dimensional subspaces:
          nred : I := n rem ((q**dim -1) quo (q-1))
          pos : I := nred
          i : I := 0
          for i in 0..dim-1 while nred >= 0 repeat
            pos := nred
            nred := nred - (q**i)
          i  := if i = 0 then 0 else i-1
          coefficients : V R := new(dim,0$R)
          coefficients.(dim-i) := 1$R
          iR : L I := wholeRagits(pos::RADIX q)
          for j in 1..(maxIndex iR) repeat
            coefficients.(dim-((#iR)::I) +j) := index((iR.j+(q::I))::PI)$R
          result : V R := new(nn,0)
          for i in 1..maxIndex coefficients repeat
            newAdd : V R := coefficients.i * basis.i
            for j in 1..nn repeat
              result.j := result.j + newAdd.j
          result