gap-libs 4r6p5-3 / usr / share / gap / lib / algebra.gi

#############################################################################
##
#W  algebra.gi                  GAP library                     Thomas Breuer
#W                                                        and Willem de Graaf
##
##
#Y  Copyright (C)  1997,  Lehrstuhl D für Mathematik,  RWTH Aachen,  Germany
#Y  (C) 1998 School Math and Comp. Sci., University of St Andrews, Scotland
#Y  Copyright (C) 2002 The GAP Group
##
##  This file contains generic methods for algebras and algebras-with-one.
##


#############################################################################
##
#M  Representative( <A> ) . . . . . . . . one element of a left operator ring
##
InstallMethod( Representative,
    "for left operator ring with known generators",
    [ IsLeftOperatorRing and HasGeneratorsOfLeftOperatorRing ],
    RepresentativeFromGenerators( GeneratorsOfLeftOperatorRing ) );


#############################################################################
##
#M  Representative( <A> ) . . .  one element of a left operator ring-with-one
##
InstallMethod( Representative,
    "for left operator ring-with-one with known generators",
    [ IsLeftOperatorRingWithOne and HasGeneratorsOfLeftOperatorRingWithOne ],
    RepresentativeFromGenerators( GeneratorsOfLeftOperatorRingWithOne ) );


#############################################################################
##
#M  FLMLORByGenerators( <R>, <gens> ) . . . .  <R>-FLMLOR generated by <gens>
#M  FLMLORByGenerators( <R>, <gens>, <zero> )
##
InstallMethod( FLMLORByGenerators,
    "for ring and collection",
    [ IsRing, IsCollection ],
    function( R, gens )
    local A;
    A:= Objectify( NewType( FamilyObj( gens ),
                            IsFLMLOR and IsAttributeStoringRep ),
                   rec() );
    SetLeftActingDomain( A, R );
    SetGeneratorsOfLeftOperatorRing( A, AsList( gens ) );

    CheckForHandlingByNiceBasis( R, gens, A, false );
    return A;
    end );

InstallOtherMethod( FLMLORByGenerators,
    "for ring, homogeneous list, and ring element",
    [ IsRing, IsHomogeneousList, IsRingElement ],
    function( R, gens, zero )
    local A;
    A:= Objectify( NewType( CollectionsFamily( FamilyObj( zero ) ),
                            IsFLMLOR and IsAttributeStoringRep ),
                   rec() );
    SetLeftActingDomain( A, R );
    SetGeneratorsOfLeftOperatorRing( A, gens );
    SetZero( A, zero );

    if IsEmpty( gens ) then
      SetIsTrivial( A, true );
      SetGeneratorsOfLeftModule( A, gens );
    fi;

    CheckForHandlingByNiceBasis( R, gens, A, zero );
    return A;
    end );


#############################################################################
##
#M  FLMLORWithOneByGenerators( <R>, <gens> )  unit. <R>-FLMLOR gen. by <gens>
#M  FLMLORWithOneByGenerators( <R>, <gens>, <zero> )
##
InstallMethod( FLMLORWithOneByGenerators,
    "for ring and collection",
    [ IsRing, IsCollection ],
    function( R, gens )
    local A;
    A:= Objectify( NewType( FamilyObj( gens ),
                            IsFLMLORWithOne and IsAttributeStoringRep ),
                   rec() );
    SetLeftActingDomain( A, R );
    SetGeneratorsOfLeftOperatorRingWithOne( A, AsList( gens ) );

    CheckForHandlingByNiceBasis( R, gens, A, false );
    return A;
    end );

InstallOtherMethod( FLMLORWithOneByGenerators,
    "for ring, homogeneous list, and ring element",
    [ IsRing, IsHomogeneousList, IsRingElement ],
    function( R, gens, zero )
    local A;
    A:= Objectify( NewType( CollectionsFamily( FamilyObj( zero ) ),
                            IsFLMLORWithOne and IsAttributeStoringRep ),
                   rec() );
    SetLeftActingDomain( A, R );
    SetGeneratorsOfLeftOperatorRingWithOne( A, AsList( gens ) );
    SetZero( A, zero );

    CheckForHandlingByNiceBasis( R, gens, A, zero );
    return A;
    end );


#############################################################################
##
#F  Algebra( <F>, <gens> )
#F  Algebra( <F>, <gens>, <zero> )
#F  Algebra( <F>, <gens>, "basis" )
#F  Algebra( <F>, <gens>, <zero>, "basis" )
##
InstallGlobalFunction( FLMLOR, function( arg )
    local A;

    # ring and list of generators
    if Length( arg ) = 2 and IsRing( arg[1] )
                         and IsList( arg[2] ) and 0 < Length( arg[2] ) then
      A:= FLMLORByGenerators( arg[1], arg[2] );

    # ring, list of generators plus zero
    elif Length( arg ) = 3 and IsRing( arg[1] )
                           and IsList( arg[2] ) then
      if arg[3] = "basis" then
        A:= FLMLORByGenerators( arg[1], arg[2] );
        UseBasis( A, arg[2] );
      else
        A:= FLMLORByGenerators( arg[1], arg[2], arg[3] );
      fi;

    # ring, list of generators plus zero
    elif Length( arg ) = 4 and IsRing( arg[1] )
                           and IsList( arg[2] )
                           and arg[4] = "basis" then
      A:= FLMLORByGenerators( arg[1], arg[2], arg[3] );
      UseBasis( A, arg[2] );

    # no argument given, error
    else
      Error( "usage: FLMLOR( <F>, <gens> ), ",
             "FLMLOR( <F>, <gens>, <zero> )" );
    fi;

    # Return the result.
    return A;
end );


#############################################################################
##
#F  Subalgebra( <A>, <gens> ) . . . . . subalgebra of <A> generated by <gens>
#F  Subalgebra( <A>, <gens>, "basis" )
##
InstallGlobalFunction( SubFLMLOR, function( arg )
    local S;
    if    Length( arg ) <= 1
       or not IsFLMLOR( arg[1] )
       or not IsHomogeneousList( arg[2] ) then
      Error( "first argument must be a FLMLOR,\n",
             "second argument must be a list of generators" );

    elif IsEmpty( arg[2] ) then

      return SubFLMLORNC( arg[1], arg[2] );

    elif     IsIdenticalObj( FamilyObj( arg[1] ),
                          FamilyObj( arg[2] ) )
         and ForAll( arg[2], v -> v in arg[1] ) then

      S:= FLMLORByGenerators( LeftActingDomain( arg[1] ), arg[2] );
      SetParent( S, arg[1] );
      if Length( arg ) = 3 and arg[3] = "basis" then
        UseBasis( S, arg[2] );
      fi;
      return S;

    fi;
    Error( "usage: SubFLMLOR( <V>, <gens> [, \"basis\"] )" );
end );


#############################################################################
##
#F  SubalgebraNC( <A>, <gens>, "basis" )
#F  SubalgebraNC( <A>, <gens> )
##
InstallGlobalFunction( SubFLMLORNC, function( arg )
    local S;
    if IsEmpty( arg[2] ) then
      S:= Objectify( NewType( FamilyObj( arg[1] ),
                                  IsFLMLOR
                              and IsTrivial
                              and IsTwoSidedIdealInParent
                              and IsAttributeStoringRep ),
                     rec() );
      SetLeftActingDomain( S, LeftActingDomain( arg[1] ) );
      SetGeneratorsOfLeftModule( S, AsList( arg[2] ) );
    else
      S:= FLMLORByGenerators( LeftActingDomain( arg[1] ), arg[2] );
    fi;
    if Length( arg ) = 3 and arg[3] = "basis" then
      UseBasis( S, arg[2] );
    fi;
    SetParent( S, arg[1] );
    return S;
end );


#############################################################################
##
#F  AlgebraWithOne( <F>, <gens> )
#F  AlgebraWithOne( <F>, <gens>, <zero> )
#F  AlgebraWithOne( <F>, <gens>, "basis" )
#F  AlgebraWithOne( <F>, <gens>, <zero>, "basis" )
##
InstallGlobalFunction( FLMLORWithOne, function( arg )
    local A;

    # ring and list of generators
    if Length( arg ) = 2 and IsRing( arg[1] )
                         and IsList( arg[2] ) and 0 < Length( arg[2] ) then
      A:= FLMLORWithOneByGenerators( arg[1], arg[2] );

    # ring, list of generators plus zero
    elif Length( arg ) = 3 and IsRing( arg[1] )
                           and IsList( arg[2] ) then
      if arg[3] = "basis" then
        A:= FLMLORWithOneByGenerators( arg[1], arg[2] );
        UseBasis( A, arg[2] );
      else
        A:= FLMLORWithOneByGenerators( arg[1], arg[2], arg[3] );
      fi;

    # ring, list of generators plus zero
    elif Length( arg ) = 4 and IsRing( arg[1] )
                           and IsList( arg[2] )
                           and arg[4] = "basis" then
      A:= FLMLORWithOneByGenerators( arg[1], arg[2], arg[3] );
      UseBasis( A, arg[2] );

    # no argument given, error
    else
      Error( "usage: FLMLORWithOne( <F>, <gens> ), ",
             "FLMLORWithOne( <F>, <gens>, <zero> )" );
    fi;

    # Return the result.
    return A;
end );


#############################################################################
##
#F  SubalgebraWithOne( <A>, <gens> )   subalg.-with-one of <A> gen. by <gens>
##
InstallGlobalFunction( SubFLMLORWithOne, function( arg )
    local S;
    if Length( arg ) <= 1 or not IsFLMLOR( arg[1] )
                          or not IsHomogeneousList( arg[2] ) then

      Error( "first argument must be a FLMLOR,\n",
             "second argument must be a list of generators" );

    elif IsEmpty( arg[2] ) then

      return SubFLMLORWithOneNC( arg[2], arg[2] );

    elif     IsIdenticalObj( FamilyObj( arg[1] ),
                          FamilyObj( arg[2] ) )
         and ForAll( arg[2], v -> v in arg[1] )
         and ( IsFLMLORWithOne( arg[1] ) or One( arg[1] ) <> fail ) then
      S:= FLMLORWithOneByGenerators( LeftActingDomain( arg[1] ), arg[2] );
      SetParent( S, arg[1] );
      if Length( arg ) = 3 and arg[3] = "basis" then
        UseBasis( S, arg[2] );
      fi;
      return S;
    fi;
    Error( "usage: SubFLMLORWithOne( <V>, <gens> [, \"basis\"] )" );
end );


#############################################################################
##
#F  SubalgebraWithOneNC( <A>, <gens> )
##
InstallGlobalFunction( SubFLMLORWithOneNC, function( arg )
    local S, gens;
    if IsEmpty( arg[2] ) then

      # Note that `S' is in general not trivial,
      # and if we call `Objectify' here then `S' does not get a special
      # representation (e.g., as a matrix algebra).
      # So the argument that special methods would catch this case
      # does not hold!
      gens:= [ One( arg[1] ) ];
      S:= FLMLORWithOneByGenerators( LeftActingDomain( arg[1] ), gens );
      UseBasis( S, gens );

    else

      S:= FLMLORWithOneByGenerators( LeftActingDomain( arg[1] ), arg[2] );
      if Length( arg ) = 3 and arg[3] = "basis" then
        UseBasis( S, arg[2] );
      fi;

    fi;

    SetParent( S, arg[1] );
    return S;
end );


#############################################################################
##
#M  LieAlgebraByDomain( <A> )
##
##  The Lie algebra of the associative algebra <A>
##
InstallMethod( LieAlgebraByDomain,
    "for an algebra",
    [ IsAlgebra ],
    function( A )
       local T, n, zero, nullvec, S, i, j, k, m, cfs, cij, cji;

       if not IsAssociative( A ) then TryNextMethod(); fi;

# We construct a structure constants table for the  Lie algebra
# corresponding to <A>. If the structure constants of <A> are given by
# d_{ij}^k, then the structure constants of the Lie algebra will be given
# by d_{ij}^k - d_{ji}^k.


       T:= StructureConstantsTable( Basis( A ) );
       n:= Dimension( A );
       zero:= Zero( LeftActingDomain( A ) );
       nullvec:= List( [1..n], x -> zero );
       S:= EmptySCTable( n, zero, "antisymmetric" );
       for i in [1..n] do
         for j in [i+1..n] do
           cfs:= ShallowCopy( nullvec );
           cij:= T[i][j]; cji:= T[j][i];
           for m in [1..Length(cij[1])] do
             k:= cij[1][m];
             cfs[k]:= cfs[k] + cij[2][m];
           od;
           for m in [1..Length(cji[1])] do
             k:= cji[1][m];
             cfs[k]:= cfs[k] - cji[2][m];
           od;
           
           cij:= [ ];
           
           for m in [1..n] do
             if cfs[m] <> zero then
                 Add( cij, cfs[m] );
                 Add( cij, m );
             fi;
           od;
           SetEntrySCTable( S, i, j, cij );
         od;
       od;

       return LieAlgebraByStructureConstants( LeftActingDomain( A ), S );
  end );


#############################################################################
##
#F  LieAlgebra( <A> )
#F  LieAlgebra( <F>, <gens> )
#F  LieAlgebra( <F>, <gens>, <zero> )
#F  LieAlgebra( <F>, <gens>, "basis" )
#F  LieAlgebra( <F>, <gens>, <zero>, "basis" )
##
InstallGlobalFunction( LieAlgebra, function( arg )
#T check that the families have the same characteristic?
#T `CharacteristicFamily' ?
    local A,gens;

    # In the case of one domain argument,
    # construct the isomorphic Lie algebra.
    if Length( arg ) = 1 and IsDomain( arg[1] ) then
      A:= LieAlgebraByDomain( arg[1] );

    # division ring and list of generators
    elif Length( arg ) >= 2 and IsList( arg[2] ) then

       gens:= List( arg[2], x -> LieObject( x ) );

       if Length( arg ) = 2 and IsDivisionRing( arg[1] )
                            and 0 < Length( arg[2] ) then
         A:= AlgebraByGenerators( arg[1], gens );

       # division ring, list of generators plus zero
       elif Length( arg ) = 3 and IsDivisionRing( arg[1] ) then
          if arg[3] = "basis" then
            A:= AlgebraByGenerators( arg[1], gens );
            UseBasis( A, gens );
          else
            A:= AlgebraByGenerators( arg[1], gens, arg[3] );
          fi;

       # division ring, list of generators plus zero
       elif Length( arg ) = 4 and IsDivisionRing( arg[1] )
                              and arg[4] = "basis" then
          A:= AlgebraByGenerators( arg[1], gens, arg[3] );
          UseBasis( A, gens );

       else
         Error( "usage: LieAlgebra( <F>, <gens> ), ",
              "LieAlgebra( <F>, <gens>, <zero> ), LieAlgebra( <D> )");

       fi;
    # no argument given, error
    else
      Error( "usage: LieAlgebra( <F>, <gens> ), ",
             "LieAlgebra( <F>, <gens>, <zero> ), LieAlgebra( <D> )");
    fi;

    # Return the result.
    return A;
end );


#############################################################################
##
#F  EmptySCTable( <dim>, <zero> )
#F  EmptySCTable( <dim>, <zero>, \"symmetric\" )
#F  EmptySCTable( <dim>, <zero>, \"antisymmetric\" )
##
InstallGlobalFunction( EmptySCTable, function( arg )
    local dim, T, entry, i;

    if     2 <= Length( arg )
       and IsInt( arg[1] ) and IsZero( arg[2] )
       and ( Length( arg ) = 2 or IsString( arg[3] ) ) then

      dim:= arg[1];
      T:= [];
      entry:= Immutable( [ [], [] ] );
      for i in [ 1 .. dim ] do
        T[i]:= List( [ 1 .. dim ], x -> entry );
      od;

      # Store the symmetry flag.
      if Length( arg ) = 3 then
        if   arg[3] = "symmetric" then
          Add( T, 1 );
        elif arg[3] = "antisymmetric" then
          Add( T, -1 );
        else
          Error("third argument must be \"symmetric\" or \"antisymmetric\"");
        fi;
      else
        Add( T, 0 );
      fi;

      # Store the zero coefficient.
      Add( T, arg[2] );

    else
      Error( "usage: EmptySCTable( <dim>, <zero> [,\"symmetric\"] )" );
    fi;

    return T;
end );


#############################################################################
##
#F  SetEntrySCTable( <T>, <i>, <j>, <list> )
##
InstallGlobalFunction( SetEntrySCTable, function( T, i, j, list )
    local range, zero, Fam, entry, k, val, pos;

    # Check that `i' and `j' are admissible.
    range:= [ 1 .. Length( T ) - 2 ];
    if   not i in range then
      Error( "<i> must lie in ", range );
    elif not j in range then
      Error( "<j> must lie in ", range );
    fi;

    # Check `list', and construct the table entry.
    zero:= T[ Length( T ) ];
    Fam:= FamilyObj( zero );
    entry:= [ [], [] ];
    for k in [ 1, 3 .. Length( list ) -1 ] do

      val:= list[k];
      pos:= list[k+1];

      # Check that `pos' is inside the table,
      # and that its entry is assigned only once.
      if not pos in range then
        Error( "list entry ", list[k+1], " must lie in ", range );
      elif pos in entry[1] then
        Error( "position ", pos, " must occur at most once in <list>" );
      fi;

      # Check that the coefficients either fit to the zero element
      # or are rationals (with suitable denominators).
      if FamilyObj( val ) = Fam then
        if val <> zero then
          Add( entry[1], pos );
          Add( entry[2], val );
        fi;
      elif IsRat( val ) then
        if val <> 0 then
          Add( entry[1], pos );
          Add( entry[2], val * One( zero ) );
        fi;
      else
        Error( "list entry ", list[k], " does not fit to zero element" );
      fi;

    od;

    # Set the table entry.
    SortParallel( entry[1], entry[2] );
    T[i][j]:= Immutable( entry );

    # Add the value `T[j][i]' in the case of (anti-)symmetric tables.
    if   T[ Length(T) - 1 ] =  1 then
      T[j][i]:= T[i][j];
    elif T[ Length(T) - 1 ] = -1 then
      T[j][i]:= Immutable( [ entry[1], -entry[2] ] );
    fi;
end );


#############################################################################
##
#F  ReducedSCTable( <T>, <one> )
##
InstallGlobalFunction( ReducedSCTable, function( T, one )
    local new, n, i, j, entry;

    new:= [];
    n:= Length( T ) - 2;

    # Reduce the entries.
    for i in [ 1 .. n ] do
      new[i]:= [];
      for j in [ 1 .. n ] do
        entry:= T[i][j];
        entry:= [ Immutable( entry[1] ), entry[2] * one ];
        MakeImmutable( entry );
        new[i][j]:= entry;
      od;
    od;

    # Store zero coefficient and symmetry flag.
    new[ n+1 ]:= T[ n+1 ];
    new[ n+2 ]:= T[ n+2 ] * one;

    # Return the immutable new table.
    MakeImmutable( new );
    return new;
end );


#############################################################################
##
#F  GapInputSCTable( <T>, <varnam> )
##
InstallGlobalFunction( GapInputSCTable, function( T, varnam )
    local dim, str, lower, i, j, entry, k;

    # Initialize, and set the ranges for the loops.
    dim:= Length( T ) - 2;
    str:= Concatenation( varnam, ":= EmptySCTable( ",
                         String( dim ), ", ", String( T[ Length( T ) ] ) );
    lower:= [ 1 .. dim ];
    if   T[ dim+1 ] =  1 then
      Append( str, ", \"symmetric\"" );
    elif T[ dim+1 ] = -1 then
      Append( str, ", \"antisymmetric\"" );
    else
      lower:= ListWithIdenticalEntries( dim, 1 );
    fi;
    Append( str, " );\n" );

    # Fill up the table.
    for i in [ 1 .. dim ] do
      for j in [ lower[i] .. dim ] do
        entry:= T[i][j];
        if not IsEmpty( entry[1] ) then
          Append( str, "SetEntrySCTable( " );
          Append( str, varnam );
          Append( str, ", " );
          Append( str, String(i) );
          Append( str, ", " );
          Append( str, String(j) );
          Append( str, ", [" );
          for k in [ 1 .. Length( entry[1] )-1 ] do
            Append( str, String( entry[2][k] ) );
            Add( str, ',' );
            Append( str, String( entry[1][k] ) );
            Add( str, ',' );
          od;
          k:= Length( entry[1] );
          Append( str, String( entry[2][k] ) );
          Add( str, ',' );
          Append( str, String( entry[1][k] ) );
          Append( str, "] );\n" );
        fi;
      od;
    od;

    ConvertToStringRep( str );
    return str;
end );


#############################################################################
##
#F  IdentityFromSCTable( <T> )
##
InstallGlobalFunction( IdentityFromSCTable, function( T )
    local n,       # dimension of the underlying algebra
          equ,     # equation system to solve
          zero,    # zero of the field
          zerovec, # zero vector
          vec,     # right hand side of the equation system
          one,     # identity of the field
          i, j, k, # loop over rows of `equ'
          row,     # one row of the equation system
          Tpos,    #
          Tval,    #
          p,       #
          sol,
          sum;

    n:= Length( T ) - 2;
    zero:= T[ n+2 ];

    # If the table belongs to a trivial algebra,
    # the identity is equal to the zero.
    if n = 0 then
      return EmptyRowVector( FamilyObj( zero ) );
    fi;

    # Set up the equation system,
    # in row $i$ and column $(k-1)*n + j$ we have $c_{ijk}$.
    equ:= [];
    zerovec:= ListWithIdenticalEntries( n^2, zero );
    vec:= ShallowCopy( zerovec );
    one:= One( zero );

    for i in [ 1 .. n ] do
      row:= ShallowCopy( zerovec );
      for j in [ 1 .. n ] do
        Tpos:= T[i][j][1];
        Tval:= T[i][j][2];
        p:= (j-1)*n;
        for k in [ 1 .. Length( Tpos ) ] do
          row[ p + Tpos[k] ]:= Tval[k];
        od;
      od;
      Add( equ, row );
      vec[ (i-1)*n + i ]:= one;
    od;

    sol:= SolutionMat( equ, vec );

    # If we have a candidate and if the algebra is not known
    # to be commutative then check whether the candidate
    # acts trivially also from the right.
    if sol <> fail and T[ n+1 ] <> 1 then

      for j in [ 1 .. n ] do
        for k in [ 1 .. n ] do
          sum:= zero;
          for i in [ 1 .. n ] do
            Tpos:= T[j][i];
            p:= Position( Tpos[1], k );
#T cheaper !!!
            if p <> fail then
              sum:= sum + sol[i] * Tpos[2][p];
            fi;
          od;
          if ( j = k and sum <> one ) or ( j <> k and sum <> zero ) then
            return fail;
          fi;
        od;
      od;

    fi;

    # Return the result.
    return sol;
end );


#############################################################################
##
#F  QuotientFromSCTable( <T>, <num>, <den> )
##
##  We solve the equation system $<num> = x <den>$.
##  If no solution exists, `fail' is returned.
##
##  In terms of the basis $B$ with vectors $b_1, \ldots, b_n$ this means
##  for $<num> = \sum_{i=1}^n a_i b_i$,
##      $<den> = \sum_{i=1}^n c_i b_i$,
##      $x     = \sum_{i=1}^n x_i b_i$ that
##  $a_k = \sum_{i,j} c_i x_j c_{ijk}$ for all $k$.
##  Here $c_{ijk}$ denotes the structure constants w.r.t. $B$.
##  This means $a = x M$ with $M_{ik} = \sum_{j=1}^n c_{ijk} c_j$.
##
InstallGlobalFunction( QuotientFromSCTable, function( T, x, c )
    local M,        # matrix of the equation system
          n,        # dimension of the algebra
          zero,     # zero vector
          i, j,     # loop variables
          row,      # one row of `M'
          entry,    #
          val;      #

    M:= [];
    n:= Length( c );

    # If the algebra is zero dimensional,
    # the zero is also the identity and thus also its inverse.
    if n = 0 then
      return c;
    fi;

    zero:= ListWithIdenticalEntries( n, T[ Length( T ) ] );
    for i in [ 1 .. n ] do
      row:= ShallowCopy( zero );
      for j in [ 1 .. n ] do
        entry:= T[i][j];
        val:= c[j];
        row{ entry[1] }:= row{ entry[1] } + val * entry[2];
#T better!
      od;
      Add( M, row );
    od;

    # Return the quotient, or `fail'.
    return SolutionMat( M, x );
end );


#############################################################################
##
#F  TestJacobi( <T> )
##
##  We check whether for all $1 \leq m \leq n$ the equality
##  $\sum_{l=1}^n c_{jkl} c_{ilm} + c_{kil} c_{jlm} + c_{ijl} c_{klm} = 0$
##  holds.
##
InstallGlobalFunction( TestJacobi, function( T )
    local zero,           # the zero of the field
          n,              # dimension of the algebra
          i, j, k, m,     # loop variables
          cij, cki, cjk,  # structure constant vectors
          sum,
          t;

    zero:= T[ Length( T ) ];
    n:= Length( T ) - 2;

    for i in [ 1 .. n ] do
      for j in [ i+1 .. n ] do
        cij:= T[i][j];
        for k in [ j+1 .. n ] do
          cki:= T[k][i];
          cjk:= T[j][k];
          for m in [ 1 .. n ] do
            sum:= zero;
            for t in [ 1 .. Length( cjk[1] ) ] do
              sum:= sum + cjk[2][t] * SCTableEntry( T, i, cjk[1][t], m );
            od;
            for t in [ 1 .. Length( cki[1] ) ] do
              sum:= sum + cki[2][t] * SCTableEntry( T, j, cki[1][t], m );
            od;
            for t in [ 1 .. Length( cij[1] ) ] do
              sum:= sum + cij[2][t] * SCTableEntry( T, k, cij[1][t], m );
            od;
            if sum <> zero then
              return [ i, j, k ];
            fi;
          od;
        od;
      od;
    od;

    return true;
end );


#############################################################################
##
#M  MultiplicativeNeutralElement( <A> )
##
##  is the multiplicative neutral element of <A> if this exists,
##  otherwise is `fail'.
##
##  Let $(b_1, b_2, \ldots, b_n)$ be a basis of $A$, and $e$ the result of
##  `MultiplicativeNeutralElement( <A> )'.
##  Then $e = \sum_{i=1}^n a_i b_i$, and for $1 \leq k \leq n$ we have
##  $e \cdot b_j = b_j$, or equivalently
##  $\sum_{i=1}^n a_i b_i \cdot b_j = b_j$.
##  Define the structure constants by
##  $b_i \cdot b_j = \sum_{k=1}^n c_{ijk} b_k$.
##
##  Then $\sum_{i=1}^n a_i c_{ijk} = \delta_{jk}$ for $1 \leq k \leq n$.
##
##  This yields $n^2$ linear equations for the $n$ indeterminates $a_i$,
##  and a solution is a left identity.
##  For this we have to test whether it is also a right identity.
##
InstallMethod( MultiplicativeNeutralElement,
    [ IsFLMLOR and IsFiniteDimensional ],
    function( A )
    local B,       # basis of `A'
          one;     # result

    B:= Basis( A );
    one:= IdentityFromSCTable( StructureConstantsTable( B ) );
    if one <> fail then
      one:= LinearCombination( B, one );
    fi;
    return one;
    end );


#############################################################################
##
#M  IsAssociative( <A> )
##
##  We check whether the vectors of a basis satisfy the associativity law.
##  (Bilinearity of the multiplication is of course assumed.)
##
##  If $b_i \cdot b_j = \sum_{l=1}^n c_{ijl} b_l$ then we have
##  $b_i \cdot ( b_j \cdot b_k ) = ( b_i \cdot b_j ) \cdot b_k$
##  if and only if
##  $\sum_{l=1}^n c_{jkl} c_{ilm} = \sum_{l=1}^n c_{ijl} c_{lkm}$ for all
##  $1 \leq m \leq n$.
##
##  We check this equality for all $1 \leq i, j, k \leq n$.
##
InstallMethod( IsAssociative,
    "generic method for a (finite dimensional) FLMLOR",
    [ IsFLMLOR ],
    function( A )
    local T,            # structure constants table w.r.t. a basis of `A'
          zero,
          range,
          i, j, k, l, m,
          Ti,
          Tj,
          cijpos,
          cijval,
          cjkpos,
          cjkval,
          sum,
          x,
          pos;

    if not IsFiniteDimensional( A ) then
      TryNextMethod();
    fi;

    T:= StructureConstantsTable( Basis( A ) );
    zero:= Zero( LeftActingDomain( A ) );
    range:= [ 1 .. Length( T[1] ) ];
    for i in range do
      Ti:= T[i];
      for j in range do
        cijpos:= Ti[j][1];
        cijval:= Ti[j][2];
        Tj:= T[j];
        for k in range do
          cjkpos:= Tj[k][1];
          cjkval:= Tj[k][2];
          for m in range do
            sum:= zero;
            for l in [ 1 .. Length( cjkpos ) ] do
              x:= Ti[ cjkpos[l] ];
              pos:= Position( x[1], m );
              if pos <> fail then
                sum:= sum + cjkval[l] * x[2][ pos ];
              fi;
            od;
            for l in [ 1 .. Length( cijpos ) ] do

              x:= T[ cijpos[l] ][k];
              pos:= Position( x[1], m );
              if pos <> fail then
                sum:= sum - cijval[l] * x[2][ pos ];
              fi;
            od;
            if sum <> zero then
              # $i, j, k$ fail
              Info( InfoAlgebra, 2,
                    "IsAssociative fails for i = ", i, ", j = ", j,
                    ", k = ", k );
              return false;
            fi;
          od;
        od;
      od;
    od;
    return true;
    end );


#############################################################################
##
#M  IsAnticommutative( <A> )  . . . . . . . . . . . . .for a fin.-dim. FLMLOR
##
##  is `true' if the multiplication in <A> is anticommutative,
##  and `false' otherwise.
##
InstallMethod( IsAnticommutative,
    "generic method for a (finite dimensional) FLMLOR",
    [ IsFLMLOR ],
    function( A )
    local n,      # dimension of `A'
          T,      # table of structure constants for `A'
          zero,   # zero coefficient
          i, j;   # loop over rows and columns ot `T'

    if not IsFiniteDimensional( A ) then
      TryNextMethod();
    fi;

    n:= Dimension( A );
    T:= StructureConstantsTable( Basis( A ) );
    zero:= T[ n+2 ];
    for i in [ 2 .. n ] do
      for j in [ 1 .. i-1 ] do
        if    T[i][j][1] <> T[j][i][1]
           or ( not IsEmpty( T[i][j][1] )
                and PositionNot( T[i][j][2] + T[j][i][2], zero )
                        <= Length( T[i][j][2] ) ) then
          return false;
        fi;
      od;
    od;

    if Characteristic( A ) <> 2 then

      # The values on the diagonal must be zero.
      for i in [ 1 .. n ] do
        if not IsEmpty( T[i][i][1] ) then
          return false;
        fi;
      od;

    fi;

    return true;
    end );


#############################################################################
##
#M  IsCommutative( <A> )  . . . . . . . . . . . for finite dimensional FLMLOR
##
##  Check whether every basis vector commutes with every basis vector.
##
InstallMethod( IsCommutative,
    "generic method for a finite dimensional FLMLOR",
    [ IsFLMLOR ],
    IsCommutativeFromGenerators( GeneratorsOfVectorSpace ) );
#T use structure constants!


#############################################################################
##
#M  IsCommutative( <A> )  . . . . . . . . . . . . . for an associative FLMLOR
##
##  If <A> is associative then we can restrict the check to a smaller
##  equation system than that for arbitrary algebras, since we have to check
##  $x a = a x$ only for algebra generators $a$ and $x$, not for all vectors
##  of a basis.
##
InstallMethod( IsCommutative,
    "for an associative FLMLOR",
    [ IsFLMLOR and IsAssociative ],
    IsCommutativeFromGenerators( GeneratorsOfAlgebra ) );

InstallMethod( IsCommutative,
    "for an associative FLMLOR-with-one",
    [ IsFLMLORWithOne and IsAssociative ],
    IsCommutativeFromGenerators( GeneratorsOfAlgebraWithOne ) );


#############################################################################
##
#M  IsZeroSquaredRing( <A> )  . . . . . . . . for a finite dimensional FLMLOR
##
InstallMethod( IsZeroSquaredRing,
    "for a finite dimensional FLMLOR",
    [ IsFLMLOR ],
    function( A )

    if not IsAnticommutative( A ) then

      # Every zero squared ring is anticommutative.
      return false;

    elif ForAny( BasisVectors( Basis( A ) ),
                 x -> not IsZero( x*x ) ) then

      # If not all basis vectors are zero squared then we return `false'.
      return false;

    elif IsCommutative( LeftActingDomain( A ) ) then

      # If otherwise the left acting domain is commutative then we return
      # `true' because we know that <A> is anticommutative and the basis
      # vectors are zero squared.
      return true;

    else

      # Otherwise we give up.
      TryNextMethod();

    fi;
    end );


#############################################################################
##
#M  IsJacobianRing( <A> )
##
InstallMethod( IsJacobianRing,
    "for a (finite dimensional) FLMLOR",
    [ IsFLMLOR ],
    function( A )
    local n,   # dimension of `A'
          T,   # table of structure constants for `A'
          i;   # loop over the diagonal of `T'

    if not IsFiniteDimensional( A ) then
      TryNextMethod();
    fi;

    # In characteristic 2 we have to make sure that $a \* a = 0$.
#T really?
    T:= StructureConstantsTable( Basis( A ) );
    if Characteristic( A ) = 2 then
      n:= Dimension( A );
      for i in [ 1 .. n ] do
        if not IsEmpty( T[i][i][1] ) then
          return false;
        fi;
      od;
    fi;

    # Check the Jacobi identity $[a,[b,c]] + [b,[c,a]] + [c,[a,b]] = 0$.
    return TestJacobi( T ) = true;
    end );


#############################################################################
##
#M  Intersection2( <A1>, <A2> ) . . . . . . . . . intersection of two FLMLORs
##
InstallMethod( Intersection2,
    "generic method for two FLMLORs",
    IsIdenticalObj,
    [ IsFLMLOR, IsFLMLOR ],
    Intersection2Spaces( AsFLMLOR, SubFLMLORNC, FLMLOR ) );


#############################################################################
##
#M  Intersection2( <A1>, <A2> ) . . . .  intersection of two FLMLORs-with-one
##
InstallMethod( Intersection2,
    "generic method for two FLMLORs-with-one",
    IsIdenticalObj,
    [ IsFLMLORWithOne, IsFLMLORWithOne ],
    Intersection2Spaces( AsFLMLORWithOne, SubFLMLORWithOneNC,
                         FLMLORWithOne ) );


#############################################################################
##
#M  \/( <A>, <I> )  . . . . . . . . . . . .  factor of an algebra by an ideal
#M  \/( <A>, <relators> ) . . . . . . . . .  factor of an algebra by an ideal
##
##  is the factor algebra of the finite dimensional algebra <A> modulo
##  the ideal <I> or the ideal spanned by the collection <relators>.
##
InstallOtherMethod( \/,
    "for FLMLOR and collection",
    IsIdenticalObj,
    [ IsFLMLOR, IsCollection ],
    function( A, relators )
    if IsFLMLOR( relators ) then
      TryNextMethod();
    else
      return A / TwoSidedIdealByGenerators( A, relators );
    fi;
    end );

InstallOtherMethod( \/,
    "for FLMLOR and empty list",
    [ IsFLMLOR, IsList and IsEmpty ],
    function( A, empty )
    # `NaturalHomomorphismByIdeal( A, TrivialSubFLMLOR( A ) )' is the
    # identity mapping on `A', and `ImagesSource' of it yields `A'.
    return A;
    end );

InstallOtherMethod( \/,
    "generic method for two FLMLORs",
    IsIdenticalObj,
    [ IsFLMLOR, IsFLMLOR ],
    function( A, I )
    return ImagesSource( NaturalHomomorphismByIdeal( A, I ) );
    end );


#############################################################################
##
#M  TrivialSubadditiveMagmaWithZero( <A> )  . . . . . . . . . .  for a FLMLOR
##
InstallMethod( TrivialSubadditiveMagmaWithZero,
    "for a FLMLOR",
    [ IsFLMLOR ],
    A -> SubFLMLORNC( A, [] ) );


#############################################################################
##
#M  AsFLMLOR( <R>, <D> )  . . view a collection as a FLMLOR over the ring <R>
##
InstallMethod( AsFLMLOR,
    "for a ring and a collection",
    [ IsRing, IsCollection ],
    function( F, D )
    local A, L;

    D:= AsSSortedList( D );
    L:= ShallowCopy( D );
    A:= TrivialSubFLMLOR( AsFLMLOR( F, D ) );
    SubtractSet( L, AsSSortedList( A ) );
    while 0 < Length(L)  do
      A:= ClosureLeftOperatorRing( A, L[1] );
#T call explicit function that maintains an elements list?
      SubtractSet( L, AsSSortedList( A ) );
    od;
    if Length( AsList( A ) ) <> Length( D )  then
      return fail;
    fi;
    A:= FLMLOR( F, GeneratorsOfLeftOperatorRing( A ), Zero( D[1] ) );
    SetAsSSortedList( A, D );
    SetSize( A, Length( D ) );
    SetIsFinite( A, true );
#T ?

    # Return the FLMLOR.
    return A;
    end );


#############################################################################
##
#M  AsFLMLOR( <F>, <V> )  . . view a left module as FLMLOR over the field <F>
##
##  is an algebra over <F> that is equal (as set) to <V>.
##  For that, perhaps the field of <A> has to be changed before
##  getting the correct list of generators.
##
InstallMethod( AsFLMLOR,
    "for a division ring and a free left module",
    [ IsDivisionRing, IsFreeLeftModule ],
    function( F, V )
    local L, A;

    if   LeftActingDomain( V ) = F then

      A:= FLMLOR( F, GeneratorsOfLeftModule( V ) );
      if A <> V then
        return fail;
      fi;
      if HasBasis( V ) then
        SetBasis( A, Basis( V ) );
      fi;

    elif IsTrivial( V ) then

      # We need the zero.
      A:= FLMLOR( F, [], Zero( V ) );

    elif IsSubset( LeftActingDomain( V ), F ) then

      # Make sure that the field change does not change the elements.
      L:= BasisVectors( Basis( AsField( F, LeftActingDomain(V) ) ) );
      L:= Concatenation( List( L, x -> List( GeneratorsOfLeftModule( V ),
                                             y -> x * y ) ) );
      A:= FLMLOR( F, L );
      if A <> V then
        return fail;
      fi;

    elif IsSubset( F, LeftActingDomain( V ) ) then

      # Make sure that the field change does not change the elements.
      L:= BasisVectors( Basis( AsField( LeftActingDomain(V), F ) ) );
      if ForAny( L, x -> ForAny( GeneratorsOfLeftModule( V ),
                                 y -> not x * y in V ) ) then
        return fail;
      fi;
      A:= FLMLOR( F, GeneratorsOfLeftModule( V ) );
      if A <> V then
        return fail;
      fi;

    else

      V:= AsFLMLOR( Intersection( F, LeftActingDomain( V ) ), V );
      return AsFLMLOR( F, V );

    fi;

    UseIsomorphismRelation( V, A );
    UseSubsetRelation( V, A );

    return A;
    end );


#############################################################################
##
#M  AsFLMLOR( <F>, <A> )  . . . view an algebra as algebra over the field <F>
##
##  is an algebra over <F> that is equal (as set) to <D>.
##  For that, perhaps the field of <A> has to be changed before
##  getting the correct list of generators.
##
InstallMethod( AsFLMLOR,
    "for a division ring and an algebra",
    [ IsDivisionRing, IsFLMLOR ],
    function( F, D )
    local L, A;

    if   LeftActingDomain( D ) = F then

      return D;

    elif IsTrivial( D ) then

      # We need the zero.
      A:= FLMLOR( F, [], Zero( D ) );

    elif IsSubset( LeftActingDomain( D ), F ) then

      # Make sure that the field change does not change the elements.
      L:= BasisVectors( Basis( AsField( F, LeftActingDomain(D) ) ) );
      L:= Concatenation( List( L, x -> List( GeneratorsOfAlgebra( D ),
                                             y -> x * y ) ) );
      A:= FLMLOR( F, L );

    elif IsSubset( F, LeftActingDomain( D ) ) then

      # Make sure that the field change does not change the elements.
      L:= BasisVectors( Basis( AsField( LeftActingDomain(D), F ) ) );
      if ForAny( L, x -> ForAny( GeneratorsOfAlgebra( D ),
                                 y -> not x * y in D ) ) then
        return fail;
      fi;
      A:= FLMLOR( F, GeneratorsOfAlgebra( D ) );

    else

      D:= AsFLMLOR( Intersection( F, LeftActingDomain( D ) ), D );
      return AsFLMLOR( F, D );

    fi;

    UseIsomorphismRelation( D, A );
    UseSubsetRelation( D, A );
    return A;
    end );


#############################################################################
##
#M  AsFLMLORWithOne( <R>, <D> ) . . . . . . view a coll. as a FLMLOR-with-one
##
InstallMethod( AsFLMLORWithOne,
    "for a ring and a collection",
    [ IsRing, IsCollection ],
    function( F, D )
    return AsFLMLORWithOne( AsFLMLOR( F, D ) );
    end );


#############################################################################
##
#M  AsFLMLORWithOne( <F>, <V> ) . .  view a left module as a algebra-with-one
##
InstallMethod( AsFLMLORWithOne,
    "for a division ring and a free left module",
    [ IsDivisionRing, IsFreeLeftModule ],
    function( F, V )
    local L, A;

    # Check that `V' contains the identity.
    if One( V ) = fail then

      return fail;

    elif LeftActingDomain( V ) = F then

      A:= FLMLORWithOne( F, GeneratorsOfLeftModule( V ) );
      if A <> V then
        return fail;
      fi;

      # Left module generators and basis are maintained.
      if HasGeneratorsOfLeftModule( V ) then
        SetGeneratorsOfLeftModule( A, GeneratorsOfLeftModule( V ) );
      fi;
      if HasBasis( V ) then
        SetBasis( A, Basis( V ) );
      fi;

    elif IsSubset( LeftActingDomain( V ), F ) then

      # Make sure that the field change does not change the elements.
      L:= BasisVectors( Basis( AsField( F, LeftActingDomain(V) ) ) );
      L:= Concatenation( List( L, x -> List( GeneratorsOfLeftModule( V ),
                                             y -> x * y ) ) );
      A:= FLMLORWithOne( F, L );
      if A <> V then
        return fail;
      fi;

    elif IsSubset( F, LeftActingDomain( V ) ) then

      # Make sure that the field change does not change the elements.
      L:= BasisVectors( Basis( AsField( LeftActingDomain(V), F ) ) );
      if ForAny( L, x -> ForAny( GeneratorsOfLeftModule( V ),
                                 y -> not x * y in V ) ) then
        return fail;
      fi;
      A:= FLMLORWithOne( F, GeneratorsOfLeftModule( V ) );
      if A <> V then
        return fail;
      fi;

    else

      # Note that we need not use the isomorphism and subset relations
      # (see below) because this is the task of the calls to
      # `AsAlgebraWithOne'.
      A:= AsAlgebraWithOne( Intersection( F, LeftActingDomain( V ) ), V );
      return AsAlgebraWithOne( F, A );

    fi;

    UseIsomorphismRelation( V, A );
    UseSubsetRelation( V, A );
    return A;
    end );


#############################################################################
##
#M  AsFLMLORWithOne( <F>, <D> ) . . . . view an algebra as a algebra-with-one
##
InstallMethod( AsFLMLORWithOne,
    "for a division ring and an algebra",
    [ IsDivisionRing, IsFLMLOR ],
    function( F, D )
    local L, A;

    # Check that `D' contains the identity.
    if One( D ) = fail then

      return fail;

    elif LeftActingDomain( D ) = F then

      A:= FLMLORWithOne( F, GeneratorsOfLeftOperatorRing( D ) );

      # Left module generators and basis are maintained.
      if HasGeneratorsOfLeftModule( D ) then
        SetGeneratorsOfLeftModule( A, GeneratorsOfLeftModule( D ) );
      fi;
      if HasBasis( D ) then
        SetBasis( A, Basis( D ) );
      fi;

    elif IsSubset( LeftActingDomain( D ), F ) then

      # Make sure that the field change does not change the elements.
      L:= BasisVectors( Basis( AsField( F, LeftActingDomain(D) ) ) );
      L:= Concatenation( List( L, x -> List( GeneratorsOfAlgebra( D ),
                                             y -> x * y ) ) );
      A:= FLMLORWithOne( F, L );

    elif IsSubset( F, LeftActingDomain( D ) ) then

      # Make sure that the field change does not change the elements.
      L:= BasisVectors( Basis( AsField( LeftActingDomain(D), F ) ) );
      if ForAny( L, x -> ForAny( GeneratorsOfAlgebra( D ),
                                 y -> not x * y in D ) ) then
        return fail;
      fi;
      A:= FLMLORWithOne( F, GeneratorsOfLeftOperatorRing( D ) );

    else

      # Note that we need not use the isomorphism and subset relations
      # (see below) because this is the task of the calls to
      # `AsAlgebraWithOne'.
      D:= AsAlgebraWithOne( Intersection( F, LeftActingDomain( D ) ), D );
      return AsAlgebraWithOne( F, D );

    fi;

    UseIsomorphismRelation( D, A );
    UseSubsetRelation( D, A );
    return A;
    end );


#############################################################################
##
#M  AsFLMLORWithOne( <F>, <D> ) . . view an alg.-with-one as an alg.-with-one
##
InstallMethod( AsFLMLORWithOne,
    "for a division ring and a algebra-with-one",
    [ IsDivisionRing, IsFLMLORWithOne ],
    function( F, D )
    local L, A;

    if   LeftActingDomain( D ) = F then

      return D;

    elif IsSubset( LeftActingDomain( D ), F ) then

      # Make sure that the field change does not change the elements.
      L:= BasisVectors( Basis( AsField( F, LeftActingDomain(D) ) ) );
      L:= Concatenation( List( L, x -> List( GeneratorsOfAlgebra( D ),
                                             y -> x * y ) ) );
      A:= AlgebraWithOne( F, L );

    elif IsSubset( F, LeftActingDomain( D ) ) then

      # Make sure that the field change does not change the elements.
      L:= BasisVectors( Basis( AsField( LeftActingDomain(D), F ) ) );
      if ForAny( L, x -> ForAny( GeneratorsOfAlgebra( D ),
                                 y -> not x * y in D ) ) then
        return fail;
      fi;
      A:= AlgebraWithOne( F, GeneratorsOfAlgebra( D ) );

    else

      # Note that we need not use the isomorphism and subset relations
      # (see below) because this is the task of the calls to
      # `AsAlgebraWithOne'.
      D:= AsAlgebraWithOne( Intersection( F, LeftActingDomain( D ) ), D );
      return AsAlgebraWithOne( F, D );

    fi;

    UseIsomorphismRelation( D, A );
    UseSubsetRelation( D, A );
    return A;
    end );


#############################################################################
##
#M  ClosureLeftOperatorRing( <A>, <a> ) . . . . . . . closure with an element
##
InstallMethod( ClosureLeftOperatorRing,
    "for a FLMLOR and a ring element",
#T why not a general function for any left operator ring?
#T (need `LeftOperatorRingByGenerators' ?)
    IsCollsElms,
    [ IsFLMLOR, IsRingElement ],
    function( A, a )

    # if possible test if the element lies in the ring already,
    if     HasGeneratorsOfLeftOperatorRing( A )
       and a in GeneratorsOfLeftOperatorRing( A ) then
      return A;

    # otherwise make a new left operator ring
    else
      return FLMLOR( LeftActingDomain( A ),
                 Concatenation( GeneratorsOfLeftOperatorRing( A ), [ a ] ) );
    fi;
    end );

InstallMethod( ClosureLeftOperatorRing,
    "for an FLMLOR with basis, and a ring element",
    IsCollsElms,
    [ IsFLMLOR and HasBasis, IsRingElement ],
    function( A, a )

    # test if the element lies in the FLMLOR already,
    if a in A then
      return A;

    # otherwise make a new FLMLOR
    else
      return FLMLOR( LeftActingDomain( A ),
#T FLMLORByGenerators?
                 Concatenation( BasisVectors( Basis( A ) ), [ a ] ),
                 "basis" );
    fi;
    end );

InstallMethod( ClosureLeftOperatorRing,
    "for a FLMLOR-with-one and a ring element",
    IsCollsElms,
    [ IsFLMLORWithOne, IsRingElement ],
    function( A, a )

    # if possible test if the element lies in the ring already,
    if a in GeneratorsOfLeftOperatorRingWithOne( A ) then
      return A;

    # otherwise make a new left operator ring-with-one
    else
      return FLMLORWithOne( LeftActingDomain( A ),
                 Concatenation( GeneratorsOfLeftOperatorRingWithOne( A ),
                                [ a ] ) );
    fi;
    end );

InstallMethod( ClosureLeftOperatorRing,
    "for a FLMLOR-with-one with basis, and a ring element",
    IsCollsElms,
    [ IsFLMLORWithOne and HasBasis, IsRingElement ],
    function( A, a )

    # test if the element lies in the FLMLOR already,
    if a in A then
        return A;

    # otherwise make a new FLMLOR-with-one
    else
      return FLMLORWithOne( LeftActingDomain( A ),
                 Concatenation( BasisVectors( Basis( A ) ), [ a ] ),
                 "basis" );
    fi;
    end );

InstallMethod( ClosureLeftOperatorRing,
    "for a FLMLOR containing the whole family, and a ring element",
    IsCollsElms,
    [ IsFLMLOR and IsWholeFamily, IsRingElement ],
    SUM_FLAGS, # this is better than everything else
    function( A, a )
    return A;
    end );


#############################################################################
##
#M  ClosureLeftOperatorRing( <A>, <U> ) .  closure of two left operator rings
##
InstallMethod( ClosureLeftOperatorRing,
    "for two left operator rings",
    IsIdenticalObj,
    [ IsLeftOperatorRing, IsLeftOperatorRing ],
    function( A, S )
    local   g;          # one generator

    for g in GeneratorsOfLeftOperatorRing( S ) do
      A := ClosureLeftOperatorRing( A, g );
    od;
    return A;
    end );

InstallMethod( ClosureLeftOperatorRing,
    "for two left operator rings-with-one",
    IsIdenticalObj,
    [ IsLeftOperatorRingWithOne, IsLeftOperatorRingWithOne ],
    function( A, S )
    local   g;          # one generator

    for g in GeneratorsOfLeftOperatorRingWithOne( S ) do
      A := ClosureLeftOperatorRing( A, g );
    od;
    return A;
    end );

InstallMethod( ClosureLeftOperatorRing,
    "for a left op. ring cont. the whole family, and a collection",
    IsIdenticalObj,
    [ IsLeftOperatorRing and IsWholeFamily, IsCollection ],
    SUM_FLAGS, # this is better than everything else
    function( A, S )
    return A;
    end );


#############################################################################
##
#M  ClosureLeftOperatorRing( <A>, <list> )  . . . .  closure of left op. ring
##
InstallMethod( ClosureLeftOperatorRing,
    "for left operator ring and list of elements",
    IsIdenticalObj,
    [ IsLeftOperatorRing, IsCollection ],
    function( A, list )
    local   g;          # one generator
    for g in list do
      A:= ClosureLeftOperatorRing( A, g );
    od;
    return A;
    end );


#############################################################################
##
#F  MutableBasisOfClosureUnderAction( <F>, <Agens>, <from>, <init>, <opr>,
#F                                    <zero>, <maxdim> )
##
##  This function is used to compute bases of finite dimensional ideals $I$
##  in *associative* algebras that are given by ideal generators $J$ and
##  (generators of) the acting algebra $A$.
##  An important special case is that of the algebra $A$ itself, given by
##  algebra generators.
##
##  The algorithm assumes that it is possible to deal with mutable bases of
##  vector spaces generated by elements of $A$.
##  It proceeds as follows.
##
##  Let $A$ be a finite dimensional algebra over the ring $F$,
##  and $I$ a two-sided ideal in $A$ that is generated (as a two-sided ideal)
##  by the set $J$.
##  (For the cases of one-sided ideals, see below.)
##
##  Let $S$ be a set of algebra generators of $A$.
##  The identity of $A$, if exists, need not be contained in $S$.
##
##  Each element $x$ in $I$ can be written as a linear combination of
##  products $j a_1 a_2 \cdots a_n$, with $j \in J$ and $a_i \in S$ for
##  $1 \leq i \leq n$.
##  Let $l(x)$ denote the minimum of the largest $n$ for involved words
##  $a_1 a_2 \cdots a_n$, taken over all possible expressions for $x$.
##
##  Define $I_i = \{ x \in I \mid l(x) \leq i \}$.
##  Then $I_i$ is an $F$-space, $A_0 = \langle J \rangle_F$,
##  and $I_0 \< I_1 \< I_2 \< \cdots I_k \< I_{k+1} \< \ldots$
##  is an ascending chain that eventually reaches $I$.
##  For $i > 0$ we have
##  $I_{i+1} = \langle I_i\cup\bigcup_{s\in S} ( I_i s\cup s I_i )\rangle_F$.
##
##  (*Note* that the computation of the $I_i$ gives us the smallest value $k$
##  such that every element is a linear combination of words in terms of the
##  algebra generators, of maximal length $k$.)
##
InstallGlobalFunction( MutableBasisOfClosureUnderAction,
    function( F, Agens, from, init, opr, zero, maxdim )
    local MB,        # mutable basis, result
          gen,       # loop over generators
          v,         #
          dim,       # dimension of the actual left module
          right,     # `true' if we have to multiply from the right
          left;      # `true' if we have to multiply from the left

    # Get the side(s) from where to multiply.
    left  := true;
    right := true;
    if   from = "left" then
      right:= false;
    elif from = "right" then
      left:= false;
    fi;

    # $I_0$
    MB  := MutableBasis( F, init, zero );
    dim := 0;

    while dim < NrBasisVectors( MB ) and dim < maxdim do

      # `MB' is a mutable basis of $I_i$.
      dim:= NrBasisVectors( MB );

      if right then

        # Compute $I^{\prime}_i = I_i + \sum_{s \in S} I_i s$.
        for gen in Agens do
          for v in BasisVectors( MB ) do
            CloseMutableBasis( MB, opr( v, gen ) );
          od;
        od;

      fi;
      if left then

        # Compute $I_i + \sum_{s \in S} s I_i$
        # resp. $I^{\prime}_i + \sum_{s \in S} s I_i$.
        for gen in Agens do
          for v in BasisVectors( MB ) do
            CloseMutableBasis( MB, opr( gen, v ) );
          od;
        od;

      fi;

    od;

    # Return the mutable basis.
    return MB;
end );


#############################################################################
##
#F  MutableBasisOfNonassociativeAlgebra( <F>, <Agens>, <zero>, <maxdim> )
##
InstallGlobalFunction( MutableBasisOfNonassociativeAlgebra,
    function( F, Agens, zero, maxdim )
    local MB,        # mutable basis, result
          dim,       # dimension of the current left module
          bv,        # current basis vectors
          v, w;      # loop over basis vectors found already

    MB  := MutableBasis( F, Agens, zero );
    dim := 0;

    while dim < NrBasisVectors( MB ) and dim < maxdim do

      dim := NrBasisVectors( MB );
      bv  := BasisVectors( MB );

      for v in bv do
        for w in bv do
          CloseMutableBasis( MB, v * w );
          CloseMutableBasis( MB, w * v );
        od;
      od;

    od;

    # Return the mutable basis.
    return MB;
end );


#############################################################################
##
#F  MutableBasisOfIdealInNonassociativeAlgebra( <F>, <Vgens>, <Igens>,
#F                                              <zero>, <from>, <maxdim> )
##
InstallGlobalFunction( MutableBasisOfIdealInNonassociativeAlgebra,
    function( F, Vgens, Igens, zero, from, maxdim )
    local MB,        # mutable basis, result
          dim,       # dimension of the current left module
          bv,        # current basis vectors
          right,     # `true' if we have to multiply from the right
          left,      # `true' if we have to multiply from the left
          v, gen;    # loop over basis vectors found already

    # Get the side(s) from where to multiply.
    left  := true;
    right := true;
    if   from = "left" then
      right:= false;
    elif from = "right" then
      left:= false;
    fi;

    dim := 0;
    MB  := MutableBasis( F, Igens, zero );

    while dim < NrBasisVectors( MB ) and dim < maxdim do

      dim := NrBasisVectors( MB );
      bv  := BasisVectors( MB );

      for v in bv do
        for gen in Vgens do

          if left then
            CloseMutableBasis( MB, gen * v );
          fi;
          if right then
            CloseMutableBasis( MB, v * gen );
          fi;
        od;
      od;

    od;

    # Return the mutable basis.
    return MB;
end );


#############################################################################
##
#M  IsSubset( <A>, <B> )  . . . . . . . . . . . .  test for subset of FLMLORs
##
##  These methods are preferable to that for free left modules because they
##  use algebra generators.
##
##  We assume that generators of an extension of the left acting domains can
##  be computed if they are fields; note that infinite field extensions do
##  not (yet) occur in {\GAP} as `LeftActingDomain' values.
##
InstallMethod( IsSubset,
    "for two FLMLORs",
    IsIdenticalObj,
    [ IsFLMLOR, IsFLMLOR ],
    function( D1, D2 )
    local F1, F2;

    F1:= LeftActingDomain( D1 );
    F2:= LeftActingDomain( D2 );
    if not ( HasIsDivisionRing( F1 ) and IsDivisionRing( F1 ) and 
             HasIsDivisionRing( F2 ) and IsDivisionRing( F2 ) ) then
      TryNextMethod();
    fi;

    # catch trivial case
    if IsSubset(GeneratorsOfLeftOperatorRing(D1),
                  GeneratorsOfLeftOperatorRing(D2)) then
	return true;
    fi;
    return IsSubset( D1, GeneratorsOverIntersection( D2,
                             GeneratorsOfLeftOperatorRing( D2 ),
                             F2, F1 ) );
    end );

InstallMethod( IsSubset,
    "for two FLMLORs-with-one",
    IsIdenticalObj,
    [ IsFLMLORWithOne, IsFLMLORWithOne ],
    function( D1, D2 )
    local F1, F2;

    F1:= LeftActingDomain( D1 );
    F2:= LeftActingDomain( D2 );
    if not ( HasIsDivisionRing( F1 ) and IsDivisionRing( F1 ) and 
             HasIsDivisionRing( F2 ) and IsDivisionRing( F2 ) ) then
      TryNextMethod();
    fi;
    return IsSubset( D1, GeneratorsOverIntersection( D2,
                             GeneratorsOfLeftOperatorRingWithOne( D2 ),
                             F2, F1 ) );
    end );


#############################################################################
##
#M  ViewObj( <A> )  . . . . . . . . . . . . . . . . . . . . . . view a FLMLOR
##
##  print left acting domain, if known also dimension or no. of generators
##
InstallMethod( ViewObj,
    "for a FLMLOR",
    [ IsFLMLOR ],
    function( A )
    Print( "<free left module over ", LeftActingDomain( A ), ", and ring>" );
    end );

InstallMethod( ViewObj,
    "for a FLMLOR with known dimension",
    [ IsFLMLOR and HasDimension ], 1,   # override method requiring gens.
    function( A )
    Print( "<free left module of dimension ", Dimension( A ),
           " over ", LeftActingDomain( A ), ", and ring>" );
    end );

InstallMethod( ViewObj,
    "for a FLMLOR with known generators",
    [ IsFLMLOR and HasGeneratorsOfAlgebra ],
    function( A )
    Print( "<free left module over ", LeftActingDomain( A ),
           ", and ring, with ",
           Length( GeneratorsOfFLMLOR( A ) ), " generators>" );
    end );


#############################################################################
##
#M  PrintObj( <A> ) . . . . . . . . . . . . . . . . . . . . .  print a FLMLOR
##
InstallMethod( PrintObj,
    "for a FLMLOR",
    [ IsFLMLOR ],
    function( A )
    Print( "FLMLOR( ", LeftActingDomain( A ), ", ... )" );
    end );

InstallMethod( PrintObj,
    "for a FLMLOR with known generators",
    [ IsFLMLOR and HasGeneratorsOfFLMLOR ],
    function( A )
    if IsEmpty( GeneratorsOfFLMLOR( A ) ) then
      Print( "FLMLOR( ", LeftActingDomain( A ), ", [], ", Zero( A ), " )" );
    else
      Print( "FLMLOR( ", LeftActingDomain( A ), ", ",
             GeneratorsOfFLMLOR( A ), " )" );
    fi;
    end );


#############################################################################
##
#M  ViewObj( <A> )  . . . . . . . . . . . . . . . . .  view a FLMLOR-with-one
##
##  print left acting domain, if known also dimension or no. of generators
##
InstallMethod( ViewObj,
    "for a FLMLOR-with-one",
    [ IsFLMLORWithOne ],
    function( A )
    Print( "<free left module over ", LeftActingDomain( A ),
           ", and ring-with-one>" );
    end );

InstallMethod( ViewObj,
    "for a FLMLOR-with-one with known dimension",
    [ IsFLMLORWithOne and HasDimension ], 1,   # override method requ. gens.
    function( A )
    Print( "<free left module of dimension ", Dimension( A ),
           " over ", LeftActingDomain( A ), ", and ring-with-one>" );
    end );

InstallMethod( ViewObj,
    "for a FLMLOR-with-one with known generators",
    [ IsFLMLORWithOne and HasGeneratorsOfFLMLORWithOne ],
    function( A )
    Print( "<free left module over ", LeftActingDomain( A ),
           ", and ring-with-one, with ",
           Length( GeneratorsOfAlgebraWithOne( A ) ), " generators>" );

    end );


#############################################################################
##
#M  PrintObj( <A> ) . . . . . . . . . . . . . . . . . print a FLMLOR-with-one
##
InstallMethod( PrintObj,
    "for a FLMLOR-with-one",
    [ IsFLMLORWithOne ],
    function( A )
    Print( "FLMLORWithOne( ", LeftActingDomain( A ), ", ... )" );
    end );

InstallMethod( PrintObj,
    "for a FLMLOR-with-one with known generators",
    [ IsFLMLORWithOne and HasGeneratorsOfFLMLOR ],
    function( A )
    if IsEmpty( GeneratorsOfFLMLORWithOne( A ) ) then
      Print( "FLMLORWithOne( ", LeftActingDomain( A ), ", [], ",
             Zero( A ), " )" );
    else
      Print( "FLMLORWithOne( ", LeftActingDomain( A ), ", ",
             GeneratorsOfFLMLORWithOne( A ), " )" );
    fi;
    end );


#############################################################################
##
#M  ViewObj( <A> )  . . . . . . . . . . . . . . . . . . . . . view an algebra
##
##  print left acting domain, if known also dimension or no. of generators
##
InstallMethod( ViewObj,
    "for an algebra",
    [ IsAlgebra ],
    function( A )
    Print( "<algebra over ", LeftActingDomain( A ), ">" );
    end );

InstallMethod( ViewObj,
    "for an algebra with known dimension",
    [ IsAlgebra and HasDimension ], 1,   # override method requiring gens.
    function( A )
    Print( "<algebra of dimension ", Dimension( A ),
           " over ", LeftActingDomain( A ), ">" );
    end );

InstallMethod( ViewObj,
    "for an algebra with known generators",
    [ IsAlgebra and HasGeneratorsOfAlgebra ],
    function( A )
    Print( "<algebra over ", LeftActingDomain( A ), ", with ",
           Length( GeneratorsOfAlgebra( A ) ), " generators>" );
    end );


#############################################################################
##
#M  PrintObj( <A> ) . . . . . . . . . . . . . . . . . . . .  print an algebra
##
InstallMethod( PrintObj,
    "for an algebra",
    [ IsAlgebra ],
    function( A )
    Print( "Algebra( ", LeftActingDomain( A ), ", ... )" );
    end );

InstallMethod( PrintObj,
    "for an algebra with known generators",
    [ IsAlgebra and HasGeneratorsOfAlgebra ],
    function( A )
    if IsEmpty( GeneratorsOfAlgebra( A ) ) then
      Print( "Algebra( ", LeftActingDomain( A ), ", [], ", Zero( A ), " )" );
    else
      Print( "Algebra( ", LeftActingDomain( A ), ", ",
             GeneratorsOfAlgebra( A ), " )" );
    fi;
    end );


#############################################################################
##
#M  ViewObj( <A> )  . . . . . . . . . . . . . . . .  view an algebra-with-one
##
##  print left acting domain, if known also dimension or no. of generators
##
InstallMethod( ViewObj, "for an algebra-with-one", [ IsAlgebraWithOne ],
function( A )
  if IsIdenticalObj(A,LeftActingDomain(A)) then
    Print( "<algebra-with-one over itself>" );
  else
    Print( "<algebra-with-one over ", LeftActingDomain( A ), ">" );
  fi;
end );

InstallMethod( ViewObj,
    "for an algebra-with-one with known dimension",
    [ IsAlgebraWithOne and HasDimension ], 1,   # override method requ. gens.
    function( A )
    Print( "<algebra-with-one of dimension ", Dimension( A ),
           " over ", LeftActingDomain( A ), ">" );
    end );

InstallMethod( ViewObj,
    "for an algebra-with-one with known generators",
    [ IsAlgebraWithOne and HasGeneratorsOfAlgebraWithOne ],
    function( A )
    Print( "<algebra-with-one over ", LeftActingDomain( A ), ", with ",
           Length( GeneratorsOfAlgebraWithOne( A ) ), " generators>" );
    end );


#############################################################################
##
#M  PrintObj( <A> ) . . . . . . . . . . . . . . . . print an algebra-with-one
##
InstallMethod( PrintObj,
    "for an algebra-with-one",
    [ IsAlgebraWithOne ],
    function( A )
    Print( "AlgebraWithOne( ", LeftActingDomain( A ), ", ... )" );
    end );

InstallMethod( PrintObj,
    "for an algebra-with-one with known generators",
    [ IsAlgebraWithOne and HasGeneratorsOfAlgebra ],
    function( A )
    if IsEmpty( GeneratorsOfAlgebraWithOne( A ) ) then
      Print( "AlgebraWithOne( ", LeftActingDomain( A ), ", [], ",
             Zero( A ), " )" );
    else
      Print( "AlgebraWithOne( ", LeftActingDomain( A ), ", ",
             GeneratorsOfAlgebraWithOne( A ), " )" );
    fi;
    end );


#############################################################################
##
#M  ViewObj( <A> )  . . . . . . . . . . . . . . . . . . view a Lie algebra
##
##  print left acting domain, if known also dimension or no. of generators
##
InstallMethod( ViewObj,
    "for a Lie algebra",
    [ IsLieAlgebra ],
    function( A )
    Print( "<Lie algebra over ", LeftActingDomain( A ), ">" );
    end );

InstallMethod( ViewObj,
    "for a Lie algebra with known dimension",
    [ IsLieAlgebra and HasDimension ], 1,       # override method requ. gens.
    function( A )
    Print( "<Lie algebra of dimension ", Dimension( A ),
           " over ", LeftActingDomain( A ), ">" );
    end );

InstallMethod( ViewObj,
    "for a Lie algebra with known generators",
    [ IsLieAlgebra and HasGeneratorsOfAlgebra ],
    function( A )
    Print( "<Lie algebra over ", LeftActingDomain( A ), ", with ",
           Length( GeneratorsOfAlgebra( A ) ), " generators>" );
end );


#############################################################################
##
#M  AsSubalgebra(<A>, <U>) . view an algebra as subalgebra of another algebra
##
InstallMethod( AsSubalgebra,
    "for two algebras",
    IsIdenticalObj,
    [ IsAlgebra, IsAlgebra ],
    function( A, U )
    local samecoeffs, S;
    if not IsSubset( A, U ) then
      return fail;
    fi;

    # Construct the generators list.
    samecoeffs:= LeftActingDomain( A ) = LeftActingDomain( U );
    if not samecoeffs then
      U:= AsAlgebra( LeftActingDomain( A ), U );
    fi;

    # Construct the subalgebra.
    S:= SubalgebraNC( A, GeneratorsOfAlgebra( U ) );

    # Maintain useful information.
    UseIsomorphismRelation( U, S );
    UseSubsetRelation( U, S );
    if samecoeffs and HasDimension( U ) then
      SetDimension( S, Dimension( U ) );
    fi;

    # Return the subalgebra.
    return S;
    end );

InstallMethod( AsSubalgebra,
    "for an algebra and an algebra-with-one",
    IsIdenticalObj,
    [ IsAlgebra, IsAlgebraWithOne ],
    function( A, U )
    local samecoeffs, S;
    if not IsSubset( A, U ) then
      return fail;
    fi;

    # Construct the generators list.
    samecoeffs:= LeftActingDomain( A ) = LeftActingDomain( U );
    if not samecoeffs then
      U:= AsAlgebraWithOne( LeftActingDomain( A ), U );
    fi;

    # Construct the subalgebra.
    S:= SubalgebraWithOneNC( A, GeneratorsOfAlgebraWithOne( U ) );

    # Maintain useful information.
    UseIsomorphismRelation( U, S );
    UseSubsetRelation( U, S );
    if samecoeffs and HasDimension( U ) then
      SetDimension( S, Dimension( U ) );
    fi;

    # Return the subalgebra.
    return S;
    end );


#############################################################################
##
#M  AsSubalgebraWithOne(<A>, <U>) . . . view algebra as subalgebra of another
##
InstallMethod( AsSubalgebraWithOne,
    "for two algebras",
    IsIdenticalObj,
    [ IsAlgebra, IsAlgebra ],
    function( A, U )
    local S;
    if not IsSubset( A, U ) or One( U ) = fail then
      return fail;
    fi;

    if LeftActingDomain( A ) <> LeftActingDomain( U ) then
      U:= AsAlgebraWithOne( LeftActingDomain( A ), U );
    fi;

    # Construct and return the subalgebra.
    S:= SubalgebraWithOneNC( A, GeneratorsOfAlgebra( U ) );

    UseIsomorphismRelation( U, S );
    UseSubsetRelation( U, S );

    return S;
    end );

InstallMethod( AsSubalgebraWithOne,
    "for an algebra and a algebra-with-one",
    IsIdenticalObj,
    [ IsAlgebra, IsAlgebraWithOne ],
    function( A, U )
    local S;
    if not IsSubset( A, U ) then
      return fail;
    fi;

    if LeftActingDomain( A ) <> LeftActingDomain( U ) then
      U:= AsAlgebraWithOne( LeftActingDomain( A ), U );
    fi;

    # Construct and return the subalgebra.
    S:= SubalgebraWithOneNC( A, GeneratorsOfAlgebraWithOne( U ) );

    UseIsomorphismRelation( U, S );
    UseSubsetRelation( U, S );

    return S;
    end );


#############################################################################
##
#F  CentralizerInFiniteDimensionalAlgebra( <A>, <S>, <issubset> )
##
##  Let $(b_1, \ldots, b_n)$ be a basis of <A>, and $(s_1, \ldots, s_m)$
##  be a basis of <S>, with $s_j = \sum_{l=1}^n v_{jl} b_l$.
##  The structure constants of <A> are $c_{ijk}$ with
##  $b_i b_j = \sum_{k=1}^n c_{ijk} b_k$.
##  Then we compute a basis of the solution space of the system
##  $\sum_{i=1}^n a_i \sum_{l=1}^n v_{jl} ( c_{ilk} - c_{lik} )$ for
##  $1 \leq j \leq m$ and $1 \leq k \leq n$.
##
##  (left null space of an $n \times (nm)$ matrix)
##
##  If the multiplication in <A> is known to be anticommutative this is used.
##  (Note that the case of commutative multiplication is handled in a more
##  general way.)
##
#T We will have problems with this approach if centralizers of algebras over
#T non-commutative coefficients domains are considered.
#T In this case, <S> might stand for generators w.r.t. a larger coefficients
#T domain that also must be centralized ...
##
InstallGlobalFunction( CentralizerInFiniteDimensionalAlgebra,
    function( A, S, issubset )

    local B,           # basis of `A'
          T,           # structure constants table w. r. to `B'
          n,           # dimension of `A'
          m,           # length of `S'
          M,           # matrix of the equation system
          v,           # coefficients of basis vectors of `S' w.r. to `B'
          zerovector,  # initialize one row of `M'
          row,         # one row of `M'
          i, j, k, l,  # loop variables
          cil, cli,    #
          offset,
          pos;

    # Handle the case that `S' may be not contained in `A'.
    # (If `A' knows a basis and `S' is a subset of `A' then
    # there are methods that return `A' itself as closure.)
    if not issubset then
      M:= ClosureAlgebra( A, S );
      return Intersection2( A,
                 CentralizerInFiniteDimensionalAlgebra( M, S, true ) );
    fi;

    # Now `S' is known to be contained in `A'.
    B:= Basis( A );
    T:= StructureConstantsTable( B );
    n:= Dimension( A );
    m:= Length( S );
    M:= [];
    v:= List( S, x -> Coefficients( B, x ) );

    zerovector:= [ 1 .. n*m ] * Zero( LeftActingDomain( A ) );

    if HasIsAnticommutative( A ) and IsAnticommutative( A ) then

      # Column $(j-1)*n + k$ contains in row $i$ the value
      # $\sum_{l=1}^n v_{jl} c_{ilk}$

      for i in [ 1 .. n ] do
        row:= ShallowCopy( zerovector );
        for l in [ 1 .. n ] do
          cil:= T[i][l];
          for j in [ 1 .. m ] do
            offset:= (j-1)*n;
            for k in [ 1 .. Length( cil[1] ) ] do
              pos:= cil[1][k] + offset;
              row[ pos ]:= row[ pos ] + v[j][l] * cil[2][k];
            od;
          od;
        od;
        Add( M, row );
      od;

    else

      # Column $(j-1)*n + k$ contains in row $i$ the value
      # $\sum_{l=1}^n v_{jl} ( c_{ilk} - c_{lik} )$

      for i in [ 1 .. n ] do
        row:= ShallowCopy( zerovector );
        for l in [ 1 .. n ] do
          cil:= T[i][l];
          cli:= T[l][i];
          for j in [ 1 .. m ] do
            offset:= (j-1)*n;
            for k in [ 1 .. Length( cil[1] ) ] do
              pos:= cil[1][k] + offset;
              row[ pos ]:= row[ pos ] + v[j][l] * cil[2][k];
            od;
            for k in [ 1 .. Length( cli[1] ) ] do
              pos:= cli[1][k] + offset;
              row[ pos ]:= row[ pos ] - v[j][l] * cli[2][k];
            od;
          od;
        od;
        Add( M, row );
      od;

    fi;

    # Solve the equation system.
    M:= NullspaceMat( M );

    # Construct the generators from the coefficient vectors.
    M:= List( M, x -> LinearCombination( B, x ) );

    # Return the subalgebra.
    if IsFLMLORWithOne( A ) then
      return SubalgebraWithOneNC( A, M, "basis" );
    else
      return SubalgebraNC( A, M, "basis" );
    fi;
end );


#############################################################################
##
#M  CentralizerOp( <A>, <S> ) . . . . . cent. of a vector space in an algebra
##
InstallMethod( CentralizerOp,
    "for a finite dimensional algebra and a vector space with parent",
    IsIdenticalObj,
    [ IsAlgebra, IsVectorSpace and HasParent ],
    function( A, S )
    if    not IsIdenticalObj( A, Parent( S ) )
       or not IsFiniteDimensional( A ) then
      TryNextMethod();
    fi;
    return CentralizerInFiniteDimensionalAlgebra( A,
               BasisVectors( Basis( S ) ),
               true );
    end );


#############################################################################
##
#M  CentralizerOp( <A>, <S> ) . . .  cent. of an algebra in an assoc. algebra
##
InstallMethod( CentralizerOp,
    "for a fin. dim. assoc. algebra and an algebra with parent",
    IsIdenticalObj,
    [ IsAlgebra and IsAssociative, IsAlgebra and HasParent ],
    function( A, S )
    if    not IsIdenticalObj( A, Parent( S ) )
       or not IsFiniteDimensional( A ) then
      TryNextMethod();
    fi;
    return CentralizerInFiniteDimensionalAlgebra( A,
               GeneratorsOfAlgebra( S ),
               true );
    end );


#############################################################################
##
#M  CentralizerOp( <A>, <S> ) . . . . . cent. of a vector space in an algebra
##
InstallMethod( CentralizerOp,
    "for a finite dimensional algebra and a vector space",
    IsIdenticalObj,
    [ IsAlgebra, IsVectorSpace ],
    function( A, S )
    if not IsFiniteDimensional( A ) then
      TryNextMethod();
    fi;
    return CentralizerInFiniteDimensionalAlgebra( A,
               BasisVectors( Basis( S ) ),
               false );
    end );


#############################################################################
##
#M  CentralizerOp( <A>, <S> ) . . .  cent. of an algebra in an assoc. algebra
##
InstallMethod( CentralizerOp,
    "for a fin. dim. assoc. algebra and an algebra",
    IsIdenticalObj,
    [ IsAlgebra and IsAssociative, IsAlgebra ],
    function( A, S )
    if not IsFiniteDimensional( A ) then
      TryNextMethod();
    fi;
    return CentralizerInFiniteDimensionalAlgebra( A,
               GeneratorsOfAlgebra( S ),
               false );
    end );


#############################################################################
##
#M  CentralizerOp( <A>, <elm> ) . . . . . . cent. of an element in an algebra
##
InstallMethod( CentralizerOp,
    "for an algebra and an element",
    IsCollsElms,
    [ IsAlgebra, IsObject ],
    function( A, elm )
    return Centralizer( A, Algebra( LeftActingDomain( A ), [ elm ] ) );
    end );


#############################################################################
##
#F  CentreFromSCTable( <T> )
##
##  Given a structure constants table <T> w.r.t. a basis $B$, say,
##  `CentreFromSCTable' returns a list of $B$-coefficients vectors
##  of a basis for the centre of an $F$-algebra with s.c. table <T>,
##  where $F$ is assumed to be commutative.
##
##  We have to solve the system
##  $\sum_{i=1}^n a_i ( c_{ijk} - c_{jik} ) = 0$ for $1 \leq j, k \leq n$.
##
BindGlobal( "CentreFromSCTable", function( T )
    local n,       # the dimension
          M,       # matrix of the equation system
          i, j, k, # loop variables
          row,     # one row in `M'
          offset,  # offset between entry in `T' and column in `M'
          pos,     # nonzero positions in $c_{ij}$
          val;     # loop over structure constants in $c_{ij}$

    n:= Length( T ) - 2;
    M:= NullMat( n, n*n, T[ Length( T ) ] );
    for i in [ 1 .. n ] do
      row:= M[i];
      for j in [ 1 .. n ] do

        offset:= (j-1)*n;
        row{ offset + T[i][j][1] }:= T[i][j][2];
        pos:= T[j][i][1];
        val:= T[j][i][2];
        for k in [ 1 .. Length( pos ) ] do
          row[ offset + pos[k] ]:= row[ offset + pos[k] ] - val[k];
#T cheaper!
        od;

      od;
    od;

    # Solve the equation system.
    return NullspaceMat( M );
    end );


#############################################################################
##
#M  Centre( <A> )
##
InstallMethod( Centre,
    "for a finite dimensional FLMLOR",
    [ IsFLMLOR ],
    function( A )
    local   C,        # centre of `A', result
            B,        # a basis of `A'
            M;        # matrix of the equation system

    if not IsFiniteDimensional( A ) then
      TryNextMethod();
    fi;

    # If necessary convert `A' to a FLMLOR over a commutative ring.
    if not IsCommutative( LeftActingDomain( A ) ) then
      A:= AsFLMLOR( Centre( LeftActingDomain( A ) ), A );
    fi;

    # Solve the equation system.
    # If a s.c. table is already known, we use it since this allows us to
    # avoid multiplications.
    # Only if no s.c. table is known and if the algebra is known to be
    # associative and if the number of algebra generators is less than the
    # number of basis vectors, we do not force the computation of a
    # s.c. table.
    # (Note that for associative algebras,
    # we have to check $x a = a x$ only for algebra generators $a$,
    # not for all vectors of a basis.)
    B:= Basis( A );
    if    HasStructureConstantsTable( B )
       or not ( HasIsAssociative( A ) and IsAssociative( A ) ) then
      M:= CentreFromSCTable( StructureConstantsTable( B ) );
    else
      if HasGeneratorsOfAlgebraWithOne( A ) then
        M:= GeneratorsOfAlgebraWithOne( A );
      else
        M:= GeneratorsOfAlgebra( A );
      fi;
      if Dimension( A ) <= 2 * Length( M ) then
        M:= CentreFromSCTable( StructureConstantsTable( B ) );
      else
        M:= List( BasisVectors( B ), bi -> Concatenation( List( M,
                             a -> Coefficients( B, bi * a - a * bi ) ) ) );
        M:= NullspaceMat( M );
      fi;
    fi;

    # Construct the generators from the coefficient vectors.
    M:= List( M, x -> LinearCombination( B, x ) );

    # Construct the centre.
    if IsFLMLORWithOne( A ) then
      C:= SubalgebraWithOneNC( A, M, "basis" );
    else
      C:= SubalgebraNC( A, M, "basis" );
    fi;
    Assert( 1, IsAbelian( C ) );
    SetIsAbelian( C, true );

    # Return the centre.
    return C;
    end );


#############################################################################
##
#F  MutableBasisOfProductSpace( <U>, <V> )
##
##  Once we have the basis vectors for the product space,
##  we only have to decide whether the result of `ProductSpace' is an ideal,
##  an algebra, or just a vector space.
##  This decision is left to the methods of `ProductSpace',
##  the computation of basis vectors is done by `MutableBasisOfProductSpace'.
##
BindGlobal( "MutableBasisOfProductSpace", function( U, V )
    local inter, # intersection of left acting domains
          u, v,  # loop over the bases
          MB;    # mutable basis of the commutator subspace, result

    if LeftActingDomain( U ) = LeftActingDomain( V ) then
      inter:= LeftActingDomain( U );
    else
      inter:= Intersection2( LeftActingDomain( U ), LeftActingDomain( V ) );
      U:= AsVectorSpace( inter, U );
      V:= AsVectorSpace( inter, V );
    fi;

    MB:= MutableBasis( inter, [], Zero( U ) );
    V:= BasisVectors( Basis( V ) );
    for u in BasisVectors( Basis( U ) ) do
      for v in V do
        CloseMutableBasis( MB, u * v );
      od;
    od;

    # Return the result.
    return [ MB, inter ];
end );


#############################################################################
##
#M  ProductSpace( <U>, <V> )  . . . . . . . . . . . for two free left modules
##
InstallMethod( ProductSpace,
    "for two free left modules",
    IsIdenticalObj,
    [ IsFreeLeftModule, IsFreeLeftModule ],
    function( U, V )
    local MB, vectors, C;

    # Compute basis vectors.
    MB:= MutableBasisOfProductSpace( U, V );

    vectors:= BasisVectors( MB[1] );
    if IsEmpty( vectors ) then
      return TrivialSubspace( U );
    fi;

    # Create the appropriate domain.
    if HasParent( U ) and HasParent( V )
                      and IsIdenticalObj( Parent( U ), Parent( V ) ) then
      C:= SubmoduleNC( Parent( U ), vectors, "basis" );
    else
      C:= FreeLeftModule( MB[2], vectors, "basis" );
    fi;

    # Insert the basis.
    SetBasis( C, ImmutableBasis( MB[1] ) );

    # Return the result.
    return C;
    end );


#############################################################################
##
#M  ProductSpace( <U>, <V> )  . . . . . . . . . . . . . . .  for two algebras
##
##  If $<U> = <V>$ is known to be an algebra then the product space is also
##  an algebra, moreover it is an ideal in <U>.
##  If <U> and <V> are known to be ideals in an algebra $A$
##  then the product space is known to be an algebra and an ideal in $A$.
##
InstallMethod( ProductSpace,
    "for two algebras",
    IsIdenticalObj,
    [ IsAlgebra, IsAlgebra ],
    function( U, V )
    local P, MB, C;

    # Look for the ideal relation that allows one to construct an ideal.
    if IsIdenticalObj( U, V ) then
      P:= U;
    elif HasParent( V ) and IsIdenticalObj( Parent( V ), U )
                        and HasIsTwoSidedIdealInParent( V )
                        and IsTwoSidedIdealInParent( V ) then
      P:= U;
    elif HasParent( U ) and IsIdenticalObj( Parent( U ), V )
                        and HasIsTwoSidedIdealInParent( U )
                        and IsTwoSidedIdealInParent( U ) then
      P:= V;
    else
      TryNextMethod();
    fi;

    # Compute basis vectors for the result.
    MB:= MutableBasisOfProductSpace( U, V )[1];

    # The result is an ideal in `U'.
    C:= SubalgebraNC( P, BasisVectors( MB ), "basis" );

    SetIsTwoSidedIdealInParent( C, true );
    SetBasis( C, ImmutableBasis( MB, C ) );

    # Return the result.
    return C;
    end );


#############################################################################
##
#M  ProductSpace( <U>, <V> )  . . . . . . . . for two ideals with same parent
##
InstallMethod( ProductSpace,
    "for two ideals with same parent",
    IsIdenticalObj,
    [ IsAlgebra and HasParent and IsTwoSidedIdealInParent,
      IsAlgebra and HasParent and IsTwoSidedIdealInParent ],
    function( U, V )
    local MB, C;

    if not IsIdenticalObj( Parent( U ), Parent( V ) ) then
      TryNextMethod();
    fi;

    # The result is an ideal in the parent of `U'.
    MB:= MutableBasisOfProductSpace( U, V )[1];
    C:= SubalgebraNC( Parent( U ), BasisVectors( MB ), "basis" );

    SetIsTwoSidedIdealInParent( C, true );
    SetBasis( C, ImmutableBasis( MB, C ) );

    # Return the result.
    return C;
    end );


#############################################################################
##
#M  RadicalOfAlgebra( <A> ) . . . . . . . . radical of an associative algebra
##
##  `RadicalOfAlgebra' computes the radical (maximal nilpotent ideal)
##  of an associative algebra <A> by first constructing a faithful
##  matrix representation.
##  (Note that there is a special implementation for associative matrix
##  algebras.)
##
##  If <A> contains an identity element then the adjoint representation is
##  already faithful.
##
##  Otherwise we add an identity element (Dorroh extension).
##  More precisely we consider the space $B = \{ (x,t) | x\in A, t\in F }$.
##  We let <A> act on this space via $a (x,t) = (ax+ta,0)$.
##  Then it is easily seen that this representation is faithful.
##
InstallMethod( RadicalOfAlgebra,
    "for an associative algebra",
    [ IsAlgebra ],
    function( A )
    local bb,    # list of matrices representing the basis elements of <A>
          n,     # dimension of <A>
          BA,    # basis of `A'
          bv,    # basis vectors of `BA'
          F,     # left acting domain of `A'
          M,     # (n+1) x (n+1) matrix
          i,j,   # loop variables
          col,   # a column of `M'
          B,     # the representation of <A>
          R,     # the radical of `B'
          bas,   # a basis of `B' (corresponding to the basis of <A>)
          rad;   # a basis of the radical of <A>

    # Make sure that the algebra is associative and not a matrix algebra.
    if not IsAssociative( A ) then
      TryNextMethod();
    fi;

    n:= Dimension( A );
    BA:= Basis( A );
    bv:= BasisVectors( BA );
    F:= LeftActingDomain( A );

    if One( A ) <> fail then

      bb:= List( bv, x -> AdjointMatrix( BA, x ) );

    else

      bb:= [];

      for i in [1..n] do
        M:=[];
        for j in [1..n] do
          col:= Coefficients( BA, bv[i] * bv[j] );
          col[n+1]:= Zero( F );
          Add( M, col );
        od;
        col:= [ 1 .. n+1 ] * Zero( F );
        col[i]:= One( F );
        Add( M, col );
        Add( bb, M );
      od;

    fi;

    # Compute the radical of the isomorphic matrix algebra.
    B:= Algebra( F, bb, "basis" );
    R:= RadicalOfAlgebra( B );

    # Transfer the radical back to the original algebra.
    bas:= BasisNC( B, bb );
    rad:= List( BasisVectors( Basis( R ) ),
                x -> LinearCombination( BA, Coefficients( bas, x ) ) );

    return SubalgebraNC( A, rad, "basis" );
    end );


#############################################################################
##
#M  IsTrivial( <A> )  . . . . . . . . . . . . . . . . . . . . .  for a FLMLOR
##
InstallMethod( IsTrivial,
    "for a FLMLOR",
    [ IsFLMLOR ],
    A -> ForAll( GeneratorsOfLeftOperatorRing( A ), IsZero ) );


#############################################################################
##
#M  IsTrivial( <A> )  . . . . . . . . . . . . . . . . . for a FLMLOR-with-one
##
##  A FLMLOR-with-one is trivial if and only if its identity is equal to its
##  zero.
##
InstallMethod( IsTrivial,
    "for a FLMLOR-with-one",
    [ IsFLMLORWithOne ],
    A -> IsZero( One( A ) ) );


#############################################################################
##
#M  GeneratorsOfLeftModule( <A> )
##
##  We assume that it is possible to construct a basis for the algebra <A>.
##  If <A> is finite dimensional and if we know algebra generators
##  then the process of successive closure under the action of <A> on itself
##  yields this.
##
InstallMethod( GeneratorsOfLeftModule,
    "for a FLMLOR",
    [ IsFLMLOR ],
    A -> BasisVectors( Basis( A ) ) );


#############################################################################
##
#M  Basis( <A> )  . . . . . . . . . . . .  basis from FLMLOR gens. for FLMLOR
##
InstallMethod( Basis,
    "for a FLMLOR",
    [ IsFLMLOR ],
    function( A )

    # If generators as left module are known
    # we do not need to multiply at all.
    if HasGeneratorsOfLeftModule( A ) then
      TryNextMethod();
    fi;

    return ImmutableBasis( MutableBasisOfNonassociativeAlgebra(
             LeftActingDomain( A ),
             GeneratorsOfLeftOperatorRing( A ),
             Zero( A ),
             infinity ), A );
    end );


#############################################################################
##
#M  Basis( <A> )  . . . . . .  basis from FLMLOR gens. for associative FLMLOR
##
InstallMethod( Basis,
    "for an associative FLMLOR",
    [ IsFLMLOR and IsAssociative ],
        function( A )
    local  mb;    

    # If generators as left module are known
    # we do not need to multiply at all.
    if HasGeneratorsOfLeftModule( A ) then
      TryNextMethod();
    fi;
  
    mb:= MutableBasisOfClosureUnderAction(
             LeftActingDomain( A ),
             GeneratorsOfLeftOperatorRing( A ),
             "left",
             GeneratorsOfLeftOperatorRing( A ),
             \*,
             Zero( A ),
             infinity );
    return ImmutableBasis( mb, A );
    
    end );


#############################################################################
##
#M  Basis( <A> )   .  basis from FLMLOR gens. for associative FLMLOR-with-one
##
InstallMethod( Basis,
    "for an associative FLMLOR-with-one",
    [ IsFLMLORWithOne and IsAssociative ],
    function( A )

    # If generators as left module are known
    # we do not need to multiply at all.
    if HasGeneratorsOfLeftModule( A ) then
      TryNextMethod();
    fi;

    return ImmutableBasis( MutableBasisOfClosureUnderAction(
             LeftActingDomain( A ),
             GeneratorsOfLeftOperatorRing( A ),
             "left",
             [ One( A ) ],
             \*,
             Zero( A ),
             infinity ), A );
    end );


#############################################################################
##
#M  Basis( <A> )  . . . . . . . . . . basis from FLMLOR gens. for Lie algebra
##
##  In a Lie algebra, every word (with brackets) in terms of algebra
##  generators is a linear combination of left-normed words;
##  this means that it is sufficient to multiply with generators from one
##  side.
##
InstallMethod( Basis,
    "for a Lie algebra",
    [ IsLieAlgebra ],
    function( A )
    local  mb;    

    # If generators as left module are known
    # we do not need to multiply at all.
    if HasGeneratorsOfLeftModule( A ) then
      TryNextMethod();
    fi;
    
    mb:= MutableBasisOfClosureUnderAction(
             LeftActingDomain( A ),
             GeneratorsOfLeftOperatorRing( A ),
             "left",
             GeneratorsOfLeftOperatorRing( A ),       
             \*,
             Zero( A ),
             infinity );
    return ImmutableBasis( mb, A );    
    end );


#############################################################################
##
#M  PowerSubalgebraSeries( <A> )
##
InstallOtherMethod( PowerSubalgebraSeries,
    "for an algebra",
    [ IsAlgebra ],
    function ( A )
    local   S,          # power subalgebra series of <A>, result
            D;          # power subalgebras

    # Compute the series by repeated calling of `ProductSpace'.
    S := [ A ];
    D := ProductSpace( A, A );
    while D <> S[ Length(S) ]  do
      Add( S, D );
      D:= ProductSpace( D, D );
    od;

    # Return the series when it becomes stable.
    return S;
    end );


#############################################################################
##
#M  IsNilpotentElement( <L>, <x> )  . . . . . . for an algebra and an element
##
##  <x> is nilpotent if its adjoint matrix $A$ (i.e. the linear map coming
##  from left multiplication by <x>) is nilpotent.
##  To check this, we only need to check whether $A^n$ (or a smaller power)
##  is zero, where $n$ denotes the dimension of <L>.
##
InstallMethod( IsNilpotentElement,
    "for an algebra, and an element",
    IsCollsElms,
    [ IsAlgebra, IsRingElement ],
    function( L, x )
    local B,     # a basis of `L'
          A,     # adjoint matrix of `x w.r. to `B'
          n,     # dimension of `L'
          i,     # loop variable
          zero;  # zero coefficient

    B := Basis( L );
    A := AdjointMatrix( B, x );
    n := Dimension( L );
    i := 1;
    zero:= Zero( A[1][1] );

    if ForAll( A, x -> n < PositionNot( x, zero ) ) then
      return true;
    fi;

    while i < n do
      i:= 2 * i;
      A:= A * A;
      if ForAll( A, x -> n < PositionNot( x, zero ) ) then
        return true;
      fi;
    od;

    return false;
    end );


#############################################################################
##
#M  GeneratorsOfLeftOperatorRing( <A> ) . . . . . . . . for a FLMLOR-with-one
##
InstallMethod( GeneratorsOfLeftOperatorRing,
    "for a FLMLOR-with-one",
    [ IsFLMLORWithOne ],
    A -> Concatenation( [ One( A ) ],
                        GeneratorsOfLeftOperatorRingWithOne( A ) ) );


#############################################################################
##
#M  GeneratorsOfLeftOperatorRing( <A> ) . . . .  for FLMLOR with module gens.
##
InstallMethod( GeneratorsOfLeftOperatorRing,
    "for a FLMLOR with known left module generators",
    [ IsFLMLOR and HasGeneratorsOfLeftModule ],
    GeneratorsOfLeftModule );


#############################################################################
##
#M  GeneratorsOfLeftOperatorRingWithOne( <A> ) . for FLMLOR with module gens.
##
InstallMethod( GeneratorsOfLeftOperatorRingWithOne,
    "for a FLMLOR-with-one with known left module generators",
    [ IsFLMLORWithOne and HasGeneratorsOfLeftModule ],
    GeneratorsOfLeftModule );


#############################################################################
##
#M  DirectSumOfAlgebras( <A1>, <A2> )
##
##  Construct a s.c. algebra.
##  (There are special methods for the sum of appropriate matrix algebras.)
##
#T embeddings/projections should be provided!
##
InstallOtherMethod( DirectSumOfAlgebras,
    "for two algebras",
    [ IsAlgebra, IsAlgebra ],
    function( A1, A2 )
    local n,     # The dimension of the resulting algebra.
          i,j,   # Loop variables.
          T,     # The table of structure constants of the direct sum.
          scT,   #
          n1,    # The dimension of A1.
          n2,    # The dimension of A2.
          ll,    # A list of structure constants.
          L,     # result.
          sym,   # if both products are (anti)symmetric, then the result
                 # will have the same property.
          R1,R2, # Root systems of A1,A2.
          f1,f2, # Embeddings of A1,A2 in L.
          R,     # Root system of L.
          RV,    # List of various things.
          r,
          pos;   # List of positions.

    if LeftActingDomain( A1 ) <> LeftActingDomain( A2 ) then
      Error( "<A1> and <A2> must be written over the same field" );
    fi;

    n1:= Dimension( A1 );
    n2:= Dimension( A2 );
    n:= n1+n2;
    T:= [];

    # Initialize the s.c. table.
    T:= EmptySCTable( n, Zero( LeftActingDomain( A1 ) ) );

    # Enter the structure constants for the first algebra.
    scT:= StructureConstantsTable( Basis( A1 ) );
    sym:= scT[n1+1];
    T{ [ 1 .. n1 ] }{ [ 1 .. n1 ] }:= scT{ [ 1 .. n1 ] };

    scT:= StructureConstantsTable( Basis( A2 ) );

    for i in [1..n2] do
      for j in [1..n2] do
        ll:= ShallowCopy( scT[i][j] );
        ll[1]:= ll[1] + n1;
        T[n1+i][n1+j]:= ll;
      od;
    od;


    # Set the (anti)symmetric flag
    if scT[n2 + 1] = sym  then
        T[n + 1] := sym;
    fi;
    if Characteristic( LeftActingDomain( A1 ) ) = 2 and sym in [ 1, -1 ]
        and scT[n2 + 1] in [ 1, -1 ]  then
        T[n + 1] := 1;
    fi;


    L:= AlgebraByStructureConstants( LeftActingDomain( A1 ), T );

    # Maintain useful information.
    if     HasIsLieAlgebra( A1 ) and HasIsLieAlgebra( A2 )
       and IsLieAlgebra( A1 ) and IsLieAlgebra( A2 ) then
       SetIsLieAlgebra( L, true );
       if HasRootSystem( A1 ) and HasRootSystem( A2 ) then
           # We can easily compute the root system of `L'.
           R1:= RootSystem( A1 ); R2:= RootSystem( A2 );
           f1:= LeftModuleGeneralMappingByImages( A1, L, CanonicalBasis(A1),
                        CanonicalBasis(L){[1..Dimension(A1)]} );
           f2:= LeftModuleGeneralMappingByImages( A2, L, CanonicalBasis(A2),
                        CanonicalBasis(L){[Dimension(A1)+1..Dimension(L)]} );
           R:= Objectify( NewType( NewFamily( "RootSystemFam", IsObject ),
                       IsAttributeStoringRep and IsRootSystemFromLieAlgebra ),
                       rec() );
           RV:= List( PositiveRootVectors( R1 ), x -> Image( f1, x ) );
           Append( RV, List( PositiveRootVectors( R2 ),
                                               x -> Image( f2, x ) ) );
           SetPositiveRootVectors( R, RV );
           RV:= List( NegativeRootVectors( R1 ), x -> Image( f1, x ) );
           Append( RV, List( NegativeRootVectors( R2 ),
                                               x -> Image( f2, x ) ) );
           SetNegativeRootVectors( R, RV );
           RV:= List( PositiveRoots( R1 ), ShallowCopy );
           for i in [1..Length(RV)] do
               Append( RV[i], ListWithIdenticalEntries(
                       Length( CartanMatrix( R2 ) ),
                       Zero( LeftActingDomain( A2 ) ) ) );
           od;
           for i in PositiveRoots( R2 ) do
               r:= ListWithIdenticalEntries( Length( CartanMatrix( R1 ) ),
                           Zero( LeftActingDomain( A1 ) ) );
               Append( r, i );
               Add( RV, r );
           od;
           SetPositiveRoots( R, RV );
           RV:= List( NegativeRoots( R1 ), ShallowCopy );
           for i in [1..Length(RV)] do
               Append( RV[i], ListWithIdenticalEntries(
                       Length( CartanMatrix( R2 ) ),
                       Zero( LeftActingDomain( A2 ) ) ) );
           od;
           for i in NegativeRoots( R2 ) do
               r:= ListWithIdenticalEntries( Length( CartanMatrix( R1 ) ),
                           Zero( LeftActingDomain( A1 ) ) );
               Append( r, i );
               Add( RV, r );
           od;
           SetNegativeRoots( R, RV );
           pos:= List( SimpleSystem( R1 ), x -> Position(
                         PositiveRoots(R1), x ) );
           RV:= PositiveRoots( R ){pos};
           pos:= List( SimpleSystem( R2 ), x -> Position(
                       PositiveRoots(R2), x ) + Length( PositiveRoots(R1) ));
           Append( RV, PositiveRoots( R ){pos} );
           SetSimpleSystem( R, RV );
           RV:= [ ];
           for i in [1,2,3] do
               RV[i]:= List( CanonicalGenerators( R1 )[i],
                             x -> Image( f1, x ) );
               Append( RV[i], List( CanonicalGenerators( R2 )[i],
                              x -> Image( f2, x ) ) );
           od;
           SetCanonicalGenerators( R, RV );
           SetCartanMatrix( R, DirectSumMat( CartanMatrix( R1 ),
                   CartanMatrix( R2 ) ) );
           SetUnderlyingLieAlgebra( R, L );
           SetRootSystem( L, R );
           if HasChevalleyBasis( A1 ) and HasChevalleyBasis( A2 ) then
               RV:= [ ];
               for i in [1,2,3] do
                   RV[i]:= List( ChevalleyBasis( A1 )[i],
                                 x -> Image( f1, x ) );
                   Append( RV[i], List( ChevalleyBasis( A2 )[i],
                           x -> Image( f2, x ) ) );
               od;
               SetChevalleyBasis( L, RV );
           fi;
       fi;
       
       if HasSemiSimpleType( A1 ) and HasSemiSimpleType( A2 ) then
           SetSemiSimpleType( L, Concatenation( SemiSimpleType(A1)," ",
                   SemiSimpleType( A2 ) ) );
       fi;

       if HasIsRestrictedLieAlgebra ( A1 ) and HasIsRestrictedLieAlgebra ( A2 )
         and IsRestrictedLieAlgebra ( A1 ) and IsRestrictedLieAlgebra ( A2 )
       then
          SetIsRestrictedLieAlgebra( L, true );
          if HasPthPowerImages( Basis ( A1 ) )
             and HasPthPowerImages( Basis ( A2 ) ) then
             if not IsBound (f1) then
                f1:= LeftModuleGeneralMappingByImages( A1, L,
                     CanonicalBasis(A1),
                     CanonicalBasis(L){[1..Dimension(A1)]} );
                f2:= LeftModuleGeneralMappingByImages( A2, L,
                     CanonicalBasis(A2),
                     CanonicalBasis(L){[Dimension(A1)+1..Dimension(L)]});
             fi;

             SetPthPowerImages( Basis ( L ),
               Concatenation ( List (PthPowerImages( Basis ( A1 ) ), x->x^f1),
                          List (PthPowerImages( Basis ( A2 ) ), x->x^f2) ) );

          fi;
       fi;

    fi;
    if     HasIsAssociative( A1 ) and HasIsAssociative( A2 )
       and IsAssociative( A1 ) and IsAssociative( A2 ) then
      SetIsAssociative( L, true );
    fi;

    # Return the result.
    return L;
    end );


#############################################################################
##
#M  DirectSumOfAlgebras( <list> ) . . . . . . .  for a dense list of algebras
##
InstallMethod( DirectSumOfAlgebras,
    "for list of algebras",
    [ IsDenseList ],
    function( list )
    local R, A, i;

    if IsEmpty( list ) then
      Error( "<list> must be nonempty" );
    fi;

    R:= LeftActingDomain( list[1] );
    for A in list do
      if not IsFLMLOR( A ) or LeftActingDomain( A ) <> R then
        Error( "all entries in <list> must be FLMLORs over <R>" );
      fi;
    od;

    A:= list[1];
    for i in [ 2 .. Length( list ) ] do
      A:= DirectSumOfAlgebras( A, list[i] );
    od;

    return A;
    end );


#############################################################################
##
#O  IsCentral( <A>, <U> )  . . . . . . . .  test if <U> is centralized by <A>
##
##  Check whether every basis vector of <A> commutes with every basis vector
##  of the subset <U>.
##
##  For associative algebras, we have to check $u a = a u$ only for algebra
##  generators $a$ and $u$ of $A$ and $U$, respectively,
##  not for all vectors of a basis.
##
InstallMethod( IsCentral,
    "for two FLMLORs",
    IsIdenticalObj,
    [ IsFLMLOR, IsFLMLOR ],
    IsCentralFromGenerators( GeneratorsOfLeftModule,
                             GeneratorsOfLeftModule ) );

InstallMethod( IsCentral,
    "for two associative FLMLORs",
    IsIdenticalObj,
    [ IsFLMLOR and IsAssociative, IsFLMLOR and IsAssociative ],
    IsCentralFromGenerators( GeneratorsOfAlgebra,
                             GeneratorsOfAlgebra ) );

InstallMethod( IsCentral,
    "for two associative FLMLORs-with-one",
    IsIdenticalObj,
    [ IsFLMLORWithOne and IsAssociative,
      IsFLMLORWithOne and IsAssociative ],
    IsCentralFromGenerators( GeneratorsOfAlgebraWithOne,
                             GeneratorsOfAlgebraWithOne ) );


#############################################################################
##
#F  FreeAlgebraConstructor( <name>, <magma> )
##
##  is used for a uniform treatment of free (associative) algebras(-with-one)
##
BindGlobal( "FreeAlgebraConstructor", function( name, magma )
    return function( arg )
    local   R,          # coefficients ring
            names,      # names of the algebra generators
            M,          # magma
            A,          # algebra
            F,          # family
            i;

    # Check the argument list.
    if Length( arg ) = 0 or not IsRing( arg[1] ) then
      Error( "first argument must be a ring" );
    fi;

    R:= arg[1];

    # Construct names of generators.
    if   Length( arg ) = 2 and IsInt( arg[2] ) then
      names:= List( [ 1 .. arg[2] ],
                    i -> Concatenation( "x.", String(i) ) );
      MakeImmutable( names );
    elif     Length( arg ) = 2
         and IsList( arg[2] )
         and ForAll( arg[2], IsString ) then
      names:= arg[2];
    elif Length( arg ) = 3 and IsInt( arg[2] ) and IsString( arg[3] ) then
      names:= List( [ 1 .. arg[2] ],
                    x -> Concatenation( arg[3], ".", String(x) ) );
      MakeImmutable( names );
    elif ForAll( arg{ [ 2 .. Length( arg ) ] }, IsString ) then
      names:= arg{ [ 2 .. Length( arg ) ] };
    else
      Error( "usage: ", name, "( <R>, <rank> )\n",
                 "or ", name, "( <R>, <name1>, ... )" );
    fi;
    
    M := magma( names );
    
    # Construct the algebra as free magma algebra of a free magma over `R'.
    A := FreeMagmaRing( R, M );

    # Store the names.
    F := ElementsFamily( FamilyObj( A ) );
    F!.names:= names;
    
    # Install grading
    if HasOne(M) then i := 0; else i := 1; fi;
    SetGrading( A, rec(min_degree := i,
                                     max_degree := infinity,
                                     source := Integers,
                                     hom_components := function(degree)
        local i, d, B, x, y;
        if HasOne(M) then
            B := [[One(M)],GeneratorsOfMagmaWithOne(M)];
        else
            B := [[],GeneratorsOfMagma(M)];
        fi;
        for d in [2..degree] do
            Add(B,[]);
            if IsAssociative(M) then
                for x in B[2] do for y in B[d] do
                    Add(B[d+1],x*y);
                od; od;
            else
                for i in [2..d] do
                    for x in B[i] do for y in B[d+2-i] do
                        Add(B[d+1],x*y);
                    od; od;
                od;
            fi;
        od;
        x := Zero(R);
        y := [One(R)];
        return VectorSpace(R, List(B[degree+1],
                       p->ElementOfMagmaRing( F, x, y, [p] )));
    end));

    # Return the result.
    return A;
    end;
end );


#############################################################################
##
#F  FreeAlgebra( <R>, <rank> ) . . . . . . . . . . free algebra of given rank
#F  FreeAlgebra( <R>, <rank>, <name> )
#F  FreeAlgebra( <R>, <name1>, <name2>, ... )
##
InstallGlobalFunction( FreeAlgebra,
    FreeAlgebraConstructor( "FreeAlgebra", FreeMagma ) );


#############################################################################
##
#F  FreeAlgebraWithOne( <R>, <rank> ) . . free algebra-with-one of given rank
#F  FreeAlgebraWithOne( <R>, <rank>, <name> )
#F  FreeAlgebraWithOne( <R>, <name1>, <name2>, ... )
##
InstallGlobalFunction( FreeAlgebraWithOne,
    FreeAlgebraConstructor( "FreeAlgebraWithOne", FreeMagmaWithOne ) );


#############################################################################
##
#F  FreeAssociativeAlgebra( <R>, <rank> )
#F  FreeAssociativeAlgebra( <R>, <rank>, <name> )
#F  FreeAssociativeAlgebra( <R>, <name1>, <name2>, ... )
##
InstallGlobalFunction( FreeAssociativeAlgebra,
    FreeAlgebraConstructor( "FreeAssociativeAlgebra", FreeSemigroup ) );


#############################################################################
##
#F  FreeAssociativeAlgebraWithOne( <R>, <rank> )
#F  FreeAssociativeAlgebraWithOne( <R>, <rank>, <name> )
#F  FreeAssociativeAlgebraWithOne( <R>, <name1>, <name2>, ... )
##
InstallGlobalFunction( FreeAssociativeAlgebraWithOne,
    FreeAlgebraConstructor( "FreeAssociativeAlgebraWithOne", FreeMonoid ) );


#############################################################################
##
#M  \.( <F>, <n> )  . . . . . . . . .  access to generators of a free algebra
##
InstallAccessToGenerators( IsMagmaRingModuloRelations,
                           "magma ring containing the whole family",
                           GeneratorsOfAlgebra );

InstallAccessToGenerators( IsMagmaRingModuloRelations and IsRingWithOne,
                           "magma ring-with-one containing the whole family",
                           GeneratorsOfAlgebraWithOne );


#############################################################################
##
#M  CentralIdempotentsOfAlgebra( <A> )
##
##   Let A be an associative algebra with one. We construct a maximal
##   system of orthogonal primitive idemoptents in the centre of A.
##   First we let B be the centre of A and Q the
##   the semisimple commutative associative algebra A/Rad(A).
##   We calculate a complete set of orthogonal idempotents in `Q'
##   and then lift them to A.
##   The orthogonal idempotents in `Q' correspond to the decomposition
##   of `Q' as a direct sum of simple ideals. Now `ideals' will contain
##   a list of ideals of `Q' such that the direct sum of these equals
##   `Q'. The variable `ids' will contain the idempotents corresponding
##   to the ideals in `ids'.
##   The algorithms has two parts: one for small fields (of size less than
##   `2*Dimension( Q )', and one for big fields.
##   If the field is big, then using a Las Vegas algorithm we find a
##   splitting element (this is an element that generates `Q'). By
##   factoring the minimal polynomial of such element we can find a
##   complete set of orthogonal idempotents in one step.
##   However, if the field is small splitting elements might not exist.
##   In this case we use decomposable elements (of which the minimum
##   polynomial factors into two (or more) relatively prime factors.
##   Then using the same procedure as for splitting elements we can
##   find some idempotents. But in this case the corresponding ideals
##   might split further. So we have to find decomposable elements in
##   these and so on.
##   Decomposable elements are found as follows: first we calculate
##   the subalgebra of all elements x such that x^q=x
##   (where `q=Size( F )').
##   This subalgebra is a number of copies of the ground field. So any
##   element independent from 1 of this subalgebra will have a minimum
##   polynomial that splits completely. On the other hand, if 1 is the
##   only basis vector of this subalgebra than the original algebra was
##   simple.
##   For a more elaborate description we refer to "W. Eberly and M.
##   Giesbrecht, Efficient Decomposition of Associative Algebras,
##   Proceedings of ISSAC 1996."
##

InstallMethod( CentralIdempotentsOfAlgebra,
    "for an associative algebra",
    [ IsAlgebra ],
    function( A )
    local F,B,Rad,Q,bQ,ids,ideals,id,i,j,k,l,set,cf,e,vv,sp,x,f,q,sol,
          eq,facs,hlist,c,p,g,gcd,bb,E,ei,ni,hom,qq;


      if One( A ) = fail then TryNextMethod(); fi;

      F:=LeftActingDomain(A);
      B:=Centre(A);

      Rad:= RadicalOfAlgebra( B );
      hom:= NaturalHomomorphismByIdeal( B, Rad );
      Q:= ImagesSource( hom );
      bQ:= BasisVectors( Basis( Q ) );
      ids:= [ One( Q ) ];
      ideals:= [ Q ];

      # The variable `k' will point to the first element of `ideals' that
      # still has to be decomposed.

      k:=1;

      if Size(F) > 2*Dimension( Q )^2 then
        set:= [ 0 .. 2*Dimension(Q)^2 ]*One( F );
      else
        set:= [ ];
      fi;

      repeat

        if Length( set ) > 1 then

          # We are in the case of a big field.

          repeat

            # We try to find an element of `Q' that generates it.
            # If we take the coefficients of such an element randomly
            # from a set of `2*Dimension(Q)^2' elements,
            # then this element generates `Q' with probability > 1/2

            bQ:= BasisVectors( Basis( ideals[k] ) );
            cf:= List( [ 1 .. Length(bQ) ], x -> Random( set ) );
            e:= LinearCombination( bQ, cf );

            # Now we calculate the minimum polynomial of `e'.

            vv:= [ MultiplicativeNeutralElement( ideals[k] ) ];
            sp:= MutableBasis( F, vv );
            x:= ShallowCopy( e );

            while not IsContainedInSpan( sp, x ) do
              Add( vv, x );
              CloseMutableBasis( sp, x );
              x:= x*e;
            od;
            sp:= UnderlyingLeftModule( ImmutableBasis( sp ) );
            cf:= ShallowCopy(
                   - Coefficients( BasisNC( sp, vv ), x )
                 );
            Add( cf, One( F ) );
            f:= ElementsFamily( FamilyObj( F ) );
            f:= LaurentPolynomialByCoefficients( f, cf, 0 );

          until DegreeOfLaurentPolynomial( f ) = Dimension( Q );

        else

          # Here the field is small.

          q:= Size( F );

        # `sol' will be a basis of the subalgebra of the k-th ideal
        # consisting of all elements x such that x^q=x.
        # If the length of this list is 1,
        # then the ideal is simple and we proceed to the next one. If all
        # ideals are simple then we quit the loop.

          sol:= [ ];
          while Length( sol ) < 2 and k <= Length( ideals ) do
            bQ:= BasisVectors( Basis( ideals[k] ) );
            eq:= [ ];
            for i in [1..Dimension( ideals[k] )] do
              Add( eq, Coefficients( Basis( ideals[k] ), bQ[i]^q-bQ[i] ) );
            od;
            sol:= List( NullspaceMat( eq ),
                        x -> LinearCombination( bQ, x ) );
            if Length(sol) = 1 then k:=k+1; fi;
          od;

          if k>Length(ideals) then break; fi;

          vv:= [ MultiplicativeNeutralElement( ideals[k] ) ];
          sp:= MutableBasis( F, vv );

          e:= sol[1];
          if IsContainedInSpan( sp, e ) then e:=sol[2]; fi;

        # We calculate the minimum polynomial of `e'.

          x:= ShallowCopy( e );
          while not IsContainedInSpan( sp, x ) do
            Add( vv, x );
            CloseMutableBasis( sp, x );
            x:= x*e;
          od;
          sp:= UnderlyingLeftModule( ImmutableBasis( sp ) );
          cf:=  ShallowCopy(
                  - Coefficients( BasisNC( sp, vv ), x )
                );
          Add( cf, One( F ) );

          f:= ElementsFamily( FamilyObj( F ) );
          f:= LaurentPolynomialByCoefficients( f, cf, 0 );

        fi;

        facs:= Factors( PolynomialRing( F ), f );

      # Now we find elements h1,...,hs such that `hi = 1 mod facs[i]' and
      # `hi = 0 mod facs[j]' if `i<>j'.
      # This is possible due to the Chinese remainder theorem.

        hlist:= [ ];
        for i in [1..Length( facs )] do
          cf:= List( [ 1..Length( facs ) ], x -> Zero( F ) );
          cf[i]:= One(F);
          j:= 1;
          c:= cf[1];
          p:= facs[1];
          while j < Length(facs) do
            j:= j + 1;
            g:= GcdRepresentation( p, facs[j] );
            gcd:= g[1]*p+g[2]*facs[j];
            qq:= g[1]*(cf[j]-c)/gcd;
            if qq<>0*qq then
              c:= p*EuclideanRemainder( qq, facs[j] )
                              + c;
            fi;
            p:= p*facs[j] / gcd;
          od;

          Add( hlist, EuclideanRemainder( c*facs[i]^0 , p ) );

        od;

      # Now a set of orthogonal idempotents is given by `hi(e)'.
      # We evaluate `hi(e)' in a rather strange way; this in order to make
      # sure that the one is the one of `ideals[ k ]' ('e^0' will be the
      # one of the big algebra `Q').

        id:= List( hlist, x -> Value( x, e,
                        MultiplicativeNeutralElement( ideals[k] ) ) );

        if Length(set) = 0 then

        # We are in the case of a small field;
        # so we append the new idempotents and ideals
        # (and erase the old ones). (If `E' is an idempotent,
        # then the corresponding ideal is given by `E*Q*E'.)

          Append(ids,id);

          for l in [1..Length(id)] do
            bb:=List(BasisVectors(Basis(ideals[k])),x->id[l]*x*id[l]);
            Add(ideals,Subalgebra(Q,bb));
          od;

          ideals:=Filtered(ideals,x->x<>ideals[k]);
          ids:=Filtered(ids,x->x<>ids[k]);
        else

        # Here the field is big so we found the complete list of idempotents
        # in one step.

          ids:= id;
          k:=Length(ideals)+1;
        fi;

        while k<=Length(ideals) and Dimension( ideals[k] ) = 1 do k:=k+1; od;

      until k>Length(ideals);

      id:= List( ids, e -> PreImagesRepresentative( hom, e ) );

      # Now we lift the idempotents to the big algebra `A'. The
      # first idempotent is lifted as follows:
      # We have that `id[1]^2-id[1]' is an element of `Rad'.
      # We construct the sequences e_{i+1} = e_i + n_i - 2e_in_i,
      # and n_{i+1}=e_{i+1}^2-e_{i+1}, starting with e_0=id[1].
      # It can be proved by induction that e_q is an idempotent in `A'
      # because n_0^{2^q}=0.
      # Now `E' will be the sum of all idempotents lifted so far.
      # Then the next lifted idempotent is obtained by setting
      # `ei:=id[i]-E*id[i]-id[i]*E+E*id[i]*E;'
      # and lifting as above. It can be proved that in this manner we
      # get an orthogonal system of primitive idempotents.

      E:= Zero( F )*id[1];

      for i in [1..Length(id)] do
        ei:= id[i]-E*id[i]-id[i]*E+E*id[i]*E;
        q:= 0;
        while 2^q <= Dimension( Rad ) do
          q:= q+1;
        od;
        ni:= ei*ei-ei;
        for j in [1..q] do
          ei:= ei+ni-2*ei*ni;
          ni:= ei*ei-ei;
        od;
        id[i]:= ei;
        E:= E+ei;
      od;

      return AsSSortedList(id);
end );


##############################################################################
##
#M  IsSimpleAlgebra( <A> )  . . . . . . . . . . . .for an associative algebra
##
##  A test whether <A> is simple.
##
InstallMethod( IsSimpleAlgebra,
    "for an associative algebra",
    [ IsAlgebra ],
    function( A )

    if not IsAssociative( A ) then
      TryNextMethod();
    fi;

    if Dimension( RadicalOfAlgebra( A ) ) <> 0 then
      return false;
    else
      return Length( CentralIdempotentsOfAlgebra( A ) ) = 1;
    fi;
    end );


###############################################################################
##
#M  LeviMalcevDecomposition( <L> )
##
##  A Levi-Malcev subalgebra of `L' is a semisimple subalgebra complementary to
##  the radical `R'. We find a Levi-Malcev subalgebra of `L' by first
##  computing a complementary subspace to `R'. This subspace is a Levi-Malcev
##  subalgebra modulo `R'. Then we change the basis vectors such that they
##  form a basis of a Levi-Malcev subalgebra modulo the second term of the
##  derived series of `R' after that we consider the third term of the
##  derived series, and so on.
##
InstallMethod( LeviMalcevDecomposition,
    "for an associative or a Lie algebra",
    [ IsAlgebra ],
    function( L )
    local R,             # The solvable radical of `L'.
          s,             # The dimension of the Levi subalgebra.
          F,             # coefficients domain of `L'
          bas,bb,        # Lists of basiselements.
          sp,            # A vector space.
          subalg,        # Boolean.
          a,i,j,k,l,m,   # Loop variables.
          x,             # Element of `L'.
          ser,           # The derived series of `R'.
          p,             # The length of the derived series.
          Rbas,          # A special basis of `R'.
          levi,          # Basis of a Levi complement.
          T,             # Structure constants table of `L', w.r.t. a
                         # particular basis.
          cf,cf1,cf2,    # Coefficient vectors.
          klist,         # List of integers.
          comp,          # List of basis vectors of a complement.
          dim,           # The length of `comp'.
          B,             # A basis.
          cij,           # Entry of the table of structure constants.
          eqs,           # Matrix of equation set.
          rl,            # Right hand side of the equation system.
          eqno,          # Number of the equation.
          sol,           # Solution set to the equations.
          r,offset;      # Integers.

    if IsLieAlgebra( L ) then
      R:= LieSolvableRadical( L );
      offset:= 1;
    elif IsAssociative( L ) then
      R:= RadicalOfAlgebra( L );
      offset:=0;
    fi;

    if Dimension( R ) = 0 then
      return [ L, R ];
    elif Dimension( R ) = Dimension( L ) then
      return [ TrivialSubalgebra( L ), R ];
    fi;

    s:= Dimension( L ) - Dimension( R );

    # `bb' will be a basis of a complement to `R' in `L'.

    bas:= ShallowCopy( BasisVectors( Basis( R ) ) );
    F:= LeftActingDomain( L );
    sp:= MutableBasis( F, bas );
    bb:= [ ];
    for k in BasisVectors( Basis( L ) ) do
      if Length( bb ) = s then
        break;
      elif not IsContainedInSpan( sp, k ) then
        Add( bb, k );
        CloseMutableBasis( sp, k );
      fi;
    od;

    sp:= MutableBasis( F, bb );
    subalg:= true;
    for i in [1..Length(bb)] do
      for j in [offset*i+1..Length(bb)] do
        if not IsContainedInSpan( sp, bb[i]*bb[j] ) then
          subalg:= false;
          break;
        fi;
      od;
    od;
    if subalg then
      Info( InfoAlgebra, 1,
            "LeviDecomposition: subalgebra test successful" );
      return [ SubalgebraNC( L, bb, "basis" ), R ];
    fi;

    ser:= PowerSubalgebraSeries( R );

    # We now calculate a basis of `R' such that the first k1 elements
    # form a basis of the last nonzero term of the derived series `ser',
    # the first k2 ( k2>k1 ) elements form a basis of the next to last
    # element of the derived series, and so on.

    p:= Length( ser );
    i:= p-1;
    Rbas:= ShallowCopy( BasisVectors( Basis( ser[p-1] ) ) );
    sp:= MutableBasis( F, Rbas );
    while Length(Rbas) < Dimension(R) do
      if Length(Rbas) = Dimension(ser[i]) then
        i:= i-1;
        k:= 1;
      else
        x:= BasisVectors( Basis( ser[i] ) )[k];
        if not IsContainedInSpan( sp, x ) then
          Add( Rbas, x );
          CloseMutableBasis( sp, x );
        fi;
        k:= k+1;
      fi;
    od;

    levi:= ShallowCopy( bb );
    Append(bb,Rbas);

    # So now `bb' is a list of basis vectors of `L' such that
    # the first elements form a basis of a complement to `R'
    # and the remaining elements are a basis for `R' of the form
    # described above.
    # We now calculate a structure constants table of `L' w.r.t. this basis.

    sp:= VectorSpace( F, bb );
    B:= BasisNC( sp, bb );
    T:= List([1..s],x->[]);
    for i in [1..s] do
      for j in [offset*i+1..s] do
        cf:= Coefficients( B, levi[i]*levi[j] ){[1..s]};
        klist:= Filtered([1..s],i->cf[i]<>0);
        cf:= Filtered(cf,x->x<>0);
        T[i][j]:= [klist,cf];
      od;
    od;

    # Now `levi' is a Levi-Malcev subalgebra modulo `R'.
    # The loop modifies this modulo statement.
    # After the first round `levi' will be a Levi-Malcev subalgebra modulo
    # the second element of the derived series.
    # After the second step `levi' will be a Levi-Malcev subalgebra modulo
    # the third element of the derived series.
    # And so on.

    for a in [1..p-1] do

      # `comp' will be a basis of the complement of the `a+1'-th element
      # of the derived series in the `a'-th element of the derived series.
      # `B' will be a basis of the `a'-th term of the derived series,
      # such that the basis elements of the complement come first.
      # So if we have an element v of the `a'-th term of the derived series,
      # then by taking the coefficients w.r.t. `B', we can easily find
      # the part that belongs to `comp'.
      # The equations we have are vector equations in the space `comp',
      # i.e., in the quotient of two elements of the derived series.
      # But we do not want to work with this quotient directly.

      comp:= Rbas{ [ Dimension(ser[a+1])+1 .. Dimension(ser[a]) ] };
      dim:= Length(comp);
      bb:= ShallowCopy( comp );
      for i in [1..Dimension(ser[a+1])] do
        Add(bb,Rbas[i]);
      od;
      sp:= VectorSpace( F, bb );
      B:= BasisNC( sp, bb );

      cf:= List( comp, x -> Coefficients( B, x ){[1..dim]} );

      eqs:= NullMat( s*dim, dim*s*(s-offset)/(offset+1), F );
      rl:= [];
      for i in [1..s] do
        for j in [offset*i+1..s] do
          cij:= T[i][j];
          for k in [1..dim] do

            cf1:= Coefficients(B,levi[i]*comp[k]){[1..dim]};
            cf2:= Coefficients(B,comp[k]*levi[j]){[1..dim]};

            for l in [1..dim] do
              if IsAssociative( L ) then
                eqno:= (i-1)*s*dim+(j-1)*dim+l;
              else
                eqno:= (i-1)*(2*s-i)*dim/2+(j-i-1)*dim+l;
              fi;

              eqs[(j-1)*dim+k][eqno]:= eqs[(j-1)*dim+k][eqno]+cf1[l];
              eqs[(i-1)*dim+k][eqno]:= eqs[(i-1)*dim+k][eqno]+cf2[l];

              for m in [1..Length(cij[1])] do
                r:=cij[1][m];
                if r <= s then
                  eqs[(r-1)*dim+k][eqno]:= eqs[(r-1)*dim+k][eqno]-
                                           cij[2][m]*cf[k][l];
                fi;
              od;
            od;
          od;

          x:= Zero(L);
          for m in [1..Length(cij[1])] do
            if cij[1][m] <= s then
              x:= x+cij[2][m]*levi[cij[1][m]];
            fi;
          od;
          x:= x-levi[i]*levi[j];
          Append(rl,Coefficients(B,x){[1..dim]});

        od;
      od;

      sol:= SolutionMat( eqs, rl );

      if sol = fail then return sol; fi;

      for i in [1..s] do
        for j in [1..dim] do
          levi[i]:=levi[i]+sol[(i-1)*dim+j]*comp[j];
        od;
      od;
    od;

    return [ SubalgebraNC( L, levi, "basis" ), R ];
    end );


############################################################################
##
#M  DirectSumDecomposition( <A> )   ........direct sum decomposition of <A>
##
InstallMethod( DirectSumDecomposition,
        "for semisimple associative algebras",
        [ IsAlgebra and IsAssociative ],
    function( A )
    local R;

    R:= RadicalOfAlgebra( A );
    if Dimension( R ) > 0 then
        TryNextMethod();
    fi;
    return List( CentralIdempotentsOfAlgebra( A ), x -> Ideal( A, [ x ] ) );
end );


#############################################################################
##
#E