Ziele Untersuchung
mit Columbo Integrität von
Datenbanken Interaktion und
Portierbarkeit Ergonomie der
Schnittstellen

Angebot Produkte Projekt Beratung

Mittel Analytik Modellierung Sprachen Algebra Logik Hardware Denken Kreativität

Zusammenhänge Gesellschaft Wirtschaft Branche Firma

Benutzer


products/Sources/formale Sprachen/GAP/lib/ (GAP Algebra Version 4.15.1^©) Datei vom 18.9.2025 mit Größe 32 kB

Impressum alglie.gi Sprache: Shell

#############################################################################
##
## This file is part of GAP, a system for computational discrete algebra.
## This file's authors include Thomas Breuer, and Willem de Graaf.
##
## Copyright of GAP belongs to its developers, whose names are too numerous
## to list here. Please refer to the COPYRIGHT file for details.
##
## SPDX-License-Identifier: GPL-2.0-or-later
##
## This file contains methods for Lie algebras.
##

#############################################################################
##
#M LieUpperCentralSeries( <L> ) . . . . . . . . . . for a Lie algebra
##
InstallMethod( LieUpperCentralSeries,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )

 local S, # upper central series of <L>, result
 C, # Lie centre
 hom; # homomorphisms of <L> to `<L>/<C>'

 S := [ TrivialSubalgebra( L ) ];
 C := LieCentre( L );
 while C <> Last(S) do

 # Replace `L' by `L / C', compute its centre, and get the preimage
 # under the natural homomorphism.
 Add( S, C );
 hom:= NaturalHomomorphismByIdeal( L, C );
 C:= PreImages( hom, LieCentre( Range( hom ) ) );
#T we would like to get ideals!
#T is it possible to teach the hom. that the preimage of an ideal is an ideal?

 od;

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

#############################################################################
##
#M LieLowerCentralSeries( <L> ) . . . . . . . . . . for a Lie algebra
##
InstallMethod( LieLowerCentralSeries,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )

 local S, # lower central series of <L>, result
 C; # commutator subalgebras

 # Compute the series by repeated calling of `ProductSpace'.
 S := [ L ];
 C := LieDerivedSubalgebra( L );
 while C <> Last(S) do
 Add( S, C );
 C:= ProductSpace( L, C );
 od;

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

#############################################################################
##
#M LieDerivedSubalgebra( <L> )
##
## is the (Lie) derived subalgebra of the Lie algebra <L>.
## This is the ideal/algebra/subspace (equivalent in this case)
## generated/spanned by all products $uv$
## where $u$ and $v$ range over a basis of <L>.
##
InstallMethod( LieDerivedSubalgebra,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 L -> ProductSpace( L, L ) );

#############################################################################
##
#M LieDerivedSeries( <L> )
##
InstallMethod( LieDerivedSeries,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function ( L )

 local S, # (Lie) derived series of <L>, result
 D; # (Lie) derived subalgebras

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

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

#############################################################################
##
#M IsLieSolvable( <L> ) . . . . . . . . . . . . . . . for a Lie algebra
##
InstallMethod( IsLieSolvable,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )

 local D;

 D:= LieDerivedSeries( L );
 return Dimension( Last(D) ) = 0;
 end );

InstallTrueMethod( IsLieSolvable, IsLieNilpotent );

#############################################################################
##
#M IsLieNilpotent( <L> ) . . . . . . . . . . . . . . . for a Lie algebra
##
InstallMethod( IsLieNilpotent,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )

 local D;

 D:= LieLowerCentralSeries( L );
 return Dimension( Last(D) ) = 0;
 end );

InstallTrueMethod( IsLieNilpotent, IsLieAbelian );

#############################################################################
##
#M IsLieAbelian( <L> ) . . . . . . . . . . . . . . for a Lie algebra
##
## It is of course sufficient to check products of algebra generators,
## no basis and structure constants of <L> are needed.
## But if we have already a structure constants table we use it.
##
InstallMethod( IsLieAbelian,
 "for a Lie algebra with known basis",
 true,
 [ IsAlgebra and IsLieAlgebra and HasBasis ], 0,
 function( L )

 local B, # basis of `L'
 T, # structure constants table w.r.t. `B'
 i, # loop variable
 j; # loop variable

 B:= Basis( L );
 if not HasStructureConstantsTable( B ) then
 TryNextMethod();
 fi;

 T:= StructureConstantsTable( B );
 for i in T{ [ 1 .. Length( T ) - 2 ] } do
 for j in i do
 if not IsEmpty( j[1] ) then
 return false;
 fi;
 od;
 od;
 return true;
 end );

InstallMethod( IsLieAbelian,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )

 local i, # loop variable
 j, # loop variable
 zero, # zero of `L'
 gens; # algebra generators of `L'

 zero:= Zero( L );
 gens:= GeneratorsOfAlgebra( L );
 for i in [ 1 .. Length( gens ) ] do
 for j in [ 1 .. i-1 ] do
 if gens[i] * gens[j] <> zero then
 return false;
 fi;
 od;
 od;

 # The algebra multiplication is trivial, and the algebra does
 # not know about a basis.
 # Here we know at least that the algebra generators are space
 # generators.
 if not HasGeneratorsOfLeftModule( L ) then
 SetGeneratorsOfLeftModule( L, gens );
 fi;

 # Return the result.
 return true;
 end );

InstallTrueMethod( IsLieAbelian, IsAlgebra and IsZeroMultiplicationRing );

##############################################################################
##
#M LieCentre( <L> ) . . . . . . . . . . . . . . . . . . . for a Lie algebra
##
## We solve the system
## $\sum_{i=1}^n a_i c_{ijk} = 0$ for $1 \leq j, k \leq n$
## (instead of $\sum_{i=1}^n a_i ( c_{ijk} - c_{jik} ) = 0$).
##
## Additionally we know that the centre of a Lie algebra is an ideal.
##
InstallMethod( LieCentre,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( A )

 local R, # left acting domain of `A'
 C, # Lie centre of `A', result
 B, # a basis of `A'
 T, # structure constants table w.r. to `B'
 n, # dimension of `A'
 M, # matrix of the equation system
 zerovector, #
 i, j, # loop over ...
 row; # one row of `M'

 R:= LeftActingDomain( A );

 if Characteristic( R ) <> 2 and HasCentre( A ) then

 C:= Centre( A );
#T change it to an ideal!

 else

 # Catch the trivial case.
 n:= Dimension( A );
 if n = 0 then
 return A;
 fi;

 # Construct the equation system.
 B:= Basis( A );
 T:= StructureConstantsTable( B );
 M:= [];
 zerovector:= [ 1 .. n*n ] * Zero( R );
 for i in [ 1 .. n ] do
 row:= ShallowCopy( zerovector );
 for j in [ 1 .. n ] do
 if IsBound( T[i][j] ) then
 row{ (j-1)*n + T[i][j][1] }:= T[i][j][2];
 fi;
 od;
 M[i]:= row;
 od;

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

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

 # Construct the Lie centre.
 C:= IdealNC( A, M, "basis" );

 fi;

 # Return the Lie centre.
 return C;
 end );

##############################################################################
##
#M LieCentralizer( <A>, <S> ) . . . . . for a Lie algebra and a vector space
##
## 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}^m 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}^m v_{jl} c_{ilk} = 0$ for
## $1 \leq j \leq m$ and $1 \leq k \leq n$.
##
## (left null space of an $n \times (nm)$ matrix)
##
InstallMethod( LieCentralizer,
 "for an abelian Lie algebra and a vector space",
 IsIdenticalObj,
 [ IsAlgebra and IsLieAlgebra and IsLieAbelian,
 IsVectorSpace ], 0,
 function( A, S )

 if IsSubset( A, S ) then
 return A;
 else
 TryNextMethod();
 fi;
 end );

InstallMethod( LieCentralizer,
 "for a Lie algebra and a vector space",
 IsIdenticalObj,
 [ IsAlgebra and IsLieAlgebra, IsVectorSpace ], 0,
 function( A, S )

 local R, # left acting domain of `A'
 B, # basis of `A'
 T, # structure constants table w. r. to `B'
 n, # dimension of `A'
 m, # dimension 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, #
 offset,
 vjl,
 pos;

 # catch trivial case
 if Dimension(S) = 0 then
 return A;
 fi;

 R:= LeftActingDomain( A );
 B:= Basis( A );
 T:= StructureConstantsTable( B );
 n:= Dimension( A );
 m:= Dimension( S );
 M:= [];
 v:= List( BasisVectors( Basis( S ) ),
 x -> Coefficients( B, x ) );

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

 # Column $(j-1)*n + k$ contains in row $i$ the value
 # $\sum_{l=1}^m 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;
 vjl := v[j][l];
 for k in [ 1 .. Length( cil[1] ) ] do
 pos:= cil[1][k] + offset;
 row[ pos ]:= row[ pos ] + vjl * cil[2][k];
 od;
 od;
 od;
 Add( M, row );
 od;

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

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

 # Return the subalgebra.

 return SubalgebraNC( A, M, "basis" );

 end );

##############################################################################
##
#M LieNormalizer( <L>, ) . . . . . . for a Lie algebra and a vector space
##
## If $(x_1, \ldots, x_n)$ is a basis of $L$ and $(u_1, \ldots, u_s)$ is
## a basis of $U$, then $x = \sum_{i=1}^n a_i x_i$ is an element of $N_L(U)$
## iff $[x,u_j] = \sum_{k=1}^s b_{j,k} u_k$ for $j = 1, \ldots, s$.
## This leads to a set of $n s$ equations for the $n + s^2$ unknowns $a_i$
## and $b_{jk}$.
## If $u_k= \sum_{l=1}^n v_{kl} x_l$, then these equations can be written as
## $\sum_{i=1}^n (\sum_{j=1}^n v_{lj} c_{ijk} a_i -
## \sum_{i=1}^s v_{ik} b_{li} = 0$,
## for $1 \leq k \leq n$ and $1 \leq j \leq s$,
## where the $c_{ilp}$ are the structure constants of $L$.
## From the solution we only need the "normalizer part" (i.e.,
## the $a_i$ part).
##
InstallMethod( LieNormalizer,
 "for a Lie algebra and a vector space",
 IsIdenticalObj,
 [ IsAlgebra and IsLieAlgebra, IsVectorSpace ], 0,
 function( L, U )

 local R, # left acting domain of `L'
 B, # a basis of `L'
 T, # the structure constants table of `L' w.r.t. `B'
 n, # the dimension of `L'
 s, # the dimension of `U'
 A, # the matrix of the equation system
 i, j, k, l, # loop variables
 v, # the coefficients of the basis of `U' wrt `B'
 cij,
 bas,
 b,
 pos;

 # catch trivial case
 if Dimension(U) = 0 then
 return L;
 fi;

 # We need not work if `U' knows to be an ideal in its parent `L'.
 if HasParent( U ) and IsIdenticalObj( L, Parent( U ) )
 and HasIsLeftIdealInParent( U ) and IsLeftIdealInParent( U ) then
 return L;
 fi;

 R:= LeftActingDomain( L );
 B:= Basis( L );
 T:= StructureConstantsTable( B );
 n:= Dimension( L );
 s:= Dimension( U );

 if s = 0 or n = 0 then
 return L;
 fi;

 v:= List( BasisVectors( Basis( U ) ),
 x -> Coefficients( B, x ) );

 # The equations.
 # First the normalizer part, \ldots

 A:= NullMat( n + s*s, n*s, R );
 for i in [ 1..n ] do
 for j in [ 1..n ] do
 cij:= T[i][j];
 for l in [ 1..s ] do
 for k in [ 1..Length( cij[1] ) ] do
 pos:= (l-1)*n+cij[1][k];
 A[i][pos]:= A[i][pos]+v[l][j]*cij[2][k];
 od;
 od;
 od;
 od;

 # \ldots and then the "superfluous" part.

 for k in [1..n] do
 for l in [1..s] do
 for i in [1..s] do
 A[ n+(l-1)*s+i ][ (l-1)*n+k ]:= -v[i][k];
 od;
 od;
 od;

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

 # Extract the `normalizer part' of the solution.
 l:= Length(b);
 bas:= NullMat( l, n, R );
 for i in [ 1..l ] do
 for j in [ 1..n ] do
 bas[i][j]:= b[i][j];
 od;
 od;

 # Construct the generators from the coefficients list.
 bas:= List( bas, x -> LinearCombination( B, x ) );

 # Return the subalgebra.
 return SubalgebraNC( L, bas, "basis" );
 end );

##############################################################################
##
#M KappaPerp( <L>, ) . . . . . . . . for a Lie algebra and a vector space
##
#T Should this better be `OrthogonalSpace( <F>, )' where <F> is a
#T bilinear form?
#T How to represent forms in GAP?
#T (Clearly the form must know about the space <L>.)
##
## If $(x_1,\ldots, x_n)$ is a basis of $L$ and $(u_1,\ldots, u_s)$ is a
## basis of $U$ such that $u_k = \sum_{j=1}^n v_{kj} x_j$ then an element
## $x = \sum_{i=1}^n a_i x_i$ is an element of $U^{\perp}$ iff the $a_i$
## satisfy the equations
## $\sum_{i=1}^n ( \sum_{j=1}^n v_{kj} \kappa(x_i,x_j) ) a_i = 0$ for
## $k = 1, \ldots, s$.
##
InstallMethod( KappaPerp,
 "for a Lie algebra and a vector space",
 IsIdenticalObj,
 [ IsAlgebra and IsLieAlgebra, IsVectorSpace ], 0,
 function( L, U )

 local R, # left acting domain of `L'
 B, # a basis of L
 kap, # the matrix of the Killing form w.r.t. `B'
 A, # the matrix of the equation system
 n, # the dimension of L
 s, # the dimension of U
 v, # coefficient list of the basis of U w.r.t. the basis of L
 i,j,k, # loop variables
 bas; # the basis of the solution space

 R:= LeftActingDomain( L );
 B:= Basis( L );
 n:= Dimension( L );
 s:= Dimension( U );

 if s = 0 or n = 0 then
 return L;
 fi;

 v:= List( BasisVectors( Basis( U ) ),
 x -> Coefficients( B, x ) );
 A:= NullMat( n, s, R );
 kap:= KillingMatrix( B );

 # Compute the equations that define the subspace.
 for k in [ 1..s ] do
 for i in [ 1..n ] do
 for j in [ 1..n ] do
 A[i][k]:= A[i][k] + v[k][j] * kap[i][j];
 od;
 od;
 od;

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

 # Extract the generators.
 bas:= List( bas, x -> LinearCombination( B, x ) );

 return SubspaceNC( L, bas, "basis" );

 end );

#############################################################################
##
#M AdjointMatrix( , <x> )
##
## If the basis vectors are $(b-1, b_2, \ldots, b_n)$, and
## $x = \sum_{i=1}^n x_i b_i$ then $b_j$ is mapped to
## $[ x, b_j ] = \sum_{i=1}^n x_i [ b_i b_j ]
## = \sum_{k=1}^n ( \sum_{i=1}^n x_i c_{ijk} ) b_k$,
## so the entry in the $k$-th row and the $j$-th column of the adjoint
## matrix is $\sum_{i=1}^n x_i c_{ijk}$.
##
## Note that $ad_x$ is a left multiplication, so also the action of the
## adjoint matrix is from the left (i.e., on column vectors).
##
InstallMethod( AdjointMatrix,
 "for a basis of a Lie algebra, and an element",
 IsCollsElms,
 [ IsBasis, IsRingElement ], 0,
 function( B, x )

 local n, # dimension of the algebra
 T, # structure constants table w.r. to `B'
 zerovector, # zero of the field
 M, # adjoint matrix, result
 j, i, l, # loop variables
 cij, # structure constants vector
 k, # one position in structure constants vector
 row; # one row of `M'

 x:= Coefficients( B, x );
 n:= Length( BasisVectors( B ) );
 T:= StructureConstantsTable( B );
 zerovector:= [ 1 .. n ] * Last(T);
 M:= [];
 for j in [ 1 .. n ] do
 row:= ShallowCopy( zerovector );
 for i in [ 1 .. n ] do
 cij:= T[i][j];
 for l in [ 1 .. Length( cij[1] ) ] do
 k:= cij[1][l];
 row[k]:= row[k] + x[i] * cij[2][l];
 od;
 od;
 M[j]:= row;
 od;

 return TransposedMat( M );
 end );

#T general function for arbitrary algebras? (right/left multiplication)
#T RegularRepresentation: right multiplication satisfies M_{xy} = M_x M_y
#T is just the negative of the adjoint ...

#############################################################################
##
#M RightDerivations( )
##
## Let $n$ be the dimension of $A$.
## We start with $n^2$ indeterminates $D = [ d_{i,j} ]_{i,j}$ which
## means that $D$ maps $b_i$ to $\sum_{j=1}^n d_{ij} b_j$.
##
## (Note that this is row convention.)
##
## This leads to the following linear equation system in the $d_{ij}$.
## $\sum_{k=1}^n ( c_{ijk} d_{km} - c_{kjm} d_{ik} - c_{ikm} d_{jk} ) = 0$
## for all $1 \leq i, j, m \leq n$.
## The solution of this system gives us a vector space basis of the
## algebra of derivations.
##
InstallMethod( RightDerivations,
 "method for a basis of an algebra",
 true,
 [ IsBasis ], 0,
 function( B )

 local T, # structure constants table w.r. to 'B'
 L, # underlying Lie algebra
 R, # left acting domain of 'L'
 n, # dimension of 'L'
 eqno,offset,
 A,
 i, j, k, m,
 M; # the Lie algebra of derivations

 if not IsAlgebra( UnderlyingLeftModule( B ) ) then
 Error( " must be a basis of an algebra" );
 fi;

 if IsLieAlgebra( UnderlyingLeftModule( B ) ) then
 offset:= 1;
 else
 offset:= 0;
 fi;

 T:= StructureConstantsTable( B );
 L:= UnderlyingLeftModule( B );
 R:= LeftActingDomain( L );
 n:= Dimension( L );

 if n = 0 then
 return NullAlgebra( R );
 fi;

 # The rows in the matrix of the equation system are indexed
 # by the $d_{ij}$; the $((i-1) n + j)$-th row belongs to $d_{ij}$.

 # Construct the equation system.
 if offset = 1 then
 A:= NullMat( n^2, (n-1)*n*n/2, R );
 else
 A:= NullMat( n^2, n^3, R );
 fi;
 eqno:= 0;
 for i in [ 1 .. n ] do
 for j in [ offset*i+1 .. n ] do
 for m in [ 1 .. n ] do
 eqno:= eqno+1;
 for k in [ 1 .. n ] do
 A[ (k-1)*n+m ][eqno]:= A[ (k-1)*n+m ][eqno] +
 SCTableEntry( T,i,j,k );
 A[ (i-1)*n+k ][eqno]:= A[ (i-1)*n+k ][eqno] -
 SCTableEntry( T,k,j,m );
 A[ (j-1)*n+k ][eqno]:= A[ (j-1)*n+k ][eqno] -
 SCTableEntry( T,i,k,m );
 od;
 od;
 od;
 od;

 # Solve the equation system.
 # Note that if `L' is a Lie algebra and $n = 1$ the matrix is empty.

 if n = 1 and offset = 1 then
 A:= [ [ One( R ) ] ];
 else
 A:= NullspaceMatDestructive( A );
 fi;

 # Construct the generating matrices from the vectors.
 A:= List( A, v -> List( [ 1 .. n ],
 i -> v{ [ (i-1)*n + 1 .. i*n ] } ) );

 # Construct the Lie algebra.
 if IsEmpty( A ) then
 M:= AlgebraByGenerators( R, [],
 LieObject( Immutable( NullMat( n, n, R ) ) ) );
 else
 A:= List( A, LieObject );
 M:= AlgebraByGenerators( R, A );
 UseBasis( M, A );
 fi;

 # Return the derivations.
 return M;
end );

#############################################################################
##
#M LeftDerivations( )
##
## Let $n$ be the dimension of $A$.
## We start with $n^2$ indeterminates $D = [ d_{i,j} ]_{i,j}$ which
## means that $D$ maps $b_i$ to $\sum_{j=1}^n d_{ji} b_j$.
##
## (Note that this is column convention.)
##
InstallMethod( LeftDerivations,
 "method for a basis of an algebra",
 true,
 [ IsBasis ], 0,
 function( B )

 local T, # structure constants table w.r. to 'B'
 L, # underlying Lie algebra
 R, # left acting domain of 'L'
 n, # dimension of 'L'
 eqno,offset,
 A,
 i, j, k, m,
 M; # the Lie algebra of derivations

 if not IsAlgebra( UnderlyingLeftModule( B ) ) then
 Error( " must be a basis of an algebra" );
 fi;

 if IsLieAlgebra( UnderlyingLeftModule( B ) ) then
 offset:= 1;
 else
 offset:= 0;
 fi;

 T:= StructureConstantsTable( B );
 L:= UnderlyingLeftModule( B );
 R:= LeftActingDomain( L );
 n:= Dimension( L );

 if n = 0 then
 return NullAlgebra( R );
 fi;

 # The rows in the matrix of the equation system are indexed
 # by the $d_{ij}$; the $((i-1) n + j)$-th row belongs to $d_{ij}$.

 # Construct the equation system.
 if offset = 1 then
 A:= NullMat( n^2, (n-1)*n*n/2, R );
 else
 A:= NullMat( n^2, n^3, R );
 fi;
 eqno:= 0;
 for i in [ 1 .. n ] do
 for j in [ offset*i+1 .. n ] do
 for m in [ 1 .. n ] do
 eqno:= eqno+1;
 for k in [ 1 .. n ] do
 A[ (m-1)*n+k ][eqno]:= A[ (m-1)*n+k ][eqno] +
 SCTableEntry( T,i,j,k );
 A[ (k-1)*n+i ][eqno]:= A[ (k-1)*n+i ][eqno] -
 SCTableEntry( T,k,j,m );
 A[ (k-1)*n+j ][eqno]:= A[ (k-1)*n+j ][eqno] -
 SCTableEntry( T,i,k,m );
 od;
 od;
 od;
 od;

 # Solve the equation system.
 # Note that if `L' is a Lie algebra and $n = 1$ the matrix is empty.

 if n = 1 and offset = 1 then
 A:= [ [ One( R ) ] ];
 else
 A:= NullspaceMatDestructive( A );
 fi;

 # Construct the generating matrices from the vectors.
 A:= List( A, v -> List( [ 1 .. n ],
 i -> v{ [ (i-1)*n + 1 .. i*n ] } ) );

 # Construct the Lie algebra.
 if IsEmpty( A ) then
 M:= AlgebraByGenerators( R, [],
 LieObject( Immutable( NullMat( n, n, R ) ) ) );
 else
 A:= List( A, LieObject );
 M:= AlgebraByGenerators( R, A );
 UseBasis( M, A );
 fi;

 # Return the derivations.
 return M;
end );

#############################################################################
##
#M KillingMatrix( )
##
## We have $\kappa_{i,j} = \sum_{k,l=1}^n c_{jkl} c_{ilk}$ if $c_{ijk}$
## are the structure constants w.r. to .
##
## (The matrix is symmetric, no matter whether the multiplication is
## (anti-)symmetric.)
##
InstallMethod( KillingMatrix,
 "for a basis of a Lie algebra",
 true,
 [ IsBasis ], 0,
 function( B )

 local T, # s.c. table w.r. to `B'
 L, # the underlying algebra
 R, # left acting domain of `L'
 kappa, # the matrix of the killing form, result
 n, # dimension of `L'
 zero, # the zero of `R'
 i, j, k, t, # loop variables
 row, # one row of `kappa'
 val, # one entry of `kappa'
 cjk; # `T[j][k]'

 T:= StructureConstantsTable( B );
 L:= UnderlyingLeftModule( B );
 R:= LeftActingDomain( L );
 n:= Dimension( L );
 kappa:= [];
 zero:= Zero( R );

 for i in [ 1 .. n ] do
 row:= [];
 for j in [ 1 .. i ] do

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

 od;
 kappa[i]:= row;

 od;

 # Return the result.
 return kappa;
 end );

##############################################################################
##
#M AdjointBasis( )
##
## The input is a basis of a (Lie) algebra $L$.
## This function returns a particular basis $C$ of the matrix space generated
## by $ad L$, namely a basis consisting of elements of the form $ad x_i$
## where $x_i$ is a basis element of .
## An extra component `indices' is added to this space.
## This is a list of integers such that `ad .basisVectors[ indices[i] ]'
## is the `i'-th basis vector of <C>, for i in [1..Length(indices)].
## (This list is added in order to be able to identify the basis element of
## with the property that its adjoint matrix is equal to a given basis
## vector of <C>.)
##
InstallMethod( AdjointBasis,
 "for a basis of a Lie algebra",
 true,
 [ IsBasis ], 0,
 function( B )

 local bb, # the basis vectors of `B'
 n, # the dimension of `B'
 F, # the field over which the algebra is defined
 adL, # a list of matrices that form a basis of adLsp
 adLsp, # the matrix space spanned by ad L
 inds, # the list of indices
 i, # loop variable
 adi, # the adjoint matrix of the i-th basis vector of `B'
 adLbas; # the basis of `adLsp' compatible with `adL'

 bb:= BasisVectors( B );
 n:= Length( bb );
 F:= LeftActingDomain( UnderlyingLeftModule( B ) );
 adL:= [];
 adLsp:= MutableBasis( F, [ NullMat(n,n,F) ] );
#T better declare the zero ?
 inds:= [];
 for i in [1..n] do
 adi:= AdjointMatrix( B, bb[i] );
 if not IsContainedInSpan( adLsp, adi ) then
 Add( adL, adi );
 Add( inds, i );
 CloseMutableBasis( adLsp, adi );
 fi;
 od;

 if adL = [ ] then
 adLbas:= Basis( VectorSpace( F, [ ], NullMat( n, n, F ) ) );
 else
 adLbas:= Basis( VectorSpace( F, adL ), adL );
 fi;

 SetIndicesOfAdjointBasis( adLbas, inds );

 return adLbas;
 end );

##############################################################################
##
#M IsRestrictedLieAlgebra( <L> ) . . . . . . . . . . . . . for a Lie algebra
##
## A Lie algebra <L> is defined to be {\em restricted} when it is defined
## over a field of characteristic $p \neq 0$, and for every basis element
## $x$ of <L> there exists $y\in <L>$ such that $(ad x)^p = ad y$
## (see Jacobson, p. 190).
##
InstallMethod( IsRestrictedLieAlgebra,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )

 local F, # the field over which L is defined
 B, # the basis of L
 p, # the characteristic of F
 adL, # basis for the (matrix) vector space generated by ad L
 v; # loop over basis vectors of adL

 F:= LeftActingDomain( L );
 p:= Characteristic( F );
 if p = 0 then
 return false;
 fi;

 B:= Basis( L );

 adL:= AdjointBasis( B );

 # Check if ad(L) is closed under the p-th power operation.
 for v in BasisVectors( adL ) do
 if not v^p in UnderlyingLeftModule( adL ) then
 return false;
 fi;
 od;

 return true;
 end );

#############################################################################
##
#F PowerSi( <F>, )
##
InstallGlobalFunction( PowerSi, function( F, i )

 local si, # a function of two arguments: seqs, a list of sequences,
 # and l, a list containing the two arguments of the
 # function s_i. The list seqs contains
 # all possible sequences of i-1 1's and p-2-i+1 2's
 # This function returns the value of s_i(l[1],l[2])
 j,k, # loop variables
 p, # the characteristic of F
 combs, # the list of all i-1 element subsets of {1,2,\ldots, p-2}
 # it serves to make the list seqs
 v, # a vector of 1's and 2's
 seqs; # the list of all sequences of 1's and 2's of length p-2,
 # with i-1 1's and p-2 2's serving as input for si
 # for example, the sequence [1,1,2] means the element
 # [[[[x,y],x],x],y] (the first element [x,y] is present
 # 1 1 2 in all terms of the sum s_i(x,y)).

 si:= function( seqs, l )

 local x,
 j,k,
 sum;

 for j in [1..Length(seqs)] do
 x:= l[1]*l[2];
 for k in seqs[j] do
 x:= x*l[k];
 od;
 if j=1 then
 sum:= x;
 else
 sum:= sum+x;
 fi;
 od;
 return ( i * One( F ) )^(-1)*sum;
 end;

 p:= Characteristic( F );
 combs:= Combinations( [1..p-2], i-1 );

 # Now all sequences of 1's and 2's of length p and containing i-1 1's
 # are constructed the 1's in the jth sequence are put at the positions
 # contained in combs[j].
 seqs:=[];
 for j in [1..Length( combs )] do
 v:= List( [1..p-2], x -> 2);
 for k in combs[j] do
 v[k]:= 1;
 od;
 Add( seqs, v );
 od;
 return arg -> si( seqs, arg );
end );

#############################################################################
##
#F PowerS( <L> )
##
InstallMethod( PowerS,
 "for a Lie algebra",
 true,
 [ IsLieAlgebra ], 0,
 function( L )

 local F, # the coefficients domain
 p; # the characteristic of `F'

 F:= LeftActingDomain( L );
 p:= Characteristic( F );
 return List( [ 1 .. p-1 ], i -> PowerSi( F , i ) );
 end );

##############################################################################
##
#F PthPowerImage( , <x> )
##
BindGlobal("PTHPOWERIMAGE_PPI_VEC", function(L,zero,p,bv,pmap,cf,x)
 local
 n, # the dimension of L
 s, # the list of s_i functions
 im, # the image of x under the p-th power map
 i,j; # loop variables

 n:= Dimension( L );
 s:= PowerS( L );
 im:= Zero( L );

 # First the sum of all $\alpha_i^p x_i^{[p]}$ is calculated.
 for i in [1..n] do
 im:= im + cf[i]^p * pmap[i];
 od;

 # To this the double sum of all
 # $s_i(\alpha_j x_j, \sum_{k=j+1}^n \alpha_k x_k)$
 # is added.
 for j in [1..n-1] do
 if cf[j] <> zero then
 x:= x - cf[j] * bv[j];
 for i in [1..p-1] do
 im:= im + s[i]( [cf[j]*bv[j],x] );
 od;
 fi;
 od;

 return im;
end);

InstallMethod( PthPowerImage,
 "for a basis of an algebra, and a ring element",
 IsCollsElms,
 [ IsBasis, IsRingElement ], 0,
 function( B, x )

 local L, # the Lie algebra of which B is a basis
 F, # the coefficients domain of `L'
 n, # the dimension of L
 p, # the characteristic of the ground field
 s, # the list of s_i functions
 pmap, # the list containing x_i^{[p]}
 cf, # the coefficients of x wrt the basis of L
 im, # the image of x under the p-th power map
 i,j, # loop variables
 zero, # zero of `F'
 bv, # basis vectors of `B'
 adx, # adjoint matrix of x
 adL; # a basis of the matrix space ad L

 L:= UnderlyingLeftModule( B );
 if not IsLieAlgebra( L ) then
 TryNextMethod();
 fi;

 F:= LeftActingDomain( L );
 p:= Characteristic( F );

 if Dimension( LieCentre( L ) ) = 0 then

 # We calculate the inverse image $ad^{-1} ((ad x)^p)$.
 adx:= AdjointMatrix( B, x );
 adL:= AdjointBasis( B );
 adx:= adx^p;
 cf:= Coefficients( adL, adx );
 return LinearCombination( B, cf );

 else
 return PTHPOWERIMAGE_PPI_VEC(L,Zero(F),p,BasisVectors(B),PthPowerImages(B),Coefficients(B,x),x);
 fi;
 end );

InstallMethod( PthPowerImage, "for an element of a restricted Lie algebra",
 [ IsJacobianElement ], # weaker filter, we maybe only discovered later
 # that the algebra is restricted
 function(x)
 local fam;
 fam := FamilyObj(x);
 if not IsBound(fam!.pMapping) then TryNextMethod(); fi;
 return PTHPOWERIMAGE_PPI_VEC(fam!.fullSCAlgebra,fam!.zerocoeff,Characteristic(fam),fam!.basisVectors,fam!.pMapping,ExtRepOfObj(x),x);
end);

InstallMethod( PthPowerImage, "for an element of a restricted Lie algebra and an integer",
 [ IsJacobianElement, IsInt ],
 function(x,n)
 local fam;
 fam := FamilyObj(x);
 if not IsBound(fam!.pMapping) then TryNextMethod(); fi;
 while n>0 do
 x := PTHPOWERIMAGE_PPI_VEC(fam!.fullSCAlgebra,fam!.zerocoeff,Characteristic(fam),fam!.basisVectors,fam!.pMapping,ExtRepOfObj(x),x);
 n := n-1;
 od;
 return x;
end);

InstallMethod( PClosureSubalgebra, "for a subalgebra of restricted jacobian elements",
 [ IsLieAlgebra and IsJacobianElementCollection ],
 function(A)
 local i, oldA;

 repeat
 oldA := A;
 for i in Basis(oldA) do
 A := ClosureLeftModule(A,PthPowerImage(i));
 od;
 until A=oldA;
 return A;
end);

#############################################################################
##
#M PthPowerImages( ) . . . . . . . . . . . for a basis of a Lie algebra
##
InstallMethod( PthPowerImages,
 "for a basis of a Lie algebra",
 true,
 [ IsBasis ], 0,
 function( B )

 local L, # the underlying algebra
 p, # the characteristic of `L'
 adL, # a basis of the matrix space spanned by ad L
 basL; # the list of basis vectors `b' of `B' such that
 # `ad b' is a basis vector of `adL'

 L:= UnderlyingLeftModule( B );
 if not IsRestrictedLieAlgebra( L ) then
 Error( "<L> must be a restricted Lie algebra" );
 fi;

 p:= Characteristic( LeftActingDomain( L ) );

 adL:= AdjointBasis( B );

 if Dimension( UnderlyingLeftModule( adL ) ) = 0 then

 # The algebra is abelian.
 return List( BasisVectors( B ), x -> Zero( L ) );

 fi;

 # Now `IndicesOfAdjointBasis( adL )' is a list of indices with `i'-th
 # entry the position of the basis vector of `B'
 # whose adjoint matrix is the `i'-th basis vector of `adL'.
 basL:= BasisVectors( B ){ IndicesOfAdjointBasis( adL ) };

 # We calculate the coefficients of $x_i^{[p]}$ wrt the basis basL.
 return List( BasisVectors( B ),
 x -> Coefficients( adL, AdjointMatrix( B, x ) ^ p ) * basL );
#T And why do you compute the adjoint matrices again?
#T Aren't they stored as basis vectors in adL ?
 end );

############################################################################
##
#M CartanSubalgebra( <L> )
##
## A Cartan subalgebra of the Lie algebra <L> is by definition a nilpotent
## subalgebra equal to its own normalizer in <L>.
##
## By definition, an Engel subalgebra of <L> is the generalized eigenspace
## of a non nilpotent element, corresponding to the eigenvalue 0.
## In a restricted Lie algebra of characteristic p we have that every Cartan
## subalgebra of an Engel subalgebra of <L> is a Cartan subalgebra of <L>.
## Hence in this case we construct a decreasing series of Engel subalgebras.
## When it becomes stable we have found a Cartan subalgebra.
## On the other hand, when <L> is not restricted and is defined over a field
## $F$ of cardinality greater than the dimension of <L> we can proceed as
## follows.
## Let $a$ be a non nilpotent element of <L> and $K$ the corresponding
## Engel subalgebra. Furthermore, let $b$ be a non nilpotent element of $K$.
## Then there is an element $c \in F$ such that $a + c ( b - a )$ has an
## Engel subalgebra strictly contained in $K$
## (see Humphreys, proof of Lemma A, p 79).
##
InstallMethod( CartanSubalgebra,
 "for a Lie algebra",
 true,
 [ IsLieAlgebra ], 0,
 function( L )

 local n, # the dimension of L
 F, # coefficients domain of `L'
 root, # prim. root of `F' if `F' is finite
 K, # a subalgebra of L (on termination a Cartan subalg.)
 a,b, # (non nilpotent) elements of L
 A, # matrix of the equation system (ad a)^n(x)=0
 bas, # basis of the solution space of Ax=0
 sp, # the subspace of L generated by bas
 found,ready, # boolean variables
 c, # an element of `F'
 newelt, # an element of L of the form a+c*(b-a)
 i; # loop variable

 n:= Dimension(L);
 F:= LeftActingDomain( L );

 if IsRestrictedLieAlgebra( L ) then

 K:= L;
 while true do

 a:= NonNilpotentElement( K );

 if a = fail then
 # `K' is a nilpotent Engel subalgebra, hence a Cartan subalgebra.
 return K;
 fi;

 # `a' is a non nilpotent element of `K'.
 # We construct the generalized eigenspace of this element w.r.t.
 # the eigenvalue 0. This is a subalgebra of `K' and of `L'.
 A:= TransposedMat( AdjointMatrix( Basis( K ), a));
 A:= A ^ Dimension( K );
 bas:= NullspaceMat( A );
 bas:= List( bas, x -> LinearCombination( Basis( K ), x ) );
 K:= SubalgebraNC( L, bas, "basis");

 od;

 elif n < Size( F ) then

 # We start with an Engel subalgebra. If it is nilpotent
 # then it is a Cartan subalgebra and we are done.
 # Otherwise we make it smaller.

 a:= NonNilpotentElement( L );

 if a = fail then
 # `L' is nilpotent.
 return L;
 fi;

 # `K' will be the Engel subalgebra corresponding to `a'.

 A:= TransposedMat( AdjointMatrix( Basis( L ), a ) );
 A:= A^n;
 bas:= NullspaceMat( A );
 bas:= List( bas, x -> LinearCombination( Basis( L ), x ) );
 K:= SubalgebraNC( L, bas, "basis");

 # We locate a nonnilpotent element in this Engel subalgebra.

 b:= NonNilpotentElement( K );

 # If `b = fail' then `K' is nilpotent and we are done.
 ready:= ( b = fail );

 while not ready do

 # We locate an element $a + c*(b-a)$ such that the Engel subalgebra
 # belonging to this element is smaller than the Engel subalgebra
 # belonging to `a'.
 # We do this by checking a few values of `c'
 # (At most `n' values of `c' will not yield a smaller subalgebra.)

 sp:= VectorSpace( F, BasisVectors( Basis(K) ), "basis");
 found:= false;
 if Characteristic( F ) = 0 then
 c:= 0;
 else
 root:= PrimitiveRoot( F );
 c:= root;
 fi;
 while not found do

 if Characteristic( F ) = 0 then
 c:= c+1;
 else
 c:= c*root;
 fi;
 newelt:= a+c*(b-a);

 # Calculate the Engel subalgebra belonging to `newelt'.
 A:= TransposedMat( AdjointMatrix( Basis( K ), newelt ) );
 A:= A^Dimension( K );
 bas:= NullspaceMat( A );
 bas:= List( bas, x -> LinearCombination( Basis( K ), x ) );

 # We have found a smaller subalgebra if the dimension is smaller
 # and new basis elements are contained in the old space.

 found:= Length( bas ) < Dimension( K );
 i:= 1;
 while i <= Length( bas ) and found do
 if not bas[i] in sp then
 found:= false;
 fi;
 i:= i+1;
 od;
 od;

 a:= newelt;
 K:= SubalgebraNC( L, bas, "basis");
 b:= NonNilpotentElement( K );

 # If `b = fail' then `K' is nilpotent and we are done.
 ready:= b = fail;

 od;

 return K;

 else

 # the field over which <L> is defined is too small
 TryNextMethod();

 fi;
 end );

##############################################################################
##
#M AdjointAssociativeAlgebra( <L>, <K> )
##
## This function calculates a basis of the associative matrix algebra
## generated by ad_L K, where <K> is a subalgebra of <L>.
## If {x_1,\ldots ,x_n} is a basis of K, then this algebra is spanned
## by all words
## ad x_{i_1}\cdots ad x_{i_t}
## where t>0.
## The degree of such a word is t.
## The algorithm first calculates a maximal linearly independent set
## of words of degree 1, then of degree 2 and so on.
## Since ad x ad y -ady ad x = ad [x,y], we have that we only have
## to consider words where i_1\leq i_2\leq \cdots \leq i_t.
##
InstallMethod( AdjointAssociativeAlgebra,
 "for a Lie algebra and a subalgebra",
 true,
 [ IsAlgebra and IsLieAlgebra, IsAlgebra and IsLieAlgebra ], 0,
 function( L, K )

 local n, # the dimension of L
 F, # the field of L
 asbas, # a list containing the basis elts. of the assoc. alg.
 highdeg, # a list of the elements of the highest degree computed
 # so far
 degree1, # a list of elements of degree 1 (i.e. ad x_i)
 lowinds, # a list of indices such that lowinds[i] is the smallest
 # index in the word highdeg[i]
 hdeg, # the new highdeg constructed each step
 linds, # the new lowinds constructed each step
 i,j,k, # loop variables
 ind, # an index
 m, # a matrix
 posits, # a list of positions in matrices:
 # posits[i] is a list of the form [p,q] such that
 # the matrix asbas[i] has a nonzero entry at position
 # [p][q] and furthermore the matrices asbas[j] with j>i
 # will have a zero at that position (so the basis
 # constructed will be in `upper triangular form')
 l1,l2, # loop variables
 found; # a boolean

 F:= LeftActingDomain( L );

 if Dimension( K ) = 0 then
 return Algebra( F, [ [ [ Zero(F) ] ] ] );
 elif IsLieAbelian( L ) then
 return Algebra( F, [ AdjointMatrix( Basis( L ),
 GeneratorsOfAlgebra( K )[1] ) ] );
 fi;

 n:= Dimension( L );

 # Initialisations that ensure that the first step of the loop will select
 # a maximal linearly independent set of matrices of degree 1.

 degree1:= List( BasisVectors( Basis(K) ),
 x -> AdjointMatrix( Basis(L), x ) );
 posits := [ [ 1, 1 ] ];
 highdeg := [ IdentityMat( n, F ) ];
 asbas := [ Immutable( highdeg[1] ) ];
 lowinds := [ Dimension( K ) ];

 # If after some steps all words of degree t (say) can be reduced modulo
 # lower degree, then all words of degree >t can be reduced to linear
 # combinations of words of lower degree.
 # Hence in that case we are done.

 while not IsEmpty( highdeg ) do

 hdeg:= [];
 linds:= [];

 for i in [1..Length(highdeg)] do

 # Now we multiply all elements `highdeg[i]' with all possible
 # elements of degree 1 (i.e. elements having an index <= the lowest
 # index of the word `highdeg[i]')

 ind:= lowinds[i];
 for j in [1..ind] do

 m:= degree1[j]*highdeg[i];

 # Now we first reduce `m' on the basis computed so far
 # and then add it to the basis.

 for k in [1..Length(posits)] do
 l1:= posits[k][1];
 l2:= posits[k][2];
 m:= m-(m[l1][l2]/asbas[k][l1][l2])*asbas[k];
 od;

 if not IsZero( m ) then

 #'m' is not an element of the span of `asbas'

 Add( hdeg, m );
 Add( linds, j );
 Add( asbas, m);

 # Now we look for a nonzero entry in `m'
 # and add the position of that entry to `posits'.

 found:= false;
 l1:= 1; l2:= 1;
 while not found do
 if not IsZero( m[l1][l2] ) then
 Add( posits, [l1,l2] );
 found:= true;
 else
 if l2 = n then
 l1:= l1+1;
 l2:= 1;
 else
 l2:= l2+1;
 fi;
 fi;
 od;

 fi;

 od;

 od;

 if lowinds = [Dimension(K)] then

 # We are in the first step and hence `degree1' must be made
 # equal to the linearly independent set that we have just calculated.

 degree1:= ShallowCopy( hdeg );
 linds:= [1..Length(degree1)];

 fi;

 highdeg:= ShallowCopy( hdeg );
 lowinds:= ShallowCopy( linds );

 od;

 return Algebra( F, asbas, "basis" );
 end );

##############################################################################
##
#M LieNilRadical( <L> )
##
## Let $p$ be the characteristic of the coefficients field of <L>.
## If $p=0$ the we use the following characterisation of the LieNilRadical:
## Let $S$ be the solvable radical of <L>. And let $H$ be a Cartan subalgebra
## of $S$. Decompose $S$ as $S = H \oplus S_1(H)$, where $S_1(H)$ is the
## Fitting 1-component of the adjoint action of $H$ on $S$. Let $H*$ be the
## associative algebra generated by $ad H$, then $S_1(H)$ is the intersection
## of the spaces $H*^i( S )$ for $i>0$. Let $R$ be the radical of the
## algebra $H*$. Then the LieNilRadical of <L> consists of $S_1(H)$ together
## with all elements $x$ in $H$ such that $ad x\in R$. This last space
## is also characterised as the space of all elements $x$ such that
## $ad x$ lies in the vector space spanned by all nilpotent parts of all
## $ad h$ for $h\in H$.
##
## In the case where $p>0$ we calculate the radical of the associative
## matrix algebra $A$ generated by $ad `L'$.
## The nil radical is then equal to $\{ x\in L \mid ad x \in A \}$.
##
InstallMethod( LieNilRadical,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )

 local F, # the coefficients domain of `L'
 p, # the characteristic of `F'
 bv, # basis vectors of a basis of `L'
 S, # the solvable radical of `L'
 H, # Cartan subalgebra of `S'
 HS, # Fitting 1-component of `S' wrt `H'
 adH, # list of ad x for x in a basis of `H'
 n, # dimension of `L'
 t, # the dimension of an ideal
 eqs, # equation set
 I, # basis vectors of an ideal of `L'
 i,j,k, # loop variables
 sol, # solution set
 adL, # list of matrices ad x where x runs through a basis of
 # `L'
 A, # an associative algebra
 R, # the radical of this algebra
 B; # list of basis vectors of R

 F:= LeftActingDomain( L );
 p:= Characteristic( F );

 if p = 0 then

 # The LieNilRadical of <L> is equal to
 # the LieNilRadical of its solvable radical.

 S:= LieSolvableRadical( L );
 n:= Dimension( S );

 if n in [ 0, 1 ] then return S; fi;

 H:= CartanSubalgebra(S);

 if Dimension(H) = n then return S; fi;

# We calculate the Fitting 1-component $S_1(H)$.

 HS:= ProductSpace( H, S );
 while Dimension( HS ) + Dimension( H ) <> n do
 HS:= ProductSpace( H, HS );
 od;

 if Dimension( H ) = 1 then
 return IdealNC( L, BasisVectors(Basis(HS)), "basis" );
 fi;

# Now we compute the intersection of `R' and `<ad H>'.

 adH:= List( BasisVectors(Basis(H)), x -> AdjointMatrix(Basis(S),x));
 R:= RadicalOfAlgebra( AdjointAssociativeAlgebra( S, H ) );
 B:= BasisVectors( Basis( R ) );

 eqs:= NullMat(Dimension(H)+Dimension(R),n^2,F);
 for i in [1..n] do
 for j in [1..n] do
 for k in [1..Dimension(H)] do
 eqs[k][j+(i-1)*n]:= adH[k][i][j];
 od;
 for k in [1..Dimension(R)] do
 eqs[Dimension(H)+k][j+(i-1)*n]:= B[k][i][j];
 od;
 od;
 od;
 sol:= NullspaceMat( eqs );
 I:= List( sol, x-> LinearCombination( Basis(H), x{[1..Dimension(H)]} ) );

 Append( I, BasisVectors( Basis( HS ) ) );

 return IdealNC( L, I, "basis" );

 else

 n:= Dimension( L );
 bv:= BasisVectors( Basis(L) );
 adL:= List( bv, x -> AdjointMatrix(Basis(L),x) );
 A:= AdjointAssociativeAlgebra( L, L );
 R:= RadicalOfAlgebra( A );

 if Dimension( R ) = 0 then

 # In this case the intersection of `ad L' and `R' is the centre of L.
 return LieCentre( L );

 fi;

 B:= BasisVectors( Basis( R ) );
 t:= Dimension( R );

 # Now we compute the intersection of `R' and `<ad L>'.

 eqs:= NullMat(n+t,n*n,F);
 for i in [1..n] do
 for j in [1..n] do
 for k in [1..n] do
 eqs[k][j+(i-1)*n]:= adL[k][i][j];
 od;
 for k in [1..t] do
 eqs[n+k][j+(i-1)*n]:= -B[k][i][j];
 od;
 od;
 od;
 sol:= NullspaceMat( eqs );
 I:= List( sol, x-> LinearCombination( bv, x{[1..n]} ) );
 return SubalgebraNC( L, I, "basis" );

 fi;

 end );

##############################################################################
##
#M LieSolvableRadical( <L> )
##
## In characteristic zero, the solvable radical of the Lie algebra <L> is
## just the orthogonal complement of $[ <L> <L> ]$ w.r.t. the Killing form.
##
## In characteristic $p > 0$, the following fact is used:
## $R( <L> / NR( <L> ) ) = R( <L> ) / NR( <L> )$ where
## $R( <L> )$ denotes the solvable radical of $L$ and $NR( <L> )$ its
## nil radical).
##
InstallMethod( LieSolvableRadical,
 "for a Lie algebra",
 true,
 [ IsLieAlgebra ], 0,
 function( L )

 local LL, # the derived algebra of L
 n, # the nil radical of L
 B, # a basis of the solvable radical of L
 quo, # the quotient L/n
 r1, # the solvable radical of L/n
 hom; # the canonical map L -> L/n

 if Characteristic( LeftActingDomain( L ) ) = 0 then

 LL:= LieDerivedSubalgebra( L );
 B:= BasisVectors( Basis( KappaPerp( L, LL ) ) );

 else

 n:= LieNilRadical( L );

 if Dimension( n ) = 0 or Dimension( n ) = Dimension( L ) then
 return n;
 fi;

 hom:= NaturalHomomorphismByIdeal( L, n );
 quo:= ImagesSource( hom );
 r1:= LieSolvableRadical( quo );
 B:= BasisVectors( Basis( r1 ) );
 B:= List( B, x -> PreImagesRepresentative( hom, x ) );
 Append( B, BasisVectors( Basis( n ) ) );

 fi;

 SetIsLieSolvable( L, Length( B ) = Dimension( L ) );

 return IdealNC( L, B, "basis");

 end );

##############################################################################
##
#M DirectSumDecomposition( <L> )
##
## This function calculates a list of ideals of `L' such that `L' is equal
## to the direct sum of them.
## The existence of a decomposition of `L' is equivalent to the existence
## of a nontrivial idempotent in the centralizer of `ad L' in the full
## matrix algebra `M_n(F)'. In the general case we try to find such
## idempotents.
## In the case where the Killing form of `L' is nondegenerate we can use
## a more elegant method. In this case the action of the Cartan subalgebra
## will `identify' the direct summands.
##
InstallMethod( DirectSumDecomposition,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )

 local F, # The field of `L'.
 BL, # basis of `L'
 bvl, # basis vectors of `BL'
 n, # The dimension of `L'.
 m, # An integer.
 set, # A list of integers.
 C, # The centre of `L'.
 bvc, # basis vectors of a basis of `C'
 D, # The derived subalgebra of `L'.
 CD, # The intersection of `C' and `D'.
 H, # A Cartan subalgebra of `L'.
 BH, # basis of `H'
 B, # A list of bases of subspaces of `L'.
 cf, # Coefficient list.
 comlist, # List of commutators.
 ideals, # List of ideals.
 bb, # List of basis vectors.
 B1,B2, # Bases of the ideals.
 sp, # A vector space.
 x, # An element of `sp'.
 b, # A list of basis vectors.
 bas,res, # Bases of associative algebras.
 i,j,k,l, # Loop variables.
 centralizer, # The centralizer of `adL' in the matrix algebra.
 Rad, # The radical of `centralizer'.
 M,mat, # Matrices.
 facs, # A list of factors of a polynomial.
 f, # Polynomial.
 contained, # Boolean variable.
 adL, # A basis of the matrix space `ad L'.
 Q, # The factor algebra `centralizer/Rad'
 q, # Number of elements of the field of `L'.
 ei,ni,E, # Elements from `centralizer'
 hom, # A homomorphism.
 id, # A list of idempotents.
 vv; # A list of vectors.

 F:= LeftActingDomain( L );
 n:= Dimension( L );
 if n=0 then
 return [ L ];
 fi;

 if RankMat( KillingMatrix( Basis( L ) ) ) = n then

 # The algorithm works as follows.
 # Let `H' be a Cartan subalgebra of `L'.
 # First we decompose `L' into a direct sum of subspaces `B[i]'
 # such that the minimum polynomial of the adjoint action of an element
 # of `H' restricted to `B[i]' is irreducible.
 # If `L' is a direct sum of ideals, then each of these subspaces
 # will be contained in precisely one ideal.
 # If the field `F' is big enough then we can look for a splitting
 # element in `H'.
 # This is an element `h' such that the minimum polynomial of `ad h'
 # has degree `dim L - dim H + 1'.
 # If the size of the field is bigger than `2*m' then there is a
 # powerful randomised algorithm (Las Vegas type) for finding such an
 # element. We just take a random element from `H' and with probability
 # > 1/2 this will be a splitting element.
 # If the field is small, then we use decomposable elements instead.

 H:= CartanSubalgebra( L );
 BH:= Basis( H );
 BL:= Basis( L );

 m:= (( n - Dimension(H) ) * ( n - Dimension(H) + 2 )) / 8;

 if 2*m < Size(F) and ( not Characteristic( F ) in [2,3] ) then

 set:= [ -m .. m ];

 repeat
 cf:= List([ 1 .. Dimension( H ) ], x -> Random( set ) );
 x:= LinearCombination( BH, cf );
 M:= AdjointMatrix( BL, x );
 f:= CharacteristicPolynomial( F, F, M );
 f:= f/Gcd( f, Derivative( f ) );
 until DegreeOfLaurentPolynomial( f )
 = Dimension( L ) - Dimension( H ) + 1;

 # We decompose the action of the splitting element:

 facs:= Factors( PolynomialRing( F ), f );
 B:= [];
 for i in facs do
 Add( B, List( NullspaceMat( TransposedMat( Value( i, M ) ) ),
 x -> LinearCombination( BL, x ) ) );
 od;

 B:= Filtered( B, x -> not ( x[1] in H ) );

 else

 # Here `L' is a semisimple Lie algebra over a small field or a field
 # of characteristic 2 or 3. This means that
 # the existence of splitting elements is not assured. So we work
 # with decomposable elements rather than with splitting ones.
 # A decomposable element is an element from the associative
 # algebra `T' generated by `ad H' that has a reducible minimum
 # polynomial. Let `V' be a stable subspace (under the action of `H')
 # computed in the process. Then we proceed as follows.
 # We choose a random element from `T' and restrict it to `V'. If this
 # element has an irreducible minimum polynomial of degree equal to
 # the dimension of the associative algebra `T' restricted to `V',
 # then `V' is irreducible. On the other hand,
 # if this polynomial is reducible, then we decompose `V'.

 # `bas' will be a basis of the associative algebra generated by
 # `ad H'. The computation of this basis is facilitated by the fact
 # that we know the dimension of this algebra.

 bas:= List( BH, x -> AdjointMatrix( Basis( L ), x ) );
 sp:= MutableBasis( F, bas );

 k:=1; l:=1;
 while k<=Length(bas) do
 if Length(bas)=Dimension(L)-Dimension(H) then break; fi;
 M:= bas[ k ]*bas[ l ];
 if not IsContainedInSpan( sp, M ) then
 CloseMutableBasis( sp, M );
 Add( bas, M );
 fi;
 if l < Length(bas) then l:=l+1;
 else k:=k+1; l:=1;
 fi;
 od;
 Add( bas, Immutable( IdentityMat( Dimension( L ), F ) ) );

 # Now `B' will be a list of subspaces of `L' stable under `H'.
 # We stop once every element from `B' is irreducible.

 cf:= AsList( F );
 B:= [ ProductSpace( H, L ) ];
 k:= 1;

 while k <= Length( B ) do
 if Dimension( B[k] ) = 1 then
 k:=k+1;
 else
 b:= BasisVectors( Basis( B[k] ) );
 M:= LinearCombination( bas, List( bas, x -> Random( cf ) ) );

 # Now we restrict `M' to the space `B[k]'.

 mat:= [ ];
 for i in [1..Length(b)] do
 x:= LinearCombination( BL, M*Coefficients( BL, b[i] ) );
 Add( mat, Coefficients( Basis( B[k], b ), x ) );
 od;
 M:= TransposedMat( mat );

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

 if Length(facs)=1 then

 # We restrict the basis `bas' to the space `B[k]'. If the length
 # of the result is equal to the degree of `f' then `B[k]' is
 # irreducible.

 sp:= MutableBasis( F,
 [ Immutable( IdentityMat( Dimension(B[k]), F ) ) ] );
 for j in [1..Length(bas)] do
 mat:= [ ];
 for i in [1..Length(b)] do
 x:= LinearCombination( BL, bas[j]*Coefficients( BL, b[i] ) );
 Add( mat, Coefficients( Basis( B[k], b ), x ) );
 od;
 mat:= TransposedMat( mat );

 if not IsContainedInSpan( sp, mat ) then
 CloseMutableBasis( sp, mat );
 fi;

 od;
 res:= BasisVectors( sp );

 if Length( res ) = DegreeOfLaurentPolynomial( f ) then

 # The space is irreducible.

 k:=k+1;

 fi;
 else

 # We decompose.

 for i in facs do
 vv:= List( NullspaceMat( TransposedMat( Value( i, M ) ) ),
 x -> LinearCombination( b, x ) );
 sp:= VectorSpace( F, vv );
 if not sp in B then Add( B, sp ); fi;
 od;

 # We remove the old space from the list;

 B:= Filtered( B, x -> (x <> B[k]) );

 fi;
 fi;

 od;

 B:= List( B, x -> BasisVectors( Basis( x ) ) );
 fi;

 # Now the pieces in `B' are grouped together.

 ideals:=[];

 while B <> [ ] do

 # Check whether `B[1]' is contained in any of
 # the ideals obtained so far.

 contained := false;
 i:=1;
 while not contained and i <= Length(ideals) do
 if B[1][1] in ideals[i] then
 contained:= true;
 fi;
 i:=i+1;
 od;

 if contained then # we do not need B[1] any more

 B:= Filtered( B, x -> x<> B[1] );

 else

 # `B[1]' generates a new ideal.
 # We form this ideal by taking `B[1]' together with
 # all pieces from `B' that do not commute with `B[1]'.
 # At the end of this process, `bb' will be a list of elements
 # commuting with all elements of `B'.
 # From this it follows that `bb' will generate
 # a subalgebra that is a simple ideal. (No remaining piece of `B'
 # can be in this ideal because in that case this piece would
 # generate a smaller ideal inside this one.)

 bb:= ShallowCopy( B[1] );
 B:= Filtered( B, x -> x<> B[1] );
 i:=1;
 while i<= Length( B ) do

 comlist:= [ ];
 for j in [1..Length(bb)] do
 Append( comlist, List( B[i], y -> bb[j]*y ) );
 od;

 if not ForAll( comlist, IsZero ) then
 Append( bb, B[i] );
 B:= Filtered( B, x -> x <> B[i] );
 i:= 1;
 else
 i:=i+1;
 fi;

 od;

 Add( ideals, SubalgebraNC( L, bb ) );

 fi;

 od;

 return List( ideals,
 I -> IdealNC( L, BasisVectors( Basis( I ) ), "basis" ));

 else

 # First we try to find a central component, i.e., a decomposition
 # `L=I_1 \oplus I_2' such that `I_1' is contained in the center of `L'.
 # Such a decomposition exists if and only if the center of `L' is not
 # contained in the derived subalgebra of `L'.

 C:= LieCentre( L );
 bvc:= BasisVectors( Basis( C ) );

 if Dimension( C ) = Dimension( L ) then

 #Now `L' is abelian; hence `L' is the direct sum of `dim L' ideals.

 return List( bvc, v -> IdealNC( L, [ v ], "basis" ) );

 fi;

 BL:= Basis( L );
 bvl:= BasisVectors( BL );

 if 0 < Dimension( C ) then

 D:= LieDerivedSubalgebra( L );
 CD:= Intersection2( C, D );

 if Dimension( CD ) < Dimension( C ) then

 # The central component is the complement of `C \cap D' in `C'.

 B1:=[];
 k:=1;
 sp:= MutableBasis( F,
 BasisVectors( Basis( CD ) ), Zero( CD ) );
 while Length( B1 ) + Dimension( CD ) <> Dimension( C ) do
 x:= bvc[k];
 if not IsContainedInSpan( sp, x ) then
 Add( B1, x );
 CloseMutableBasis( sp, x );
 fi;
 k:=k+1;
 od;

 # The second ideal is a complement of the central component
 # in `L' containing `D'.
#W next statement modified:
 B2:= ShallowCopy( BasisVectors( Basis( D ) ) );
 k:= 1;
 b:= ShallowCopy( B1 );
 Append( b, B2 );
 sp:= MutableBasis( F, b );
 while Length( B2 )+Length( B1 ) <> n do
 x:= bvl[k];
 if not IsContainedInSpan( sp, x ) then
 Add( B2, x );
 CloseMutableBasis( sp, x );
 fi;
 k:= k+1;
 od;

 ideals:= Flat([
 DirectSumDecomposition(IdealNC( L, B1, "basis" )),
 DirectSumDecomposition(IdealNC( L, B2, "basis" ))
 ]);
 return ideals;

 fi;

 fi;

 # Now we assume that `L' does not have a central component
 # and compute the centralizer of `ad L' in `M_n(F)'.

 adL:= List( bvl, x -> AdjointMatrix( BL, x ) );
 centralizer:= FullMatrixAlgebraCentralizer( F, adL );
 Rad:= RadicalOfAlgebra( centralizer );
 if Dimension( centralizer ) - Dimension( Rad ) = 1 then
 return [ L ];
 fi;

 # We calculate a complete set of orthogonal primitive idempotents
 # in the Abelian algebra `centralizer/Rad'.

 hom:= NaturalHomomorphismByIdeal( centralizer, Rad );
 Q:= ImagesSource( hom );
 SetCentre( Q, Q );
 SetRadicalOfAlgebra( Q, Subalgebra( Q, [ Zero( Q ) ] ) );

 id:= List( CentralIdempotentsOfAlgebra( Q ),
 x->PreImagesRepresentative(hom,x));

 # 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;

 # For every idempotent of `centralizer' we calculate
 # a direct summand of `L'.

 ideals:= List( id, e -> List( TransposedMat( e ), v ->
 LinearCombination( BL, v ) ) );
 ideals:= List( ideals, ii -> BasisVectors(
 Basis( VectorSpace( F, ii ) ) ) );

 return List( ideals, ii ->
 IdealNC( L, ii, "basis" ) );

 fi;

 end );

##############################################################################
##
#M IsSimpleAlgebra( <L> ) . . . . . . . . . . . . . . . . for a Lie algebra
##
## A test whether <L> is simple.
## It works only over fields of characteristic 0.
##
InstallMethod( IsSimpleAlgebra,
 "for a Lie algebra in characteristic zero",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )
 if Characteristic( LeftActingDomain( L ) ) <> 0 then
 TryNextMethod();
 elif RankMat( KillingMatrix( Basis( L ) ) ) <> Dimension( L ) then
 return false;
 else
 return Length( DirectSumDecomposition( L ) ) = 1;
 fi;
 end );

##############################################################################
##
#F FindSl2( <L>, <x> )
##
InstallGlobalFunction( FindSl2, function( L, x )

 local n, # the dimension of `L'
 F, # the field of `L'
 B, # basis of `L'
 T, # the table of structure constants of `L'
 xc, # coefficient vector
 eqs, # a system of equations
 i,j,k,l, # loop variables
 cij, # the element `T[i][j]'
 b, # the right hand side of the equation system
 v, # solution of the equations
 z, # element of `L'
 h, # element of `L'
 R, # centralizer of `x' in `L'
 BR, # basis of `R'
 Rvecs, # basis vectors of `R'
 H, # the matrix of `ad H' restricted to `R'
 e0, # coefficient vector
 e1, # coefficient vector
 y; # element of `L'

 if not IsNilpotentElement( L, x ) then
 Error( "<x> must be a nilpotent element of the Lie algebra <L>" );
 fi;

 n:= Dimension( L );
 F:= LeftActingDomain( L );
 B:= Basis( L );
 T:= StructureConstantsTable( B );

 xc:= Coefficients( B, x );
 eqs:= NullMat( 2*n, 2*n, F );

 # First we try to find elements `z' and `h' such that `[x,z]=h'
 # and `[h,x]=2x' (i.e., such that two of the three defining equations
 # of sl_2 are satisfied).
 # This results in a system of `2n' equations for `2n' variables.

 for i in [1..n] do
 for j in [1..n] do
 cij:= T[i][j];
 for k in [1..Length(cij[1])] do
 l:= cij[1][k];
 eqs[i][l] := eqs[i][l] + xc[j]*cij[2][k];
 eqs[n+i][n+l]:= eqs[n+i][n+l] + xc[j]*cij[2][k];
 od;
 od;
 eqs[n+i][i]:= One( F );
 od;

 b:= [];
 for i in [1..n] do
 b[i]:= Zero( F );
 b[n+i]:= 2*One( F )*xc[i];
 od;

 v:= SolutionMat( eqs, b );

 if v = fail then
 # There is no sl_2 containing <x>.
 return fail;
 fi;

 z:= LinearCombination( B, v{ [ 1 .. n ] } );
 h:= LinearCombination( B, v{ [ n+1 .. 2*n ] } );

 R:= LieCentralizer( L, SubalgebraNC( L, [ x ] ) );
 BR:= Basis( R );
 Rvecs:= BasisVectors( BR );

 # `ad h' maps `R' into `R'. `H' will be the matrix of that map.

 H:= List( Rvecs, v -> Coefficients( BR, h * v ) );

 # By the proof of the lemma of Jacobson-Morozov (see Jacobson,
 # Lie Algebras, p. 98) there is an element `e1' in `R' such that
 # `(H+2)e1=e0' where `e0=[h,z]+2z'.
 # If we set `y=z-e1' then `x,h,y' will span a subalgebra of `L'
 # isomorphic to sl_2.

 H:= H+2*IdentityMat( Dimension( R ), F );
#T cheaper!

 e0:= Coefficients( BR, h * z + 2*z );
 e1:= SolutionMat( H, e0 );

 if e1 = fail then
 # There is no sl_2 containing <x>.
 return fail;
 fi;

 y:= z-LinearCombination(Rvecs,e1);

 return SubalgebraNC( L, [x,h,y], "basis" );
end );

#############################################################################
##
#M SemiSimpleType( <L> )
##
## This function works for Lie algebras over a field of characteristic not
## 2 or 3, having a nondegenerate Killing form. Such Lie algebras are
## semisimple. They are characterized as direct sums of simple Lie algebras,
## and these have been classified: a simple Lie algebra is either an element
## of the "great" classes of simple Lie algebas (A_n, B_n, C_n, D_n), or
## an exceptional algebra (E_6, E_7, E_8, F_4, G_2). This function finds
## the type of the semisimple Lie algebra `L'. Since for the calculations
## eigenvalues and eigenvectors of the action of a Cartan subalgebra are
## needed, we reduce the Lie algebra mod p (if it is of characteristic 0).
## The p may not divide the determinant of the matrix of the Killing form,
## nor may it divide the last nonzero coefficient of a minimum polynomial
## of an element of the basis of the Cartan subalgebra.
##
InstallMethod( SemiSimpleType,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )

 local CartanInteger, # Function that computes the Cartan integer.
 bvl, # basis vectors of a basis of `L'
 a, # Element of `L'.
 T,S,S1, # Structure constants tables.
 den, # Denominator of a structure constant.
 denoms, # List of denominators.
 i,j,k, # Loop variables.
 scal, # A scalar.
 K, # A Lie algebra.
 BK, # basis of `K'
 d, # The determinant of the Killing form of `K'.
 p, # A prime.
 H, # Cartan subalgebra.
 s, # An integer.
 mp, # List of minimum polynomials.
 F, # Field.
 bas, # List of basis vectors.
 simples, # List of simple subalgebras.
 types, # List of the types of the elements of simples.
 I, # An element of simples.
 BI, # basis of `I'
 bvi, # basis vectors of `BI'
 HI, # Cartan subalgebra of `I'.
 rk, # The rank of `I'.
 adH, # List of adjoint matrices.
 R, # Root system.
 basR, # Basis of `R'.
 posR, # List of the positive roots.
 fundR, # A fundamental system.
 r,r1,r2,rt, # Roots.
 Rvecs, # List of root vectors.
 basH, # List of basis vectors of a Cartan subalg. of `I'
 sp, # Vector space.
 h, # Element of a Cartan subalgebra of `I'.
 cf, # Coefficient.
 issum, # Boolean.
 CM, # Cartan Matrix.
 endpts; # The endpoints of the Dynkin diagram of `I'.

 if Characteristic( LeftActingDomain( L ) ) in [ 2, 3 ] then
 Info( InfoAlgebra, 1,
 "The field of <L> must not have characteristic 2 or 3." );
 return fail;
 fi;

 # The following function computes the Cartan integer of two roots
 # `r1' and `r2'.
 # If `s' and `t' are the largest integers such that `r1 - s*r2' and
 # `r1 + t*r2' are elements of the root system `R',
 # then the Cartan integer of `r1' and `r2' is `s-t'.

 CartanInteger := function( R, r1, r2 )

 local R1,s,t,rt;

 R1:= ShallowCopy( R );
 Add( R1, R[1]-R[1] );
 s:= 0;
 t:= 0;
 rt:= r1-r2;
 while rt in R1 do
 rt:= rt-r2;
 s:= s+1;
 od;

 rt:= r1+r2;
 while rt in R1 do
 rt:= rt+r2;
 t:= t+1;
 od;
 return s-t;
 end;

 # We test whether the Killing form of `L' is nondegenerate.

 if RankMat( KillingMatrix( Basis( L ) ) ) <> Dimension( L ) then
 Info( InfoAlgebra, 1,
 "The Killing form of <L> is degenerate." );
 return fail;
 fi;

 # First we produce a basis of `L' such that the first basis elements
 # form a basis of a Cartan subalgebra of `L'. Then if `L' is defined
 # over a field of characteristic 0 we do the following. We
 # multiply by an integer in order to ensure that the structure
 # constants are integers.
 # Finally we reduce modulo an appropriate prime `p'.

 H:= CartanSubalgebra( L );
 rk:= Dimension( H );
 bas:= ShallowCopy( BasisVectors( Basis( H ) ) );
 sp:= MutableBasis( LeftActingDomain( L ), bas );
 k:= 1;
 bvl:= BasisVectors( Basis( L ) );
 while Length( bas ) < Dimension( L ) do
 a:= bvl[k];
 if not IsContainedInSpan( sp, a ) then
 Add( bas, a );
 CloseMutableBasis( sp, a );
 fi;
 k:= k+1;
 od;
 T:= StructureConstantsTable( BasisNC( L, bas ) );

 p:= Characteristic( LeftActingDomain( L ) );

 if p = 0 then

 denoms:=[];
 for i in [1..Dimension(L)] do
 for j in [1..Dimension(L)] do
 for k in [1..Length(T[i][j][2])] do
 den:= DenominatorRat( T[i][j][2][k] );
 if (den <> 1) and (not den in denoms) then
 Add( denoms, den );
 fi;
 od;
 od;
 od;

 if denoms <> [] then
 S:= EmptySCTable( Dimension( L ), 0, "antisymmetric" );
 scal:= Lcm( denoms );
 for i in [1..Dimension(L)] do
 for j in [1..Dimension(L)] do
 S[i][j]:= [T[i][j][1],scal*T[i][j][2]];
 od;
 od;
 else
 S:=T;
 fi;

 K:= LieAlgebraByStructureConstants( LeftActingDomain( L ), S );

 BK:= Basis( K );
 d:= DeterminantMat( KillingMatrix( BK ) );
 F:= LeftActingDomain( L );

# `mp' will be a list of minimum polynomials of basis elements of the
# Cartan subalgebra.

 mp:= List( BasisVectors( BK ){[1..rk]},
 x -> CharacteristicPolynomial( F, F, AdjointMatrix( BK, x ) ) );
 mp:= List( mp, x -> x/Gcd( Derivative( x ), x ) );
 d:= d * Product( List( mp, p ->
 CoefficientsOfLaurentPolynomial(p)[1][1] ) );
 p:= 5;
 s:=7;

 # We determine a prime `p>5' not dividing `d' and an integer `s'
 # such that the minimum polynomials of the basis elements
 # of the Cartan subalgebra will split into linear factors
 # over the field of `p^s' elements,
 # and such that `p^s<=2^16'
 # (the maximum size of a finite field in GAP).

 while p^s > 65536 do

 while d mod p = 0 do
 p:= NextPrimeInt( p );
 od;

 F:= GF( p );

 S1:= EmptySCTable( Dimension( K ), Zero( F ), "antisymmetric" );
 for i in [1..Dimension(K)] do
 for j in [1..Dimension(K)] do
 S1[i][j]:= [S[i][j][1], One( F )*List( S[i][j][2], x -> x mod p)];
 od;
 od;

 K:= LieAlgebraByStructureConstants( F, S1 );
 BK:= Basis( K );
 mp:= List( BasisVectors( BK ){[1..rk]},
 x -> CharacteristicPolynomial( F, F, AdjointMatrix( BK, x ) ) );
 s:= Lcm( Flat( List( mp, p -> List( Factors( p ),
 DegreeOfLaurentPolynomial ) )));

 if p=65521 then p:= 1; fi;

 od;

 if p = 1 then
 Info( InfoAlgebra, 1,
 "We cannot find a small modular splitting field for <L>" );

 return fail;
 fi;

 else

 # Here `L' is defined over a field of characteristic p>0. We determine
 # an integer `s' such that the Cartan subalgebra splits over
 # `GF( p^s )'.

 F:= LeftActingDomain( L );
 K:= LieAlgebraByStructureConstants( F, T );
 BK:= Basis( K );
 mp:= List( BasisVectors( BK ){[1..rk]},
 x -> CharacteristicPolynomial( F, F, AdjointMatrix( BK, x ) ) );
 s:= Lcm( Flat( List( mp, p -> List( Factors( p ),
 DegreeOfLaurentPolynomial ) )));
 s:= s*Dimension( LeftActingDomain( L ) );
 if p^s > 2^16 then
 Info( InfoAlgebra, 1,
 "We cannot find a small modular splitting field for <L>" );

 return fail;
 fi;
 S1:= T;
 fi;

 F:= GF( p^s );
 K:= LieAlgebraByStructureConstants( F, S1 );

 # We already know a Cartan subalgebra of `K'.

 BK:= Basis( K );
 H:= SubalgebraNC( K, BasisVectors( BK ){ [ 1 .. rk ] }, "basis" );
 SetCartanSubalgebra( K, H );

 simples:= DirectSumDecomposition( K );

 types:= "";

 # Now for every simple Lie algebra in simples we have to determine
 # its type.
 # For Lie algebras not equal to B_l, C_l or E_6,
 # this is determined by the dimension and the rank.
 # In the other cases we have to examine the root system.

 for I in simples do

 if not IsEmpty( types ) then
 Append( types, " " );
 fi;

 HI:= Intersection2( H, I );
 rk:= Dimension( HI );

 if Dimension( I ) = 133 and rk = 7 then
 Append( types, "E7" );
 elif Dimension( I ) = 248 and rk = 8 then
 Append( types, "E8" );
 elif Dimension( I ) = 52 and rk = 4 then
 Append( types, "F4" );
 elif Dimension( I ) = 14 and rk = 2 then
 Append( types, "G2" );
 else
 if Dimension( I ) = rk^2 + 2*rk then
 Append( types, "A" ); Append( types, String( rk ) );
 elif Dimension( I ) = 2*rk^2-rk then
 Append( types, "D" ); Append( types, String( rk ) );
 elif Dimension( I ) = 10 then
 Append( types, "B2" );
 else

 # We now determine the list of roots and the corresponding
 # root vectors.
 # Since the minimum polynomials of the elements of the
 # Cartan subalgebra split completely,
 # after the call of DirectSumDecomposition,
 # the root vectors are contained in the basis of `I'.

 BI:= Basis( I );
 bvi:= BasisVectors( BI );
 adH:= List( BasisVectors(Basis(HI)), x->AdjointMatrix(BI,x));
#T better!
 R:=[];
 Rvecs:=[];
 for i in [ 1 .. Dimension( I ) ] do
 rt:= List( adH, x -> x[i][i] );
 if not IsZero( rt ) then
 Add( R, rt );
 Add( Rvecs, bvi[i] );
 fi;
 od;

 # A set of roots `basR' is determined such that the set
 # { [x_r,x_{-r}] | r\in basR } is a basis of `HI'.

 basH:= [ ];
 basR:= [ ];
 sp:= MutableBasis( F, [], Zero(I) );
 i:= 1;
 while Length( basH ) < Dimension( HI ) do
 r:= R[i];
 k:= Position( R, -r );
 h:= Rvecs[i] * Rvecs[k];
 if not IsContainedInSpan( sp, h ) then
 Add( basH, h );
 CloseMutableBasis( sp, h );
 Add( basR, r );
 fi;
 i:= i+1;
 od;

 # `posR' will be the set of positive roots.
 # A root `r' is called positive if in the list
 # [ < r, basR[i] >, i=1...Length(basR) ] the first nonzero
 # coefficient is positive
 # (< r_1, r_2 > is the Cartan integer of r_1 and r_2).

 posR:= [ ];
 for r in R do
 if (not r in posR) and (not -r in posR) then
 cf:= 0;
 i:= 0;
 while cf = 0 do
 i:= i+1;
 cf:= CartanInteger( R, r, basR[i] );
 od;
 if 0 < cf then
 Add( posR, r );
 else
 Add( posR, -r );
 fi;
 fi;
 od;

 # A positive root is a fundamental root if it is not
 # the sum of two other positive roots.

 fundR:= [ ];
 for r in posR do
 issum:= false;
 for r1 in posR do
 for r2 in posR do
 if r = r1+r2 then
 issum:= true;
 break;
 fi;
 od;
 if issum then break; fi;
 od;
 if not issum then
 Add( fundR, r );
 fi;
 od;

 # `CM' will be the matrix of Cartan integers
 # of the fundamental roots.

 CM:= List( fundR,
 ri -> List( fundR, rj -> CartanInteger( R, ri, rj ) ) );

 # Finally the properties of the endpoints determine
 # the type of the root system.

 endpts:= [ ];
 for i in [ 1 .. Length(CM) ] do
 if Number( CM[i], x -> x <> 0 ) = 2 then
 Add( endpts, i );
 fi;
 od;

 if Length( endpts ) = 3 then
 Append( types, "E6" );
 elif Sum( CM[ endpts[1] ] ) = 0 or Sum( CM[ endpts[2] ] ) = 0 then
 Append( types, "C" ); Append( types, String( rk ) );
 else
 Append( types, "B" ); Append( types, String( rk ) );
 fi;

 fi;
 fi;
 od;

 return types;
 end );

##############################################################################
##
#M NonNilpotentElement( <L> )
##
InstallMethod( NonNilpotentElement,
 "for a Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )

 local n, # the dimension of `L'
 F, # the field over which `L' is defined
 bvecs, # a list of the basisvectors of `L'
 D, # a list of elements of `L', forming a basis of a nilpotent
 # subspace
 sp, # the space spanned by `D'
 r, # the dimension of `sp'
 found, # a Boolean variable
 i, j, # loop variables
 b, c, # elements of `L'
 elm; #

 # First rule out some trivial cases.
 n:= Dimension( L );
 if n = 1 or n = 0 then
 return fail;
 fi;

 F:= LeftActingDomain( L );
 bvecs:= BasisVectors( Basis( L ) );

 if Characteristic( F ) <> 0 then

 # `D' will be a basis of a nilpotent subalgebra of L.
 if IsNilpotentElement( L, bvecs[1] ) then
 D:= [ bvecs[1] ];
 else
 return bvecs[1];
 fi;

 # `r' will be the dimension of the span of `D'.
 # If `r = n' then `L' is nilpotent and hence does not contain
 # non nilpotent elements.
 r:= 1;

 while r < n do

 sp:= VectorSpace( F, D, "basis" );

 # We first find an element `b' of `L' that does not lie in `sp'.
 found:= false;
 i:= 2;
 while not found do
 b:= bvecs[i];
 if b in sp then
 i:= i+1;
 else
 found:= true;
 fi;
 od;

 # We now replace `b' by `b * D[i]' if
 # `b * D[i]' lies outside `sp' in order to ensure that
 # `[b,sp] \subset sp'.
 # Because `sp' is a nilpotent subalgebra we only need
 # a finite number of replacement steps.

 i:= 1;
 while i <= r do
 c:= b*D[i];
 if c in sp then
 i:= i+1;
 else
 b:= c;
 i:= 1;
 fi;
 od;

 if IsNilpotentElement( L, b ) then
 Add( D, b );
 r:= r+1;
 else
 return b;
 fi;

 od;

 else

 # Now `char F =0'.
 # In this case either `L' is nilpotent or one of the
 # elements $L.1, \ldots , L.n, L.i + L.j; 1 \leq i < j$
 # is non nilpotent.

 for i in [ 1 .. n ] do
 if not IsNilpotentElement( L, bvecs[i] ) then
 return bvecs[i];
 fi;
 od;

 for i in [ 1 .. n ] do
 for j in [ i+1 .. n ] do
 elm:= bvecs[i] + bvecs[j];
 if not IsNilpotentElement( L, elm ) then
 return elm;
 fi;
 od;
 od;

 fi;

 # A non nilpotent element has not been found,
 # hence `L' is nilpotent.
 return fail;

 end );

############################################################################
##
#M PrintObj( <R> ) . . . . . . . . . . . . . . . . . . for a root system
##
InstallMethod( PrintObj,
 "for a root system",
 true, [ IsRootSystem ], 0,
 function( R )

 if HasCartanMatrix( R ) then
 Print("<root system of rank ",Length(SimpleSystem(R)),">");
 else
 Print("<root system>");
 fi;

end );

############################################################################
##
#M \.( <R>, <name> ) . . . . . . . record component access for a root system
##
InstallMethod( \.,
 "for a root system and a record component",
 true, [ IsRootSystem, IsObject ], 0,
 function( R, name )

 name:= NameRNam( name );
 if name = "roots" then
 return Concatenation( PositiveRoots(R), NegativeRoots(R) );
 elif name = "rootvecs" then
 return Concatenation( PositiveRootVectors(R),
 NegativeRootVectors(R) );
 elif name = "fundroots" then
 return SimpleSystem( R );
 elif name = "cartanmat" then
 return CartanMatrix(R);
 else
 TryNextMethod();
 fi;
end );

##############################################################################
##
#M RootSystem( <L> ) . . . . . . . . . . . . . . . . . . . for a Lie algebra
##
InstallMethod( RootSystem,
 "for a (semisimple) Lie algebra",
 true,
 [ IsAlgebra and IsLieAlgebra ], 0,
 function( L )

 local F, # coefficients domain of `L'
 BL, # basis of `L'
 H, # A Cartan subalgebra of `L'
 basH, # A basis of `H'
 sp, # A vector space
 B, # A list of bases of subspaces of `L' whose direct sum
 # is equal to `L'
 newB, # A new version of `B' being constructed
 i,j,l, # Loop variables
 facs, # List of the factors of `p'
 V, # A basis of a subspace of `L'
 M, # A matrix
 cf, # A scalar
 a, # A root vector
 ind, # An index
 basR, # A basis of the root system
 h, # An element of `H'
 posR, # A list of the positive roots
 fundR, # A list of the fundamental roots
 issum, # A boolean
 CartInt, # The function that calculates the Cartan integer of
 # two roots
 C, # The Cartan matrix
 S, # A list of the root vectors
 zero, # zero of `F'
 hts, # A list of the heights of the root vectors
 sorh, # The set `Set( hts )'
 sorR, # The soreted set of roots
 R, # The root system.
 Rvecs, # The root vectors.
 x,y, # Canonical generators.
 noPosR; # Number of positive roots.

 # Let `a' and `b' be two roots of the rootsystem `R'.
 # Let `s' and `t' be the largest integers such that `a-s*b' and `a+t*b'
 # are roots.
 # Then the Cartan integer of `a' and `b' is `s-t'.
 CartInt := function( R, a, b )
 local s,t,rt;
 s:=0; t:=0;
 rt:=a-b;
 while (rt in R) or (rt=0*R[1]) do
 rt:=rt-b;
 s:=s+1;
 od;
 rt:=a+b;
 while (rt in R) or (rt=0*R[1]) do
 rt:=rt+b;
 t:=t+1;
 od;
 return s-t;
 end;

 F:= LeftActingDomain( L );

 if RankMat( KillingMatrix( Basis( L ) ) ) <> Dimension( L ) then
 Info( InfoAlgebra, 1, "the Killing form of <L> is degenerate" );
 return fail;
 fi;

 # First we compute the common eigenvectors of the adjoint action of a
 # Cartan subalgebra `H'. Here `B' will be a list of bases of subspaces
 # of `L' such that `H' maps each element of `B' into itself.
 # Furthermore, `B' has maximal length w.r.t. this property.

 H:= CartanSubalgebra( L );
 BL:= Basis( L );
 B:= [ ShallowCopy( BasisVectors( BL ) ) ];
 basH:= BasisVectors( Basis( H ) );

 for i in basH do

 newB:= [ ];
 for j in B do

 V:= Basis( VectorSpace( F, j, "basis" ), j );
 M:= List( j, x -> Coefficients( V, i*x ) );
 facs:= Factors( PolynomialRing( F ), MinimalPolynomial( F, M ) );

 for l in facs do
 V:= NullspaceMat( Value( l, M ) );
 Add( newB, List( V, x -> LinearCombination( j, x ) ) );
 od;

 od;
 B:= newB;

 od;

 # Now we throw away the subspace `H'.

 B:= Filtered( B, x -> ( not x[1] in H ) );

 # If an element of `B' is not one dimensional then `H' does not split
 # completely, and hence we cannot compute the root system.

 for i in [ 1 .. Length(B) ] do
 if Length( B[i] ) <> 1 then
 Info( InfoAlgebra, 1, "the Cartan subalgebra of <L> in not split" );
 return fail;
 fi;
 od;

 # Now we compute the set of roots `S'.
 # A root is just the list of eigenvalues of the basis elements of `H'
 # on an element of `B'.

 S:= [];
 zero:= Zero( F );
 for i in [ 1 .. Length(B) ] do
 a:= [ ];
 ind:= 0;
 cf:= zero;
 while cf = zero do
 ind:= ind+1;
 cf:= Coefficients( BL, B[i][1] )[ ind ];
 od;
 for j in [1..Length(basH)] do
 Add( a, Coefficients( BL, basH[j]*B[i][1] )[ind] / cf );
 od;
 Add( S, a );
 od;

 Rvecs:= List( B, x -> x[1] );

 # A set of roots `basR' is calculated such that the set
 # { [ x_r, x_{-r} ] | r\in R } is a basis of `H'.

 basH:= [ ];
 basR:= [ ];
 sp:= MutableBasis( F, [], Zero(L) );
 i:=1;
 while Length( basH ) < Dimension( H ) do
 a:= S[i];
 j:= Position( S, -a );
 h:= B[i][1]*B[j][1];
 if not IsContainedInSpan( sp, h ) then
 CloseMutableBasis( sp, h );
 Add( basR, a );
 Add( basH, h );
 fi;
 i:=i+1;
 od;

 # A root `a' is said to be positive if the first nonzero element of
 # `[ CartInt( S, a, basR[j] ) ]' is positive.
 # We calculate the set of positive roots.

 posR:= [ ];
 i:=1;
 while Length( posR ) < Length( S )/2 do
 a:= S[i];
 if (not a in posR) and (not -a in posR) then
 cf:= zero;
 j:= 0;
 while cf = zero do
 j:= j+1;
 cf:= CartInt( S, a, basR[j] );
 od;
 if 0 < cf then
 Add( posR, a );
 else
 Add( posR, -a );
 fi;
 fi;
 i:=i+1;
 od;

 # A positive root is called simple if it is not the sum of two other
 # positive roots.
 # We calculate the set of simple roots `fundR'.

 fundR:= [ ];
 for a in posR do
 issum:= false;
 for i in [1..Length(posR)] do
 for j in [i+1..Length(posR)] do
 if a = posR[i]+posR[j] then
 issum:=true;
 fi;
 od;
 od;
 if not issum then
 Add( fundR, a );
 fi;
 od;

 # Now we calculate the Cartan matrix `C' of the root system.

 C:= List( fundR, i -> List( fundR, j -> CartInt( S, i, j ) ) );

 # Every root can be written as a sum of the simple roots.
 # The height of a root is the sum of the coefficients appearing
 # in that expression.
 # We order the roots according to increasing height.

 V:= BasisNC( VectorSpace( F, fundR ), fundR );
 hts:= List( posR, r -> Sum( Coefficients( V, r ) ) );
 sorh:= Set( hts );

 sorR:= [ ];
 for i in [1..Length(sorh)] do
 Append( sorR, Filtered( posR, r -> hts[Position(posR,r)] = sorh[i] ) );
 od;
 Append( sorR, -1*sorR );
 Rvecs:= List( sorR, r -> Rvecs[ Position(S,r) ] );

 # We calculate a set of canonical generators of `L'. Those are elements
 # x_i, y_i, h_i such that h_i=x_i*y_i, h_i*x_j = c_{ij} x_j,
 # h_i*y_j = -c_{ij} y_j for i \in {1..rank}

 x:= Rvecs{[1..Length(C)]};
 noPosR:= Length( Rvecs )/2;
 y:= Rvecs{[1+noPosR..Length(C)+noPosR]};
 for i in [1..Length(x)] do
 V:= VectorSpace( LeftActingDomain(L), [ x[i] ] );
 B:= Basis( V, [x[i]] );
 y[i]:= y[i]*2/Coefficients( B, (x[i]*y[i])*x[i] )[1];
 od;

 h:= List([1..Length(C)], j -> x[j]*y[j] );

 # Now we construct the root system, and install as many attributes
 # as possible. The roots are represented as lists [ \alpha(h_1),....
 # ,\alpha(h_l)], where the h_i form the `Cartan' part of the canonical
 # generators.

 R:= Objectify( NewType( NewFamily( "RootSystemFam", IsObject ),
 IsAttributeStoringRep and IsRootSystemFromLieAlgebra ),
 rec() );
 SetCanonicalGenerators( R, [ x, y, h ] );
 SetUnderlyingLieAlgebra( R, L );
 SetPositiveRootVectors( R, Rvecs{[1..noPosR]});
 SetNegativeRootVectors( R, Rvecs{[noPosR+1..2*noPosR]} );
 SetCartanMatrix( R, C );

 posR:= [ ];
 for i in [1..noPosR] do
 B:= Basis( VectorSpace( F, [ Rvecs[i] ] ), [ Rvecs[i] ] );
 posR[i]:= List( h, hj -> Coefficients( B, hj*Rvecs[i] )[1] );
 od;

 SetPositiveRoots( R, posR );
 SetNegativeRoots( R, -posR );
 SetSimpleSystem( R, posR{[1..Length(C)]} );

 return R;

 end );

##############################################################################
##
#M CanonicalGenerators( <R> ) . . . . for a root system from a Lie algebra
##
InstallMethod( CanonicalGenerators,
 "for a root system from a (semisimple) Lie algebra",
 true,
 [ IsRootSystemFromLieAlgebra ], 0,
 function( R )

 local L, rank, x, y, i, V, b, c;

 L:= UnderlyingLieAlgebra( R );
 rank:= Length( CartanMatrix( R ) );
 x:= PositiveRootVectors( R ){[1..rank]};
 y:= NegativeRootVectors( R ){[1..rank]};
 for i in [1..Length(x)] do
 V:= VectorSpace( LeftActingDomain(L), [ x[i] ] );
 b:= Basis( V, [x[i]] );
 c:= Coefficients( b, (x[i]*y[i])*x[i] )[1];
 y[i]:= y[i]*2/c;
 od;

 return [ x, y, List([1..rank], j -> x[j]*y[j] ) ];

end );

#############################################################################
##
#M ChevalleyBasis( <L> ) . . . . . . for a semisimple Lie algebra
##
InstallMethod( ChevalleyBasis,
 "for a semisimple Lie algebra with a split Cartan subalgebra",
 true, [ IsLieAlgebra ], 0,
 function( L )

 local R, n, cg, b1p, b1m, b2p, b2m, k, r, i, r1, pos,
 b1, b2, f, cfs, bHa, posRV, negRV, x, y, ha, cf,
 F, T, K, B, BK;

 # We first calculate an automorphism `f' of `L' such that
 # F(L_{\alpha}) = L_{-\alpha}, and f(H)=H, and f acts as multiplication
 # by -1 on H. For this we take the canonical generators of `L',
 # map its `x'-part onto its `y' part (and vice versa), and map
 # the `h'-part on minus itself. The automorphism is determined by this.

 R:= RootSystem( L );
 n:= Length( PositiveRoots( R ) );
 cg:= CanonicalGenerators( R );

 b1p:= ShallowCopy( cg[1] ); b1m:= ShallowCopy( cg[2] );
 b2p:= ShallowCopy( cg[2] ); b2m:= ShallowCopy( cg[1] );

 k:= 1;
 while k <= n do
 r:= PositiveRoots( R )[k];
 for i in [1..Length( CartanMatrix( R ) )] do
 r1:= r + SimpleSystem( R )[i];
 pos:= Position( PositiveRoots( R ), r1 );
 if pos<>fail and not IsBound( b1p[pos] ) then

 b1p[pos]:= cg[1][i]*b1p[k];
 b1m[pos]:= cg[2][i]*b1m[k];
 b2p[pos]:= cg[2][i]*b2p[k];
 b2m[pos]:= cg[1][i]*b2m[k];
 fi;
 od;
 k:= k+1;
 od;

 b1:= b1p; Append( b1, b1m ); Append( b1, cg[3] );
 b2:= b2p; Append( b2, b2m ); Append( b2, -cg[3] );
 f:= LeftModuleHomomorphismByImages( L, L, b1, b2 );

 # Now for every positive root vector `x' we set `y= -Image( f, x )'.
 # We compute a scalar `cf' such that `[x,y]=h', where `h' is the
 # canonical Cartan element corresponding to the root (unquely determined).
 # Then we have to multiply `x' and `y' by Sqrt( 2/cf ), in order to get
 # elements of a Chevalley basis.

 cfs:= [ ];
 bHa:= [ ];
 posRV:= [ ];
 negRV:= [ ];
 for i in [1..n] do
 x:= PositiveRootVectors( R )[i];
 y:= -Image( f, x );
 ha:= x*y;
 cf:= Coefficients( Basis( VectorSpace( LeftActingDomain(L),
 [x] ), [x] ), ha*x )[1];
 if i <= Length( CartanMatrix( R ) ) then Add( bHa, (2/cf)*ha ); fi;
 Add( cfs, Sqrt( 2/cf ) );
 posRV[i]:= x; negRV[i]:= y;
 od;

 # In general the `cfs' will lie in a field extension of the ground field.
 # We construct the Lie algebra over that field with the same structure
 # constants as `L'. Then we map the Chevalley basis elements into
 # this new Lie algebra. Then we take the structure constants table of
 # this new Lie algebra with respect to the Chevalley basis, and
 # form a new Lie algebra over the same field as `L' with this table.

 F:= DefaultField( cfs );
 T:= StructureConstantsTable( Basis( L ) );
 K:= LieAlgebraByStructureConstants( F, T );
 BK:= CanonicalBasis( K );
 B:= [ ];
 for i in [1..n] do
 B[i]:= LinearCombination( BK, cfs[i]*Coefficients( Basis(L),
 posRV[i] ) );
 B[n+i]:= LinearCombination( BK, cfs[i]*Coefficients( Basis(L),
 negRV[i] ) );
 od;
 for i in [1..Length(bHa)] do
 B[2*n+i]:= LinearCombination( BK, Coefficients( Basis(L),
 bHa[i] ) );
 od;

 T:= StructureConstantsTable( Basis( K, B ) );
 K:= LieAlgebraByStructureConstants( LeftActingDomain(L), T );
 B:= BasisVectors( CanonicalBasis( K ) );

 # Now the basis elements of `K' form a Chevalley basis. Furthermore,
 # `K' is isomorphic to `L'. We construct the isomorphism, and map
 # the basis elements of `K' into `L', thus getting a Chevalley basis
 # in `L'.

 b1p:= B{[1..Length(CartanMatrix(R))]};
 b1m:= B{[n+1..n+Length(CartanMatrix(R))]};
 b2p:= ShallowCopy( cg[1] ); b2m:= ShallowCopy( cg[2] );

 k:= 1;
 while k <= n do
 r:= PositiveRoots( R )[k];
 for i in [1..Length( CartanMatrix( R ) )] do
 r1:= r + SimpleSystem( R )[i];
 pos:= Position( PositiveRoots( R ), r1 );
 if pos<>fail and not IsBound( b1p[pos] ) then

 b1p[pos]:= b1p[i]*b1p[k];
 b1m[pos]:= b1m[i]*b1m[k];
 b2p[pos]:= b2p[i]*b2p[k];
 b2m[pos]:= b2m[i]*b2m[k];
 fi;
 od;
 k:= k+1;
 od;

 b1:= b1p; Append( b1, b1m ); Append( b1, B{[2*n+1..Length(B)]} );
 b2:= b2p; Append( b2, b2m ); Append( b2, cg[3] );
 f:= LeftModuleHomomorphismByImages( K, L, b1, b2 );

 return [ List( B{[1..n]}, x -> Image( f, x ) ),
 List( B{[n+1..2*n]}, y -> Image( f, y ) ),
 cg[3] ];

 end);

#############################################################################
##
#F DescriptionOfNormalizedUEAElement( <T>, <listofpairs> )
##
InstallGlobalFunction( DescriptionOfNormalizedUEAElement,
 function( T, listofpairs )

 local normalized, # ordered list of normalized coeff./monom. pairs
 indices, # list that stores at position $i$ up to what
 # position the $i$-th monomial is known to be
 # normalized
 s, i, j, k, l, # loop variables
 2i, # `2*i'
 scalar, # coefficient of the monomial under work
 mon, # monomial under work
 len, # length of the monomial under work
 head, # initial part of the monomial under work
 middle, # middle part of the monomial under work
 tail, # trailing part of the monomial under work
 index, # new value of `indices[i]'
 Tcoeffs, # one entry in `T'
 lennorm, # length of `normalized' at the moment
 zero; # zero coefficient

 normalized := [];

 while not IsEmpty( listofpairs ) do

 listofpairs:= Compacted( listofpairs );

 # `indices' is a list of positive integers $[ j_1, j_2, \ldots, j_m ]$
 # s.t. the initial part $x_{i_1}^{e_1} \cdots x_{i_{j_k}}^{e_{j_k}}$
 # of the $k$-th monomial is known to be normalized,
 # i.e., $i_1 < i_2 < \cdots < i_{j_k}$.
 # (So $j_k = 1$ for all $k$ will always be correct.)
 indices:= ListWithIdenticalEntries( Length( listofpairs )/2, 1 );

 # Loop over the monomials that shall be normalized.
 for i in [ 1, 2 .. Length( indices ) ] do

 # If the `i'-th monomial is already normalized,
 # put it into `normalized'.
 # Otherwise swap the first non-ordered generators.
 2i:= 2*i;
 scalar:= listofpairs[ 2i ];
 mon:= listofpairs[ 2i-1 ];
 len:= Length( mon );
 j:= 2 * indices[i] - 1;
 while j < len - 2 do

 if mon[j] < mon[ j+2 ] then

 # `mon' is better normalized than `indices' tells.
 j:= j+2;
 indices[i]:= indices[i] + 1;

 elif mon[j] = mon[ j+2 ] then

 # absorption
 mon[ j+1 ]:= mon[ j+1 ] + mon[ j+3 ];
 for k in [ j+2 .. len-2 ] do
 mon[k]:= mon[ k+2 ];
 od;
 Unbind( mon[ len ] );
 Unbind( mon[ len-1 ] );
 len:= len - 2;

 else

 # We must swap two generators.
 # First construct head and tail of the arising monomials.
 head:= mon{ [ 1 .. j-1 ] };

 middle:= [ mon[ j+2 ], mon[j+3], mon[j], mon[j+1] ];

 tail:= mon{ [ j+4 .. len ] };

 # Adjust `indices[i]'.
 index:= indices[i] - 1;
 if index = 0 then
 index:= 1;
 fi;

 indices[i]:= index;

 # Replace the monomial by the swapped one.
 listofpairs[ 2i-1 ]:= Concatenation( head, middle, tail );

 # Add the coeffs/monomials that are given by the commutator.
 # The part between `head' and `tail' of these listofpairs is
 # $a_{ji}=\sum_{k=1}^d c_{ijk} x_d$.
 # Here we use the following formula (which is easily proved
 # by induction):
 #
 # x_j^m x_i^n = x_i^n x_j^m + \sum_{l=0}^{m-1} \sum_{k=0}^{n-1}
 # x_j^l x_i^k a_{ji} x_i^{n-1-k} x_j^{m-1-l}
 #
 #
 # where x_jx_i = x_ix_j + a_{ji}
 #
 Tcoeffs:= T[ mon[j] ][ mon[ j+2 ] ];
 for s in [ 1 .. Length( Tcoeffs[1] ) ] do
 for l in [ 0 .. mon[j+1] - 1 ] do
 for k in [ 0 .. mon[j+3] - 1 ] do

 middle:= [ ];

 if l > 0 then
 middle:= [ mon[j], l ];
 fi;
 if k > 0 then
 Append( middle, [ mon[j+2], k ] );
 fi;
 Append( middle, [ Tcoeffs[1][s], 1 ] );
 if mon[j+3]-1-k > 0 then
 Append( middle, [ mon[j+2], mon[j+3]-1-k ] );
 fi;
 if mon[j+1]-1-l > 0 then
 Append( middle, [ mon[j], mon[j+1]-1-l ] );
 fi;

 Append( listofpairs,
 [ Concatenation( head, middle, tail ),
 scalar * Tcoeffs[2][s] ] );
 Add( indices, index );
 od;
 od;
 od;

 break;

 fi;

 od;

 # If the monomial is normalized then move it to `normalized'.
 if len - 2 <= j then

 # Find the correct position in `normalized',
 # and insert the monomial.
 lennorm:= Length( normalized );
 k:= 2;
 while k <= lennorm do
 if listofpairs[ 2i-1 ] < normalized[ k-1 ] then
 for l in [ lennorm, lennorm-1 .. k-1 ] do
 normalized[l+2]:= normalized[l];
 od;
 normalized[ k-1 ]:= listofpairs[ 2i-1 ];
 normalized[ k ]:= scalar;
 break;
 elif listofpairs[ 2i-1 ] = normalized[ k-1 ] then
 normalized[k]:= normalized[k] + scalar;
 break;
 fi;
 k:= k+2;
 od;
 if lennorm < k then
 normalized[ lennorm+1 ]:= listofpairs[ 2i-1 ];
 normalized[ lennorm+2 ]:= scalar;
 fi;

 # Remove the monomial from `listofpairs'.
 Unbind( listofpairs[ 2i-1 ] );
 Unbind( listofpairs[ 2i ] );

 fi;

 od;

 od;

 # Remove monomials with multiplicity zero;
 if not IsEmpty( normalized ) then
 zero:= Zero( normalized[2] );
 for i in [ 2, 4 .. Length( normalized ) ] do
 if normalized[i] = zero then
 Unbind( normalized[ i-1 ] );
 Unbind( normalized[ i ] );
 fi;
 od;
 normalized:= Compacted( normalized );
 fi;

 # Return the normal form.
 return normalized;
end );

#############################################################################
##
#M UniversalEnvelopingAlgebra( <L> ) . . . . . . . . . . . for a Lie algebra
##
InstallOtherMethod( UniversalEnvelopingAlgebra,
 "for a finite dimensional Lie algebra and a basis of it",
 true,
 [ IsLieAlgebra, IsBasis ], 0,
 function( L, B )

 local F, # free associative algebra
 U, # universal enveloping algebra, result
 gen, # loop over algebra generators of `U'
 Fam, # elements family of `U'
 T, # s.c. table of a basis of `L'
 FamMon, # family of monomials
 FamFree; # elements family of `F'

 # Check the argument.
 if not IsFiniteDimensional( L ) then
 Error( "<L> must be finite dimensional" );
 fi;

 # Construct the universal enveloping algebra.
 F:= FreeAssociativeAlgebraWithOne( LeftActingDomain( L ),
 Dimension( L ), "x" );
 U:= FactorFreeAlgebraByRelators( F, [ Zero( F ) ] );
#T do not cheat here!

 # Enter knowledge about `U'.
 SetDimension( U, infinity );
 for gen in GeneratorsOfLeftOperatorRingWithOne( U ) do
 SetIsNormalForm( gen, true );
 od;
 SetIsNormalForm( Zero( U ), true );

 # Enter data to handle elements.
 Fam:= ElementsFamily( FamilyObj( U ) );
 Fam!.normalizedType:= NewType( Fam,
 IsElementOfFpAlgebra
 and IsPackedElementDefaultRep
 and IsNormalForm );

 T:= StructureConstantsTable( B );
 FamMon:= ElementsFamily( FamilyObj( UnderlyingMagma( F ) ) );
 FamFree:= ElementsFamily( FamilyObj( F ) );

 SetNiceNormalFormByExtRepFunction( Fam,
 function( Fam, extrep )
 local zero, i;
 zero:= extrep[1];
 extrep:= DescriptionOfNormalizedUEAElement( T, extrep[2] );
 for i in [ 1, 3 .. Length( extrep ) - 1 ] do
 extrep[i]:= ObjByExtRep( FamMon, extrep[i] );
 od;
 return Objectify( Fam!.normalizedType,
 [ Objectify( FamFree!.defaultType, [ zero, extrep ] ) ] );
 end );

 SetOne( U, ElementOfFpAlgebra( Fam, One( F ) ) );

 # Enter `L'; it is used to set up the embedding (as a vector space).
 Fam!.liealgebra:= L;
#T is not allowed ...

 # Return the universal enveloping algebra.
 return U;
 end );

#T missing: embedding of the Lie algebra (as vector space)
#T missing: relators (only compute them if they are explicitly wanted)
#T (attribute `Relators'?)

InstallMethod( UniversalEnvelopingAlgebra,
 "for a finite dimensional Lie algebra",
 true,
 [ IsLieAlgebra ], 0,
 function( L )

 return UniversalEnvelopingAlgebra( L, Basis(L) );
end );

#############################################################################
##
#F IsSpaceOfUEAElements( <V> )
##
## If <V> is a space of elements of a universal enveloping algebra,
## then the `NiceFreeLeftModuleInfo' value of <V> is a record with the
## following components.
## \beginitems
## `family' &
## the elements family of <V>,
##
## `monomials' &
## a list of monomials occurring in the generators of <V>,
##
##
## `zerocoeff' &
## the zero coefficient of elements in <V>,
##
## `zerovector' &
## the zero row vector in the nice free left module,
##
## `characteristic' &
## the characteristic of the ground field.
## \enditems
## The `NiceVector' value of $v \in <V>$ is defined as the row vector of
## coefficients of $v$ w.r.t. the list `monomials'.
##
##
DeclareHandlingByNiceBasis( "IsSpaceOfUEAElements",
 "for free left modules of elements of a universal enveloping algebra" );

#############################################################################
##
#M NiceFreeLeftModuleInfo( <V> )
#M NiceVector( <V>, <v> )
#M UglyVector( <V>, <r> )
##
InstallHandlingByNiceBasis( "IsSpaceOfUEAElements", rec(
 detect := function( F, gens, V, zero )
 return IsElementOfFpAlgebraCollection( V );
 end,

 NiceFreeLeftModuleInfo := function( V )
 local gens,
 monomials,
 gen,
 list,
 zero,
 info;

 gens:= GeneratorsOfLeftModule( V );

 monomials:= [];

 for gen in gens do
 list:= ExtRepOfObj( gen )[2];
 UniteSet( monomials, list{ [ 1, 3 .. Length( list ) - 1 ] } );
 od;

 zero:= Zero( LeftActingDomain( V ) );
 info:= rec( monomials := monomials,
 zerocoeff := zero,
 characteristic:= Characteristic( LeftActingDomain( V ) ),
 family := ElementsFamily( FamilyObj( V ) ) );

 # For the zero row vector, catch the case of empty `monomials' list.
 if IsEmpty( monomials ) then
 info.zerovector := MakeImmutable([ zero ]);
 else
 info.zerovector := MakeImmutable(ListWithIdenticalEntries( Length( monomials ),
 zero ) );
 fi;

 return info;
 end,

 NiceVector := function( V, v )
 local info, c, monomials, i, pos;
 info:= NiceFreeLeftModuleInfo( V );
 c:= ShallowCopy( info.zerovector );
 v:= ExtRepOfObj( v )[2];
 monomials:= info.monomials;
 for i in [ 2, 4 .. Length( v ) ] do
 pos:= Position( monomials, v[ i-1 ] );
 if pos = fail then
 return fail;
 fi;
 c[ pos ]:= v[i];
 od;
 return c;
 end,

 UglyVector := function( V, r )
 local info, list, i;
 info:= NiceFreeLeftModuleInfo( V );
 if Length( r ) <> Length( info.zerovector ) then
 return fail;
 elif IsEmpty( info.monomials ) then
 if IsZero( r ) then
 return Zero( V );
 else
 return fail;
 fi;
 fi;
 list:= [];
 for i in [ 1 .. Length( r ) ] do
 if r[i] <> info.zerocoeff then
 Add( list, info.monomials[i] );
 Add( list, r[i] );
 fi;
 od;
 return ObjByExtRep( info.family, [ info.characteristic, list ] );
 end ) );

#############################################################################
##
#F FreeLieAlgebra( <R>, <rank> )
#F FreeLieAlgebra( <R>, <rank>, <name> )
#F FreeLieAlgebra( <R>, <name1>, <name2>, ... )
##
InstallGlobalFunction( FreeLieAlgebra, function( arg )

 local R, # coefficients ring
 names, # names of the algebra generators
 M, # free magma
 F, # family of magma ring elements
 one, # identity of `R'
 zero, # zero of `R'
 L; # free Lie algebra, result

 # 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: FreeLieAlgebra( <R>, <rank> )\n",
 "or FreeLieAlgebra( <R>, <name1>, ... )" );
 fi;

 # Construct the algebra as magma algebra modulo relations
 # over a free magma.
 M:= FreeMagma( names );

 # Construct the family of elements of our ring.
 F:= NewFamily( "FreeLieAlgebraObjFamily",
 IsElementOfMagmaRingModuloRelations,
 IsJacobianElement and IsZeroSquaredElement );
 SetFilterObj( F, IsFamilyElementOfFreeLieAlgebra );

 one:= One( R );
 zero:= Zero( R );

 F!.defaultType := NewType( F, IsMagmaRingObjDefaultRep );
 F!.familyRing := FamilyObj( R );
 F!.familyMagma := FamilyObj( M );
 F!.zeroRing := zero;
#T no !!
 F!.names := names;

 # Set the characteristic.
 if HasCharacteristic( R ) or HasCharacteristic( FamilyObj( R ) ) then
 SetCharacteristic( F, Characteristic( R ) );
 fi;

 # Make the magma ring object.
 L:= Objectify( NewType( CollectionsFamily( F ),
 IsMagmaRingModuloRelations
 and IsAttributeStoringRep ),
 rec() );

 # Set the necessary attributes.
 SetLeftActingDomain( L, R );
 SetUnderlyingMagma( L, M );

 # Deduce useful information.
 SetIsFiniteDimensional( L, false );
 if HasIsWholeFamily( R ) and HasIsWholeFamily( M ) then
 SetIsWholeFamily( L, IsWholeFamily( R ) and IsWholeFamily( M ) );
 fi;

 # Construct the generators.
 SetGeneratorsOfLeftOperatorRing( L,
 List( GeneratorsOfMagma( M ),
 x -> ElementOfMagmaRing( F, zero, [ one ], [ x ] ) ) );

 # Install grading
 SetGrading( L, rec( min_degree := 1,
 max_degree := infinity,
 source := Integers,
 hom_components := function(degree)
 local B, d, i, x, y, z;
 B := GeneratorsOfMagma(M);
 B := [List([1..Length(B)],i->[[i],fail,B[i]])];
 for d in [2..degree] do
 Add(B,[]);
 for i in [1..d-1] do
 for x in B[i] do for y in B[d-i] do
 z := Concatenation(x[1],y[1]);
 if z<y[1] and x[1]<y[1] and (x[2]=fail or x[2]>=y[1]) then
 Add(B[d],[z,y[1],x[3]*y[3]]);

 fi;
 od; od;
 od;
 od;
 if degree<1 then B := []; else B := B[degree]; fi;
 return FreeLeftModule( R, List( B,
 p->ElementOfMagmaRing( F, zero, [ one ], [ p[3] ] )), Zero(L));
 end) );
 # Return the ring.
 return L;
end );

#############################################################################
##
#M NormalizedElementOfMagmaRingModuloRelations( <Fam>, <descr> )
##
## <descr> is a list of the form `[ <z>, <list> ]', <z> being the zero
## coefficient of the ring, and <list> being the list of monomials and
## their coefficients.
## This function returns the element described in <descr> expanded on
## the Lyndon basis of the free Lie algebra. In order to do this we do not
## need to know this basis; we only need a test whether something is a
## Lyndon element (this is done in the function `IsLyndonT').
## For the algorithm we refer to C. Reutenauer, Free Lie Algebras, Clarendon
## Press, Oxford, 1993.
##
DeclareGlobalName( "IsLyndonT" );

BindGlobal( "IsLyndonT",
function( t )

 # This function tests whether the bracketed expression `t' is
 # a Lyndon tree.

 local w,w1,w2,b,y;

 if not IsList( t ) then return true; fi;

 w:= Flat( t );
 if IsList( t[1] ) then
 w1:= Flat( t[1] );
 b:= false;
 else
 w1:= [ t[1] ];
 b:=true;
 fi;
 if IsList( t[2] ) then
 w2:= Flat( t[2] );
 else
 w2:= [ t[2] ];
 fi;

 if w<w2 and w1<w2 then
 if not b then
 y:= Flat( [ t[1][2] ] );
 if y < w2 then return false; fi;
 fi;
 else
 return false;
 fi;

 return IsLyndonT(t[1]) and IsLyndonT(t[2]);
end);

InstallMethod( NormalizedElementOfMagmaRingModuloRelations,
 "for family of free Lie algebra elements, and list",
 true,
 [ IsFamilyElementOfFreeLieAlgebra, IsList ], 0,
 function( Fam, descr )

 local todo, #The list of elements that are to be expanded
 k,i, #Loop variables
 z,s,u,v,x,y,w, #Bracketed expressions (or `trees')
 cf, #Coefficient
 found, #Boolean
 ll, #List
 zero, #The zero element of the field
 tlist, #List of elements of the free Lie algebra
 Dcopy; #Two functions

 Dcopy:=function( l )

 if not IsList(l) then return ShallowCopy( l ); fi;
 return List( l, Dcopy );
 end;

 zero:= descr[1];
 todo:= [ ];
 i:= 1;

# Every element in `todo' has the following format: [ bool, cf, br ],
# where bool is a boolean; it is true if the br is a Lyndon tree.
# cf is the coefficient (number) and br is a bracketed expression.
# The reason for `tagging' everything is that almost anywhere in the
# list cancellations may occur. This tagging provides an efficient way
# of remembering which trees were dealt with before.

 while i+1 <= Length(descr[2]) do
 Add( todo, [ false, descr[2][i+1], ExtRepOfObj( descr[2][i] ) ] );
 i:= i+2;
 od;

 k:= 1;
 while k<=Length(todo) do

 s:= todo[k][3];
 cf:= todo[k][2];
 if not IsList( s ) then
 # `s' is a Lyndon tree
 todo[k][1]:= true;
 k:= k+1;
 elif cf = zero or s[1]=s[2] then
 # `s' is zero
 ll:=Filtered([1..Length(todo)], x -> x<> k);
 todo:= todo{ll};
 elif todo[k][1] then
 # we already dealt with `s'
 k:=k+1;
 elif IsLyndonT( s ) then
 # we do not need to expand `s'
 todo[k][1]:=true;
 k:=k+1;
 else
 #we expand `s'
 found:= false; u:=s;
 z:= Dcopy( s );
 v:= z;

 # we look for a subtree `u' such that is not a Lyndon tree
 # such that its left and right subtrees are Lyndon trees.

 while not found do
 if IsLyndonT(u[1]) then
 if IsLyndonT(u[2]) then
 found:=true;
 else
 u:= u[2]; v:=v[2];
 fi;
 else
 u:= u[1]; v:=v[1];
 fi;
 od;

 if u[1]=u[2] then
 # the whole expression `s' reduces to zero.
 ll:= Filtered([1..Length(todo)], x->x<>k);
 todo:= todo{ll};
 else
 if Flat([u[1]]) > Flat([u[2]]) then
 # interchange u[1] and u[2]; this introduces a -.
 w:=u[1];
 u[1]:=u[2];
 u[2]:=w;
 i:= 1;
 found:= false;
 while i<= Length( todo ) and not found do
 if todo[i][3] = s and k<>i then
 todo[i][2]:= todo[i][2]-cf;
 if todo[i][2] = zero then
 ll:=Filtered([1..Length(todo)],x->(x<>k and x<>i ));
 else
 ll:=Filtered([1..Length(todo)],x->x<>k);
 fi;
 todo:= todo{ll};
 found:= true;
 fi;
 i:=i+1;
 od;
 if not found then todo[k][2]:=-todo[k][2]; fi;
 else
 #use the Jacobi identity.
 x:=u[1][1];
 y:=u[1][2];
 w:= u[2];
 u[1]:=[x,w];
 u[2]:=y;

 i:= 1;
 found:= false;
 while i<= Length( todo ) and not found do
 if todo[i][3] = s and k<>i then
 todo[i][2]:= todo[i][2]+cf;
 if todo[i][2] = zero then
 ll:=Filtered([1..Length(todo)],x->(x<>k and x<>i ));
 else
 ll:=Filtered([1..Length(todo)],x->x<>k);
 fi;
 todo:= todo{ll};
 found:= true;
 fi;
 i:=i+1;
 od;

 x:=v[1][1];
 y:=v[1][2];
 w:=v[2];
 v[1]:=x;
 v[2]:=[y,w];

 i:= 1;
 found:= false;
 while i<= Length( todo ) and not found do
 if todo[i][3] = z then
 todo[i][2]:= todo[i][2]+cf;
 found:= true;
 fi;
 i:= i+1;
 od;
 if not found then
 Add( todo, [false,cf,z] );
 fi;
 fi;
 k:=1;
 fi;
 fi;
 od;

# wrap the list `todo' into an element of the free Lie algebra.

 todojava.lang.StringIndexOutOfBoundsException: Range [25, 24) out of bounds for length 35
 Sort java.lang.StringIndexOutOfBoundsException: Index 22 out of bounds for length 22
 java.lang.StringIndexOutOfBoundsException: Range [15, 14) out of bounds for length 20
 MakeConfluent
 Append( ##
 od;

([ ;

 end


##############################################################################
##
#M ImageElm( <h>, <x> )
#M ImagesRepresentative( <h>, <x> )
##
## A special method for calculating the (unique) image of an element <x>
## under an FptoSCAMorphism <h>. The fact that <h> knows the images of the
## generators together with the fact that <h> is an algebra morphism is used
## (rather than the linearity of <h>).
##
(java.lang.StringIndexOutOfBoundsException: Range [37, 36) out of bounds for length 55
,

 EvalProduct(rws)

if prodthen
 ims[prod];
 else
 return EvalProduct( prod[1], ims )*EvalProduct( prod[2], ims );
 fi;
 end;

 e:=MappingGeneratorsImages(h);
 imgs:= e[2];
 e:= ExtRepOfObj(x)[2];
 im:= 0*imgs[1];
 k:= 1;
 while k <= Length(e) do
 im:= im + e[k+1]*EvalProduct( e[k], imgs );
 k:= k+2;
 od;
 return im;
end );

InstallMethod( ImageElm,
 "for Fp to SCA mapping, and element",
 FamSourceEqFamElm,
 [ IsFptoSCAMorphism, IsElementOfFpAlgebra ], 0,
 FptoSCAMorphismImageElm );

InstallMethod( ImagesRepresentative,
 "for Fp to SCA mapping, and element",
 FamSourceEqFamElm,
 [ IsFptoSCAMorphism, IsElementOfFpAlgebra ], 0,
 FptoSCAMorphismImageElm );

###########################################################################
##
#M PreImagesRepresentative( f, x )
##
InstallMethod( PreImagesRepresentative,
 "for Fp to SCA mapping, and element",
 FamRangeEqFamElm,
 [ IsFptoSCAMorphism, IsSCAlgebraObj ], 0,

 function( f, x )

 local dim, e, gens, imgs, b1, b2, levs,
 brackets, sp, deg, newlev, newbracks, d, br1, br2,
 i, j, a, b, c, z, imz, cf;

 if not IsBound( f!.bases ) then

 # We find bases of the source and the range that are mapped to
 # each other.

 dim:= Dimension( Range(f) );
 e:=MappingGeneratorsImages(f);
 gens:= e[1];
 imgs:= e[2];
 b1:= ShallowCopy( gens );
 b2:= ShallowCopy( imgs );
 levs:= [ gens ];
 brackets:= [ [1..Length(gens)] ];
 sp:= MutableBasis( LeftActingDomain(Range(f)), b2 );
 deg:= 1;
 while Length( b1 ) <> dim do
 deg:= deg+1;
 newlev:= [ ];
 newbracks:= [ ];
 # get all Lyndon elements of degree deg:
 for d in [1..Length(brackets)] do
 if Length( b1 ) = dim then break; fi;
 br1:= brackets[d];
 br2:= brackets[deg-d];
 for i in [1..Length(br1)] do
 if Length( b1 ) = dim then break; fi;
 for j in [1..Length(br2)] do
 if Length( b1 ) = dim then break; fi;
 a:= br1[i]; b:= br2[j];
 if IsLyndonT( [a,b] ) then
 c:= [a,b];
 z:= levs[d][i]*levs[deg-d][j];
 elif IsLyndonT( [b,a] ) then
 c:= [b,a];
 z:= levs[deg-d][j]*levs[d][i];
 else
 c:= [ ];
 fi;

 if c <> [] then

 imz:= Image( f, z );
 if not IsContainedInSpan( sp, imz ) then
 CloseMutableBasis( sp, imz );
 Add( b1, z );
 Add( newlev, z );
 Add( newbracks, c );
 Add( b2, imz );
 fi;
 fi;

 od;
 od;
 od;
 Add( levs, newlev );
 Add( brackets, newbracks );
 od;

 f!.bases:= [ b1, Basis( Range(f), b2 ) ];
 fi;

 cf:= Coefficients( f!.bases[2], x );
 return cf*f!.bases[1];

end);

#############################################################################
##
#M Dimension( <FpL> )
##
## A method for the dimension of a finitely-presented Lie algebra.
##
InstallMethod( Dimension,
 "for a f.p. Lie algebra",
 true,
 [ IsLieAlgebra and IsSubalgebraFpAlgebra], 0,
 function( L )
 local h;
 h:= NiceAlgebraMonomorphism( L );
 if h <> fail then
 return Dimension( Range( h ) );
 else
 TryNextMethod();
 fi;
end);

##############################################################################
##
#M IsFiniteDimensional( <FpL> )
##
## For finitely-presented Lie algebras.
##
InstallMethod( IsFiniteDimensional,
 "for a f.p. Lie algebra",
 true,
 [ IsLieAlgebra and IsSubalgebraFpAlgebra], 0,
 function( L )
 local h;
 h:= NiceAlgebraMonomorphism( L );
 if h <> fail then
 return Dimension( Range( h ) ) < infinity;
 else
 TryNextMethod();
 fi;
end);

##############################################################################
##
## FpLieAlgebraEnumeration( <arg> ) Juergen Wisliceny
## Willem de Graaf
##
## This function calculates a homomorphism of a finitely presented
## Lie algebra onto a structure constants algebra.
## The algorithm is guaranteed to terminate when the algebra is finite
## dimensional. In full length the list <arg> contains `FpL', a
## finitely presented Lie algebra, `MAX_WEIGHT', a bound on the length
## of the monomials (used for nilpotent quotients), `weights' (a list
## of weights of the variables) and finally a boolean indicating
## whether the relations are homogeneous (if so then the nilpotent
## quotient will be graded, the grading is set as an attribute of the
## range of the homomorphism).
##
## By a straightforward application of the Jacobi identity (see also the
## comments to the sub-function `LeftNormalization'), it can be seen that
## the space of all commutators of degree `n' is spanned by all left
## normed commutators (i.e., commutators of the form [[[[a,b],c],d]...]).
## By antisymmetry we have that a and b can be chosen such that a > b.
## This is the format for elements of the free Lie algebra used in the
## program. A left-normed commutator is represented by a list
## `[a,b,c,d,..]', meaning `[[[[a,b],c],d]...]'. A monomial is such a list
## together with a coefficient, e.g., `[ [a,b,c,d..], -2/3 ]'. Finally, a
## polynomial is a list of monomials.

InstallGlobalFunction( FpLieAlgebraEnumeration,

 function( arg )

local ReductionModuloTable, #
 LeftNormalization, #
 SubsVarInRels, #
 CollectPolynomial, # Sub-functions.
 UpdateTable, #
 RemoveComm, #
 RemoveEntry, #
 SubstituteVariable, #
 Dcopy, #
 gradorder, #
 grado, #

 vg, # List of pairs (newly defined commutators)
 i,j,k,l,s, # Loop variables.
 end_reached, #
 table_init, # Booleans
 relation_found, #
 u,v,rr, # Lists of relations.
 r,u1,r1,r2,k1,k2,k11, # Polynomials, monomials etc.
 S,_T, # Structure constants tables.
 rowS, # A row of the multiplication table.
 sij,tij, # Entries of the multiplication table.
 inds, # Indices (list of integers).
 tab_pols, # List of polynomials of degree two.
 intrel, # Initial relations (after first conversion).
 pp, # Ext rep of a polynomial.
 cf, # Coefficient.
 t1,t2, # Indices.
 max, # Maximum.
 R, # Lists of commtators that have been defined.
 Rw1, # A new roe of `R'.
 one, # One of the field.
 zero, # Zero of the field.
 d, # maximum of the list of (pseudo-)generators
 e, # Flat(e) is the list of (pseudo-)generators
 w,ww, # Weights.
 temp,
 defs, # Definitions of generators in terms of other
 # generators.
 FL, # The free Lie algebra.
 rels, # Relators.
 bound, # Bound for `w'.
 gens, # Generators of `FpL'.
 imgs, # Images.
 map, # The map that is constructed.
 im, # An image.
 Fam, # Elements family of `FpL'.
 K, # Structure constants algebra.
 FpL,wts,wght,pos,fle,weight,MAX_WEIGHT,genweights,comp_grad,is_hom,
 bas,gradcomps,degs,bgc;

 FpL:= arg[1];
 if Length( arg ) >= 2 then
 MAX_WEIGHT:= arg[2];
 else
 MAX_WEIGHT:= infinity;
 fi;

 Fam:= ElementsFamily( FamilyObj( FpL ) );
 FL:= Fam!.freeAlgebra;
 rels:= Fam!.relators;

 if Length( arg ) >= 3 then
 genweights:= arg[3];
 else
 genweights:= List( GeneratorsOfAlgebra( FL ), x -> 1 );
 fi;

 if Length( arg ) = 4 then
 is_hom:= arg[4];
 else
 is_hom:= false;
 fi;

 bound:= infinity;

 _T:=[];
 one:= One( LeftActingDomain( FL ) );
 zero:= Zero( LeftActingDomain( FL ) );

# Some small functions.....

 Dcopy:= function( l )

 # Deep copying, also copying the holes...

 local m,i;

 if not IsList(l) then return ShallowCopy(l); fi;
 m:=[];
 for i in [1..Length(l)] do
 if IsBound(l[i]) then m[i]:= Dcopy(l[i]); fi;
 od;
 return m;
 end;

##############################################################
##############################################################
# v, w are associative monomials. is v>w?
##
## v > w if and only if 1) Length(v)>Length(w) or
## 2) Length(v)=Length(w) and Length(v) > 1 and
## v[2] > w[2] or
## 3) Length(v)=Length(w) and v[2]=w[2] and
## v>w alphabetically.
##

 gradorder:=function(v,w)
 local k,l;

 k:=Length(v[1]); l:=Length(w[1]);
 if k<>l then return k>l;
 elif k>1 and v[1][2]<>w[1][2] then return v[1][2]>w[1][2];
 else return v>w;
 fi;

 end;

 grado:= function( v, w )
 # tries to mimic gradorder for monomials of deg 2.

 if v[2]<>w[2] then return v[2]>w[2];
 else return v>w;
 fi;
 end;

########################################################################

 CollectPolynomial:= function( r )

 # A function that collects equal things together, and gets rid of
 # things in the polynomial r that are zero.

 local i,n,t;

 if r <> [ ] then

 # first regularize...
 for i in r do
 if Length(i[1])>1 and i[1][1]<i[1][2] then
 t:=i[1][2]; i[1][2]:=i[1][1]; i[1][1]:=t;
 i[2]:=-i[2];
 fi;
 od;

 Sort( r, gradorder );
 n:= Length( r );

 for i in [1..n-1] do
 if r[i][2]=0*r[i][2] or
 (Length(r[i][1])>1 and r[i][1][1]=r[i][1][2]) then

 #the thing is zero; get rid of it.

 Unbind(r[i]);
 elif r[i][1] = r[i+1][1] then

 #the monomials are equal; collect them together.

 r[i+1][2]:=r[i][2]+r[i+1][2];
 Unbind(r[i]);
 fi;
 od;
 if r[n][2]=0*r[n][2] or
 (Length(r[n][1])>1 and r[n][1][1]=r[n][1][2]) then

 #the thing is zero; get rid of it.

 Unbind(r[n]);
 fi;
 r:= Compacted( r );

 fi;

 return r;

 end;

 ReductionModuloTable := function( k )

 # In this function a Lie polynomial `k' in standard form is
 # reduced by one step modulo the commutators already known by the
 # table. So if [x_i,x_j]= c*z is a relation in the table, and `k'
 # contains a monomial of the form [ [i,j,k,....], cf ] then this
 # monomial is replaced by [ [z,k,...], -c*cf ]

 local i,j,k1,l,m,tst,t,s,cf,p,q,a;

 a:= Dcopy( k );
 for i in [1..Length(a)] do
 l:=Length(a[i][1]);
 if l>1 then
 s:= a[i][1][1]; t:=a[i][1][2];
 if s < t then
 p:= t; q:= s;
 else
 p:=s; q:=t;
 fi;
 if IsBound( _T[p] ) and IsBound( _T[p][q] ) then
 k1:= [ ];
 tst:= _T[p][q];
 if s <> p then cf:= -1;
 else cf:= 1;
 fi;
 for j in [1..Length(tst[1])] do
 Add( k1, [[tst[1][j]], cf*a[i][2]*tst[2][j]] );
 od;
 if l>2 then
 m:=a[i][1]{[3..l]};
 for j in [1..Length(k1)] do
 Append(k1[j][1],m);
 od;
 fi;
 Unbind(a[i]);
 Append(a,k1);
 fi;
 fi;
 od;
 a:=Compacted(a);
 a:= CollectPolynomial( a );

 if a = [ ] or a[1][2] = one then
 return a;

 else
 cf:= 1/a[1][2];
 return List( a, x -> [x[1],x[2]*cf] );
 fi;
 end;

 LeftNormalization:= function( rel )

 # a left-normed monomial is of the form
 #
 # [a,b,c,d,e,...], meaning [[[[[a,b],c],d],e],...]
 #
 # Using the Jacobi identity every commutator can be represented
 # as a linear combination of left-normed commutators.
 #
 # In this function a polynomial `rel' is left normed.
 # The Jacobi identity is applied successively to achieve this.
 # This means that an expression of the form
 #
 # [a,b,c,[d,e],f] (where a,b,c are generators (this part is already
 # `done') and [d,e] is any bracketed expression having d and e as
 # left and right subtrees,
 #
 # to a sum
 #
 # [a,b,c,d,e,f] - [a,b,c,e,d,f].
 #
 # Justification:
 #
 # [a,b,c,[d,e]]=[[[a,b],c],[d,e]] = [X,[d,e]] (with X=[[a,b],c])
 # =-[d,[e,X]]-[e,[X,d]]
 # =[[e,X],d]+[[X,d],e]
 # =-[[X,e],d]+[[X,d],e].
 #

 local i,j,s,s1,s2,t,step_occurred;

 step_occurred:= true;
 while step_occurred do

 # if there no longer occur any Jacobi steps, then we stop.

 i:=0;
 step_occurred:= false;
 while i < Length( rel ) do
 i:=i+1; j:=0;
 while j<Length(rel[i][1]) and not step_occurred do
 j:=j+1;
 if IsList(rel[i][1][j]) and Length(rel[i][1][j])=2 then
 step_occurred:= true;
 s:=rel[i][1]{[1..j-1]}; #i.e., the part already done (the X)
 s1:=Concatenation(s,rel[i][1][j]);
 s2:=Concatenation(s,[rel[i][1][j][2],rel[i][1][j][1]]);
 t:=rel[i][1]{[j+1..Length(rel[i][1])]};
 Append(s1,t); Append(s2,t);
 rel[i][1]:=Dcopy(s1);

 # If j=1 (so if the tree starts with [x,y], then we didn't do
 # much other than changing the notation ([[x,y],b] -> [x,y,b]).

 if j>1 then Add(rel,[s2,-rel[i][2]]); fi;
 fi;
 od;
 od;
 od;
 return rel;
 end;

 SubsVarInRels:= function( rels, rs )

 # Here `rs' is a relation of the form `var = othervars', and `rels' is
 # a list of Lie polynomials. This function substitutes `var'
 # everywhere in the polynomials `rels'.

 local i,j,p,s,s1,s2,result,rel,cf;

 result:= [ ];

 for rel in rels do

 i:= 1;
 while i <= Length(rel) do

 p:= Position( rel[i][1], rs[1][1][1] );
 if p <> fail then

 # s will be the polynomial that is gotten from r by substituting
 # `the rest of rs' for the variable rs[1][1][1] on the position p
 # in r.

 s:= Dcopy( rs{[2..Length(rs)]} );
 s1:= rel[i][1]{[1..p-1]};
 s2:= rel[i][1]{[p+1..Length(rel[i][1])]};
 for j in [1..Length(s)] do
 s[j][1]:=Concatenation(s1,s[j][1],s2);
 od;
 s:= List( s, x -> [ x[1], -rel[i][2]*x[2] ] );
 Append( rel, s );
 Unbind( rel[i] );
 rel:= Compacted( rel );

 else
 i:= i+1;
 fi;

 od;

 #collect the result...

 rel:= CollectPolynomial( rel );
 if rel <> [ ] and rel[1][2] <> one then
 cf:= 1/rel[1][2];
 rel:= List( rel, x -> [x[1],cf*x[2]] );
 fi;
 if rel <> [ ] then AddSet( result, rel ); fi;

 od;
 return result;
 end;

 UpdateTable:= function( i, j, p )

 # Sets the commutator [xi,xj] in the table equal to the polynomial `p'.

 local inds,cfs,k,s,t;

 inds:=[];
 cfs:=[];
 for k in [1..Length(p)] do
 inds[k]:= p[k][1][1];
 cfs[k]:= p[k][2];
 od;
 if i < j then
 s:= j; t:= i;
 else
 s:=i; t:= j;
 fi;

 if s = i then cfs:= -cfs; fi;
 if not IsBound(_T[s]) then _T[s]:=[]; fi;
 _T[s][t]:= [inds,cfs];

 end;

 RemoveEntry:= function( k )

 # Removes all occurrences of the variable xk in the commutators
 # of the table.

 local i;

 Unbind(_T[k]);
 for i in [1..Length(_T)] do
 if IsBound( _T[i] ) then Unbind(_T[i][k]); fi;
 od;
 end;

 RemoveComm:= function( k, l )

 # Removes the commutator [xk,xl] from the table.

 local s,t;
 s:= Maximum( k, l ); t:= Minimum( k, l );
 if IsBound(_T[s][t]) then Unbind(_T[s][t]); fi;
 end;

 SubstituteVariable:= function( coms, rel )

 # Here `rel' is a polynomial of the form `var = othervars'; this
 # function substitutes `var' for `othervar' in the commutators of
 # the table prescribed by `coms'.

 local var,inds,i,cfs,c,Tij,pos,cf,ii,cc,ind,s,t;

 var := rel[1][1][1];
 inds:= [ ]; cfs:= [ ];
 for i in [2..Length(rel)] do
 Add( inds, rel[i][1][1] );
 Add( cfs, -rel[i][2] );
 od;
 cfs:= cfs/rel[1][2];

 for c in coms do

 s:= Maximum( c ); t:= Minimum( c );
 Tij:= _T[s][t];
 pos:= Position( Tij[1], var );
 if pos <> fail then
 Remove( Tij[1], pos );
 cf:= Tij[2][pos];
 if s <> c[1] then cf:= -cf; fi;
 Remove( Tij[2], pos );
 Append( Tij[1], inds );
 Append( Tij[2], cf*cfs );
 ii:= [ ]; cc:= [ ];
 if Tij[1] <> [ ] then
 SortParallel( Tij[1], Tij[2] );
 ind:= Tij[1][1]; cf:= Tij[2][1];
 ii:= [ ]; cc:= [ ];
 for i in [2..Length(Tij[1])] do
 if Tij[1][i] = ind then
 cf:= Tij[2][i] + cf;
 else
 Add( ii, ind );
 Add( cc, cf );
 ind:= Tij[1][i];
 cf:= Tij[2][i];
 fi;
 od;
 Add( ii, ind ); Add( cc, cf );
 fi;
 _T[s][t]:= [ ii, cc ];
 fi;

 od;

 end;

 wght:= function( e, wts, var )

 local p,q;

 p:= PositionProperty( e, x -> var in x );
 q:= Position( e[p], var );

 return wts[p][q];
 end;

##############################################################################
#
# The program starts. First the relations are transformed into internal format.
# That is: represented as lists of lists etc., and left-normalized.
#

 # `intrel' will be the set of relations, but represented in
 # `internal form'; meaning [ [ [[1,2],3], 1 ], [[4],-1] ], instead
 # of (x1*x2)*x3-x4 etc.

 intrel:= [ ];
 for r in rels do
 pp:= Dcopy( ExtRepOfObj( r )[2] );
 Add( intrel, List( [1,3..Length(pp)-1], x -> [ pp[x], pp[x+1] ] ) );
 od;

#############################################################################
 # now we left normalize the relations, using `LeftNormalization', i.e.,
 # the relations are written as [ [ [1,2,5], -1], [.....], [......],.... ]
 # furthermore, all relations of degree at most two will go into `pr'
 # (those will be used to initialize the table). All the others go into
 # `u'.

 tab_pols:= [ ]; u:= [ ];

 for r in intrel do
 max:= 0;
 for j in [1..Length(r)] do
 if not IsList(r[j][1]) then #transform [ i, cst ] into [ [i], cst ]
 r[j][1]:= [ r[j][1] ];
 fi;
 if Length(Flat(r[j][1])) >= 2 and r[j][1][1]=r[j][1][2] then
 Unbind(r[j]);
 else
 max:= Maximum( max, Length( Flat(r[j][1]) ) );
 fi;
 od;
 r:= Compacted( r );
 r:= LeftNormalization( r );
 r:= CollectPolynomial( r );
 if not max = 0 then
 if max <= 2 then
 cf:= 1/r[1][2];
 r:= List( r, x -> [x[1],cf*x[2]] );
 Add( tab_pols, r); # So if the relation only
 # involves monomials of deg
 # at most two, then this relation
 # goes into the 'tab_pols'.
 else
 Add( u, r );
 fi;
 fi;
 od;

 e:= [ List( [1..Length( GeneratorsOfAlgebra( FL ) )], x -> x ), [ ] ];
 if MAX_WEIGHT < infinity then
 wts:= [ genweights, [] ];
 comp_grad:= true;
 else
 wts:= [ ];
 comp_grad:= false;
 fi;

 if e[1] = [ ] then
 K:= LieAlgebraByStructureConstants( LeftActingDomain( FL ),
 EmptySCTable( 0, zero, "antisymmetric" ) );
 gens:= GeneratorsOfAlgebra( FpL );
 imgs:= List( gens, x -> Zero( K ) );
 map:= Objectify( TypeOfDefaultGeneralMapping( FpL, K,
 IsSPGeneralMapping
 and IsAlgebraGeneralMapping
 and IsFptoSCAMorphism
 and IsAlgebraGeneralMappingByImagesDefaultRep ),
 rec(
 generators := gens,
 genimages := imgs
 ) );
 SetMappingGeneratorsImages(map,[Immutable(gens),Immutable(imgs)]);
 return map;
 fi;

 # `v' will be a history of relations, i.e., `v[w]' will be the relations
 # as they were when the program was dealing with weight `w'. This is
 # used to reset the relations if a collision among variables is found.

 v:= [ Dcopy( u ) ];
 d:= Maximum( e[1] );
 R:= [ [], [] ];
 end_reached:= false;
 w:= 0;
 defs:= [ ];

 while w < bound do

 table_init:= false;

 while not table_init do

#######################################################################
# Initialize the table....
# Meaning: fill in all possible commutators of generators using the
# relations, make definitions for the commutators that cannot be decided
# upon by using the relations. If this leads to a relation among the variables,
# then that relation is substituted first, and the process is started all
# over again.

 relation_found:= false;

 for r in tab_pols do

 r1:= ReductionModuloTable( r );
 if r1<>[] then

 if Length(r1[1][1])=1 then
 relation_found:= true;
 break;
 else
 for k in [2..Length(r1)] do
 if Length(r1[k][1])=2 then
 d:=d+1;
 Add(e[2],d);
 if comp_grad then
 Add(wts[2],wght(e,wts,r1[k][1][1])+
 wght(e,wts,r1[k][1][2]) );
 fi;
 UpdateTable( r1[k][1][1], r1[k][1][2], [ [[d],-one] ] );
 Add(R[2],r1[k][1]);
 r1[k][1]:=[d];
 fi;
 od;
 UpdateTable( r1[1][1][1], r1[1][1][2], r1{[2..Length(r1)]} );
 fi;
 Add(R[2],r1[1][1]);

 fi;
 od;

 if not relation_found then

 # i.e., the previous loop has been executed without breaking
 # caused by finding a relation among the generators.

 vg:=Difference( List( Combinations(e[1],2), Reversed ), R[2] );
 Append( R[2], vg );
 for i in [1..Length(vg)] do
 d:=d+1;
 Add(e[2],d);
 if comp_grad then
 Add(wts[2],wght(e,wts,vg[i][1])+wght(e,wts,vg[i][2]));
 fi;
 UpdateTable( vg[i][1], vg[i][2], [ [[d],-one] ] );
 od;

 rr:= [ ];
 for i in [1..Length(u)] do
 r:= Dcopy( u[i] );
 while true do
 r1:= ReductionModuloTable( r );
 if r1 = r then break;
 else r:= r1;
 fi;
 od;
 if r <> [ ] then
 if Length(r[1][1]) = 1 and r[1][1][1] in e[1] then
 relation_found:= true;
 break;
 else
 Add( rr, r );
 fi;
 fi;
 od;
 fi;

 if relation_found then

 # i.e., a relation among the variables has been found in the
 # previous piece of code.

 w:= Position( List( e, x -> r1[1][1][1] in x ), true );
 if w = 1 then
 if comp_grad then
 pos:= Position( e[1], r1[1][1][1] );
 Remove( wts[1], pos );
 Remove( e[1], pos);
 else
 RemoveSet( e[1], r1[1][1][1] );
 fi;

 Add( defs, r1 );
 if e[1] = [ ] then
 K:= LieAlgebraByStructureConstants( LeftActingDomain( FL ),
 EmptySCTable( 0, zero, "antisymmetric" ) );
 gens:= GeneratorsOfAlgebra( FpL );
 imgs:= List( gens, x -> Zero( K ) );
 map:= Objectify( TypeOfDefaultGeneralMapping( FpL, K,
 IsSPGeneralMapping
 and IsAlgebraGeneralMapping
 and IsFptoSCAMorphism
 and IsAlgebraGeneralMappingByImagesDefaultRep ),
 rec(
 generators := gens,
 genimages := imgs
 ) );
 SetMappingGeneratorsImages(map,[Immutable(gens),Immutable(imgs)]);
 return map;
 fi;
 e[2]:= [ ];
 if comp_grad then
 wts[2]:= [ ];
 fi;
 tab_pols:= SubsVarInRels( tab_pols, r1 );
 u:= SubsVarInRels( u, r1 );
 _T:= [ ];
 R:= [ [], [] ];
 else
 if comp_grad then
 pos:= Position( e[w], r1[1][1][1] );
 Remove( wts[w], pos );
 Remove( e[w], pos);
 else
 RemoveSet( e[w], r1[1][1][1]);
 fi;

 u:= SubsVarInRels( v[w-1], r1 );
 SubstituteVariable( R[w], r1 );
 fi;

 else
 u:= rr;
 table_init:= true;
 fi;

 od;

##########################################################################
#
# The table has been initialized, and the commutators of weight 2
# have been defined. Now the process of increasing the weight starts.
#

 w:=1;
 while w < bound do

 w:=w+1;
 Sort( R[w], grado );

 if comp_grad then
 fle:= Flat(e);
 for i in [1..Length(fle)] do
 for j in [i+1..Length(fle)] do
 if wght(e,wts,fle[i])+wght(e,wts,fle[j])>MAX_WEIGHT then
 UpdateTable( fle[i], fle[j], [] );
 fi;
 od;
 od;
 fi;

#############################################################################
# reduction modulo relations and Jacobi identity....
#
# In this function also _T is changed; but if the function
# exits with a relation among the vars, then we change `_T' back to its
# old value (the copy `S').
#

 S:= Dcopy( _T );
 rr:= Dcopy( u );
 Rw1:= [ ];
 e[w+1]:= [ ];
 if comp_grad then
 wts[w+1]:= [ ];
 fi;
 d:= Maximum( Flat( e ) );
 relation_found:= false;

 for r in R[w] do

 t1:=r[1]; t2:=r[2];
 if t1 > t2 then
 tij:= _T[t1][t2];
 else
 tij:= ShallowCopy( _T[t2][t1] );
 tij[2]:= -ShallowCopy( tij[2] );
 fi;
 r1:= List( [1..Length(tij[1])], k -> [ [tij[1][k]], tij[2][k] ] );

 for j in e[1] do

 # The Jacobi identity that will be inspected reads as
 # [ [ t1, t2 ], j ] - [ [ t1, j ], t2 ] + [ [ t2, j ], t1 ] = 0
 # This relation can be evaluated (using the partial table) to a
 # polynomial of degree <=2. This will lead to new definitions
 # (in the case of deg. = 2), or collisions (in the case of
 # deg. = 1).

 if t2 > j then

 if t1 > j then
 tij:= _T[t1][j];
 else
 tij:= ShallowCopy( _T[j][t1] );
 tij[2]:= -ShallowCopy( tij[2] );
 fi;
 k1:= List( [1..Length(tij[1])], i->[ [tij[1][i],t2],-tij[2][i] ]);
 if t2 > j then
 tij:= _T[t2][j];
 else
 tij:= ShallowCopy( _T[j][t2] );
 tij[2]:= -ShallowCopy( tij[2] );
 fi;
 k2:= List( [1..Length(tij[1])], i->[ [tij[1][i],t1],tij[2][i] ]);

 r2:= Dcopy(r1);
 for i in r2 do Add( i[1], j ); od;

 k:= Concatenation( k1, k2, r2 );
 k:= CollectPolynomial( k );
 k:= ReductionModuloTable( k );

 if k <> [ ] then

 # Produce a relation of the form a = c1*var1+c2*var2...
 # by making new definitions. (Where a is either a commutator
 # or a variable).

 i:= 2;
 while i <= Length( k ) do

 if Length(k[i][1]) = 2 then
 if comp_grad then
 weight:= wght(e,wts,k[i][1][1])+ wght(e,wts,k[i][1][2]);
 else
 weight:= 0;
 fi;

 if weight <= MAX_WEIGHT then
 if comp_grad then
 Add( wts[w+1], weight );
 fi;
 d:= d+1;
 Add( e[w+1], d );
 UpdateTable( k[i][1][1], k[i][1][2], [ [[d],-one] ] );
 Add( Rw1, k[i][1] );
 k[i][1]:= [ d ];
 else
 Remove( k, i );
 fi;
 fi;
 i:= i+1;
 od;
 k11:= k[1][1];

 if Length(k11) = 2 then

 # The `a' in the comment above is a commutator; hence a new
 # entry for the table has been found.

 UpdateTable( k11[1], k11[2], k{[2..Length(k)]} );
 Add( Rw1, k11 );
 elif Length(k11) = 1 then

 ww:= 0;
 for i in [1..Length(e)] do
 if k11[1] in e[i] then ww:=i; break; fi;
 od;

 if ww = w+1 then

 # A collision (among the new basis elements) has been found.

 if comp_grad then
 pos:= Position( e[w+1], k11[1] );
 Remove( wts[w+1], pos );
 fi;
 RemoveSet( e[w+1], k11[1] );
 RemoveEntry( k11[1] );
 SubstituteVariable( Rw1, k );
 rr:= SubsVarInRels( rr, k );
 elif ww > 0 then
 _T:=S;
 relation_found:= true;
 r1:= [ ww, k ];
 break;
 fi;
 fi;

 i:= 0;
 while i < Length(rr) do

 i:= i+1;
 # Reduce the relations modulo the table and process them.

 while true do
 u1:= ReductionModuloTable( rr[i] );
 if u1 = rr[i] then break;
 else rr[i]:=u1;
 fi;
 od;

 if rr[i]=[] then
 Unbind( rr[i] );
 elif Length(rr[i][1][1])=1 then

 ww:= 0;
 temp:= rr[i][1][1][1];
 for l in [1..Length(e)] do
 if temp in e[l] then ww:=l; break; fi;
 od;

 if ww = w+1 then
 if comp_grad then
 pos:= Position( e[w+1],rr[i][1][1][1] );
 Remove( wts[w+1], pos );
 fi;
 RemoveSet(e[w+1],rr[i][1][1][1]);
 RemoveEntry(rr[i][1][1][1]);
 SubstituteVariable( Rw1, rr[i] );
 temp:= rr[i];
 Remove( rr, i );
 rr:= SubsVarInRels( rr, temp );
 i:= i-1; # (last call removed holes...).
 elif ww > 0 then
 _T:=S;
 relation_found:= true;
 r1:= [ww,rr[i]];
 break;
 fi;
 elif Length(rr[i][1][1])=2 then
 max := 0;
 for s in rr[i] do
 ww:= 0;
 for l in [1..Length(e)] do
 if s[1][1] in e[l] then ww:=l; break; fi;
 od;
 if Length(s[1]) = 1 then
 max:= Maximum(max,ww);
 else
 # We calculate the weight of `s[1][1]' + the weight
 # of `s[1][2]' i.e., the weight of `[s[1][1], s[1][2]]'

 for l in [1..Length(e)] do
 if s[1][2] in e[l] then ww:=ww+l; break; fi;
 od;
 max:= Maximum(max,ww);
 fi;
 od;

 if max = w+1 then
 s:= 2;
 while s <= Length( rr[i] ) do
 if Length(rr[i][s][1]) = 2 then
 if wts <> [ ] then
 weight:= wght(e,wts,rr[i][s][1][1])+
 wght(e,wts,rr[i][s][1][2] );
 else
 weight:= 0;
 fi;
 if weight <= MAX_WEIGHT then
 d:= d+1;
 Add( e[w+1], d );
 if wts <> [ ] then
 Add( wts[w+1], weight );
 fi;
 UpdateTable( rr[i][s][1][1], rr[i][s][1][2],
 [ [[d], -one] ] );
 Add(Rw1,rr[i][s][1]);
 rr[i][s][1]:= [ d ];
 else
 Remove( rr[i], s );
 fi;
 fi;
 s:= s+1;
 od;
 Add(Rw1,rr[i][1][1]);
 UpdateTable( rr[i][1][1][1], rr[i][1][1][2],
 rr[i]{[2..Length(rr[i])]});
 Unbind(rr[i]);
 fi;
 fi;

 od;
 if relation_found then break; fi;
 rr:=Compacted(rr);
 fi;

 fi;
 od;
 if relation_found then break; fi;
 od;

##########################################################################

 if relation_found then

 # Here `r1[2]' is a relation among basis elements.
 # `r1[1]' is the weight of the homogeneous component containing
 # the first variable (variable of highest weight).

 w:= r1[1];

 if w = 1 then

 # A relation among the variables of weight 1 has been found.
 # We reset everything and return to the point where the table
 # is initialized.

 if comp_grad then
 pos:= Position( e[1], r1[2][1][1][1] );
 Remove( wts[1], pos );
 fi;
 RemoveSet( e[1], r1[2][1][1][1] );
 Add( defs, r1[2] );
 if e[1]=[] then
 K:= LieAlgebraByStructureConstants( LeftActingDomain( FL ),
 EmptySCTable( 0, zero, "antisymmetric" ) );
 gens:= GeneratorsOfAlgebra( FpL );
 imgs:= List( gens, x -> Zero( K ) );
 map:= Objectify( TypeOfDefaultGeneralMapping( FpL, K,
 IsSPGeneralMapping
 and IsAlgebraGeneralMapping
 and IsFptoSCAMorphism
 and IsAlgebraGeneralMappingByImagesDefaultRep ),
 rec(
 generators := gens,
 genimages := imgs
 ) );
 SetMappingGeneratorsImages(map,[Immutable(gens),Immutable(imgs)]);
 return map;
 fi;
 e:=[ e[1], [] ];
 if comp_grad then
 wts:= [ wts[1], [] ];
 fi;
 u:= SubsVarInRels( v[1], r1[2] );
 tab_pols:= SubsVarInRels( tab_pols, r1[2] );
 _T:=[];
 R:=[ [ ], List( tab_pols, x -> x[1][1] ) ];
 v[1]:= Dcopy( u );

 # We break to the principal loop.
 break;
 else

# here `r1[2]' is of the form `var=something' where `var' is of weight
# `w', and `w>1'. This means that `var' was introduced somewhere; namely
# on level `w'. Hence the definition was [x_i,x_j]=var, where w(x_i)+
# w(x_j)=w. Hence `var' only appears in tails (right hand sides) of commutators
# of weight `>= w'. Now `var' is substituted in all products of weight `w',
# and the program starts again on that level.

 if comp_grad then
 pos:= Position( e[w], r1[2][1][1][1] );
 Remove( wts[w], pos );
 fi;
 RemoveSet(e[w], r1[2][1][1][1]);
 u:= SubsVarInRels( v[w-1], r1[2]);
 v[w-1]:=u;
 w:= w-1;
 SubstituteVariable( R[w+1], r1[2] );
 for i in [w+2..Length(e)] do e[i]:=[]; od;
 for i in [w+2..Length(R)] do
 for j in [1..Length(R[i])] do
 RemoveComm( R[i][j][1], R[i][j][2] );
 od;
 R[i]:= [ ];
 od;

 fi;

 else

 # Here Jacobi identities have been applied
 # without finding collisions between variables.

 if e[w] = [ ] and not end_reached then
 bound:= 2*w; end_reached:= true;
 elif w = bound and not end_reached then
 return fail;
 fi;

 R[w+1]:= Rw1; v[w]:= rr; u:= rr;

 if Flat( e ) = [ ] then
 K:= LieAlgebraByStructureConstants( LeftActingDomain( FL ),
 EmptySCTable( 0, zero, "antisymmetric" ) );
 gens:= GeneratorsOfAlgebra( FpL );
 imgs:= List( gens, x -> Zero( K ) );
 map:= Objectify( TypeOfDefaultGeneralMapping( FpL, K,
 IsSPGeneralMapping
 and IsAlgebraGeneralMapping
 and IsFptoSCAMorphism
 and IsAlgebraGeneralMappingByImagesDefaultRep ),
 rec(
 generators := gens,
 genimages := imgs
 ) );
 SetMappingGeneratorsImages(map,[Immutable(gens),Immutable(imgs)]);
 return map;
 fi;

 d:= Maximum( Flat( e ) );
 vg:= [ ];
 for i in e[w] do
 for j in e[1] do
 if i>j then AddSet( vg, [i,j] ); fi;
 od;
 od;
 vg:= Difference( vg, R[w+1] );
 Append( R[w+1], vg );

 for i in [1..Length(vg)] do
 if comp_grad then
 weight:= wght(e,wts,vg[i][1])+wght(e,wts,vg[i][2]);
 else
 weight:= 0;
 fi;
 if weight <= MAX_WEIGHT then
 d:= d+1;
 Add( e[w+1], d );
 if comp_grad then
 Add( wts[w+1], weight );
 fi;
 Add( R[w+1], vg[i] );
 UpdateTable( vg[i][1], vg[i][2], [ [[d],-1*one] ]);
 else
 UpdateTable( vg[i][1], vg[i][2], [ ] );
 fi;
 od;

 fi;

 od; # end of the loop in which `w' is successively increased.
 od; # end of the main loop,

 # Now we construct a table of structure constants from `_T'.

 e:=Filtered(e,x->x<>[]);
 inds:=Flat(e);

 S:=[];
 for i in inds do
 rowS:= [ ];
 for j in inds do
 if i=j then
 Add( rowS, [ [], [] ] );
 else
 if i < j then
 tij:= ShallowCopy( _T[j][i] );
 tij[2]:= -ShallowCopy( tij[2] );
 else
 tij:= _T[i][j];
 fi;
 sij:=[[],[]];
 for k in [1..Length(tij[1])] do
 sij[1][k]:= Position( inds, tij[1][k] );
 sij[2][k]:= tij[2][k];
 od;
 Add( rowS, sij );
 fi;
 od;
 Add( S, rowS );
 od;
 Add( S, -1 ); Add( S, zero );

 K:= LieAlgebraByStructureConstants( LeftActingDomain( FL ), S );
 if is_hom then
 wts:= Flat( wts );
 bas:= [1..Dimension(K)];
 SortParallel( wts, bas );

 gradcomps:= [ ];
 degs:= [ ];
 k:= 1;
 while k <= Length(wts) do
 bgc:= [ Basis( K )[bas[k]] ];
 Add( degs, wts[k] );
 while k < Length( wts ) and wts[k]=wts[k+1] do
 k:= k+1;
 Add( bgc, Basis(K)[bas[k]] );
 od;
 Add( gradcomps, VectorSpace( LeftActingDomain( K ), bgc ) );
 k:= k+1;
 od;

 Add( gradcomps, Subspace( K, [ ] ) );

 SetGrading( K, rec( min_degree:= Minimum( wts ),
 max_degree:= Maximum( wts ),
 source:= Integers,
 hom_components:= function( d )
 if d in degs then
 return gradcomps[Position(degs,d)];
 else
 return Last(gradcomps);
 fi;
 end
 ) );

 fi;

 gens:= GeneratorsOfAlgebra( FpL );

 if Dimension( K ) = 0 then #trivial case
 imgs:= List( gens, x -> Zero(K) );
 else
 # We process the definitions, (of generators as linear combinations
 # of other generators).
 i:= Length( defs );
 while i > 1 do
 for j in [1..i-1] do
 for k in [1..Length(defs[j])] do
 if defs[j][k][1] = defs[i][1][1] then
 Append( defs[j], List( defs[i]{[2..Length(defs[i])]}, x ->
 [ x[1], -defs[j][k][2]*x[2] ] )
 );
 Unbind( defs[j][k] );
 fi;
 od;
 defs[j]:= Compacted( defs[j] );
 od;
 i:= i-1;
 od;

 imgs:= [ ];

 #For every generator of the Fp Lie algebra we calculate an image...

 for i in [1..Length(gens)] do
 if i in e[1] then
 Add( imgs, Basis( K )[Position( inds, i )] );
 else
 for j in [1..Length(defs)] do
 if defs[j][1][1][1] = i then break; fi;
 od;
 im:= Zero( K );
 for k in [2..Length(defs[j])] do
 im:= im + -defs[j][k][2]*Basis( K )[defs[j][k][1][1]];
 od;
 Add (imgs, im );
 fi;
 od;
 fi;

# Construct the map...

 map:= Objectify( TypeOfDefaultGeneralMapping( FpL, K,
 IsSPGeneralMapping
 and IsAlgebraGeneralMapping
 and IsFptoSCAMorphism
 and IsAlgebraGeneralMappingByImagesDefaultRep ),
 rec(
 generators := gens,
 genimages := imgs
 ) );
 SetMappingGeneratorsImages(map,[Immutable(gens),Immutable(imgs)]);

 return map;

end );

InstallMethod( NiceAlgebraMonomorphism,
 "for a f.p. Lie algebra",
 true,
 [ IsLieAlgebra and IsSubalgebraFpAlgebra], 0,

 function( FpL )

 return FpLieAlgebraEnumeration( FpL );

end );

InstallGlobalFunction( NilpotentQuotientOfFpLieAlgebra,

 function( arg )

 local FpL,L,weights,is_homogeneous,rels,weight,w,r,k,j,er,fol,maxw,N;

 # unwrapping the arguments...

 if Length( arg ) = 2 then
 FpL:= arg[1]; maxw:= arg[2];
 L:= ElementsFamily( FamilyObj( FpL ) )!.freeAlgebra;
 weights:= List( GeneratorsOfAlgebra( L ), x -> 1 );
 elif Length( arg ) = 3 then
 FpL:= arg[1]; maxw:= arg[2]; weights:= arg[3];
 else
 Error("Number of arguments must be two or three");
 fi;

 # checking whether the relations are homogeneous; if so then
 # the resulting structure constants Lie algebra will have a
 # natural grading.

 is_homogeneous:= true;
 rels:= ElementsFamily( FamilyObj( FpL ) )!.relators;

 for r in rels do

 weight:= infinity;
 er:= ExtRepOfObj( r )[2];
 for k in [1,3..Length(er)-1] do
 fol:= Flat( [ er[k] ] );
 w:= 0;
 for j in fol do
 w:= w+weights[j];
 od;
 if weight = infinity then
 weight:= w;
 elif weight <> w then
 is_homogeneous:= false;
 break;
 fi;
 od;
 if not is_homogeneous then break; fi;
 od;

 N:= FpLieAlgebraEnumeration( FpL, maxw, weights, is_homogeneous );
 SetIsLieNilpotent( Range(N), true );
 return N;

end );

##############################################################################
##
#F FpLieAlgebraByCartanMatrix( <C> )
##
##
InstallGlobalFunction( FpLieAlgebraByCartanMatrix, function( C )

 local i,j,k, # Loop variables.
 l, # The rank.
 L, # The free Lie algebra.
 g, # Generators of `L'.
 x,h,y, # Lists of generators of `L'.
 rels, # List of relations.
 rx,ry; # Relations.

 l:= Length( C );
 L:= FreeLieAlgebra( Rationals, 3*l );
 g:= GeneratorsOfAlgebra( L );
 x:= g{[1..l]};
 h:= g{[l+1..2*l]};
 y:= g{[2*l+1..3*l]};

 rels:= [ ];
 for i in [1..l] do
 for j in [i+1..l] do
 Add( rels, h[i]*h[j] );
 od;
 od;

 for i in [1..l] do
 for j in [1..l] do
 if i=j then
 Add( rels, x[i]*y[j]-h[i] );
 else
 Add( rels, x[i]*y[j] );
 fi;
 od;
 od;

 for i in [1..l] do
 for j in [1..l] do
 Add( rels, h[i]*x[j]-C[j][i]*x[j] );
 Add( rels, h[i]*y[j]+C[j][i]*y[j] );
 od;
 od;

 for i in [1..l] do
 for j in [1..l] do
 if i <> j then
 rx:= x[j]; ry:= y[j];
 for k in [1..1-C[j][i]] do
 rx:= x[i]*rx; ry:= y[i]*ry;
 od;
 Add( rels, rx ); Add( rels, ry );
 fi;
 od;
 od;

 return L/rels;

end );

#############################################################################
##
#M JenningsLieAlgebra( <G> )
##
## The Jennings Lie algebra of the p-group G.
##
##

InstallMethod( JenningsLieAlgebra,
 "for a p-group",
 true,
 [IsGroup], 0,

function ( G )

 local J, # Jennings series of G
 Homs, # Homomorphisms of J[i] onto the quotient J[i]/J[i+1]
 grades, # List of the full images of the maps in Homs
 gens, # List of the generators of the quotients J[i]/J[i+1],
 # i.e., a basis of the Lie algebra.
 pos, # list of positions: if pos[j] = p, then the element
 # gens[j] belongs to grades[p]
 i,j,k, # loop variables
 tempgens,
 t, # integer
 T, # multiplication table of the Lie algebra
 dim, # dimension of the Lie algebra
 a,b,c,f, # group elements
 e, # ext rep of a group element
 co, # entry of the multiplication table
 p, # the prime of G
 F, # ground field
 L, # the Lie algebra to be constructed
 pimgs, # pth-power images
 B, # Basis of L
 vv, x, # elements of L
 comp, # homogeneous component
 grading, # list of homogeneous components
 pcgps, # list of pc groups, isom to the elts of `grades'.
 hom_pcg, # list of isomomorphisms of `grades[i]' to `pcgps[i]'.
 enum_gens, # List of numbers of elts of `gens' in extrep.
 pp, # Position in a list.
 hm;

 # We do not know the characteristic if `G' is trivial.
 if IsTrivial( G ) then
 Error( "<G> must be a nontrivial p-group" );
 fi;

 # Construct the homogeneous components of `L':

 J:=JenningsSeries ( G );
 Homs:= List ( [1..Length(J)-1] , x ->
 NaturalHomomorphismByNormalSubgroupNC( J[x], J[x+1] ));
 grades := List ( Homs , Range );
 hom_pcg:= List( grades, IsomorphismSpecialPcGroup );
 pcgps:= List( hom_pcg, Range );
 gens := [];
 enum_gens:= [ ];
 pos := [];
 for i in [1.. Length(grades)] do
 tempgens:= GeneratorsOfGroup( pcgps[i] );
 Append ( gens , tempgens);

 # Record the number that each generator has in extrep.
 Add( enum_gens, List( tempgens, x -> ExtRepOfObj( x )[1] ) );
 Append ( pos , List ( tempgens , x-> i ) );
 od;

 # Construct the field and the multiplication table:

 dim:= Length(gens);
 p:= PrimePGroup( G );
 F:= GF( p );
 T:= EmptySCTable( dim , Zero(F) , "antisymmetric" );
 pimgs := [];
 for i in [1..dim] do
 a:= PreImagesRepresentative( Homs[pos[i]] ,
 PreImagesRepresentative( hom_pcg[pos[i]], gens[i] ) );

 # calculate the p-th power image of `a':

 if pos[i]*p <= Length(Homs) then
 Add( pimgs, Image( hom_pcg[pos[i]*p],
 Image( Homs[pos[i]*p], a^p) ) );
 else
 Add( pimgs, "zero" );
 fi;

 for j in [i+1.. dim] do
 if pos[i]+pos[j] <= Length( Homs ) then

 # Calculate the commutator [a,b], and map the result into
 # the correct homogeneous component.

 b:= PreImagesRepresentative( Homs[pos[j]],
 PreImagesRepresentative( hom_pcg[pos[j]], gens[j] ));
 c:= Image( hom_pcg[pos[i] + pos[j]],
 Image(Homs[pos[i] + pos[j]], a^-1*b^-1*a*b) );
 e:= ExtRepOfObj(c);
 co:=[];
 for k in [1,3..Length(e)-1] do
 pp:= Position( enum_gens[pos[i]+pos[j]], e[k] );
 t:= Sum( enum_gens{[1..pos[i]+pos[j]-1]}, Length )+pp;
 Add( co, One( F )*e[k+1] );
 Add( co, t );
 od;
 SetEntrySCTable( T, i, j, co );
 fi;

 od;
 od;

 L:= LieAlgebraByStructureConstants( F, T );

 B:= Basis( L );

 # Now we compute the natural grading of `L'.
 grading:= [ ];
 k:= 1;
 for i in [1..Length(enum_gens)] do
 comp:= [ ];
 for j in [1..Length(enum_gens[i])] do
 Add( comp, B[k] );
 k:= k+1;
 od;
 Add( grading, Subspace( L, comp ) );
 od;

 Add( grading, Subspace( L, [ ] ) );

 SetGrading( L, rec( min_degree:= 1,
 max_degree:= Length( grading ) - 1,
 source:= Integers,
 hom_components:= function( d )
 if d in [1..Length(grading)] then
 return grading[d];
 else
 return Last(grading);
 fi;
 end
 )
 );

 vv:= BasisVectors( B );

 # Set the pth-power images of the basis elements of `B':

 for i in [1..Length(pimgs)] do
 if pimgs[i] = "zero" then
 pimgs[i]:= Zero( L );
 else
 e:= ExtRepOfObj( pimgs[i] );
 x:= Zero( L );
 for k in [1,3..Length(e)-1] do
 pp:= Position( enum_gens[pos[i]*p], e[k] );
 t:= Sum( enum_gens{[1..pos[i]*p-1]}, Length )+pp;
 x:= x+ One( F )*e[k+1]*vv[t];
 od;
 pimgs[i]:= x;
 fi;
 od;
 SetPthPowerImages( B, pimgs );
 SetIsRestrictedLieAlgebra( L, true );
 FamilyObj(Representative(L))!.pMapping := pimgs;
 SetIsLieNilpotent( L, true );

 hm:= function( g, i )

 local h, e, x, k, pp, f, t;

 if not g in J[i] then
 Error("<g> is not an element of the i-th term of the series used to define <L>");
 fi;

 h:= Image( hom_pcg[i], Image(Homs[i], g ));
 e:= ExtRepOfObj(h);
 x:= Zero(L);
 for k in [1,3..Length(e)-1] do
 pp:= Position( enum_gens[i], e[k] );
 f:= GeneratorsOfGroup( pcgps[i] )[pp];
 t:= Position( gens, f );
 x:= x + e[k+1]*Basis(L)[t];
 od;
 return x;
 end ;

 SetNaturalHomomorphismOfLieAlgebraFromNilpotentGroup( L, hm );

 return L;

end );

#############################################################################
##
#M PCentralLieAlgebra( <G> )
##
## The p-central Lie algebra of the p-group G.
##
##
InstallMethod( PCentralLieAlgebra,
 "for a p-group",
 true,
 [IsGroup], 0,

function ( G )

 local J, # p-central series of G
 Homs, # Homomorphisms of J[i] onto the quotient J[i]/J[i+1]
 grades, # List of the full images of the maps in Homs
 gens, # List of the generators of the quotients J[i]/J[i+1],
 # i.e., a basis of the Lie algebra.
 pos, # list of positions: if pos[j] = p, then the element
 # gens[j] belongs to grades[p]
 i,j,k, # loop variables
 tempgens,
 t, # integer
 T, # multiplication table of the Lie algebra
 dim, # dimension of the Lie algebra
 a,b,c,f, # group elements
 e, # ext rep of a group element
 co, # entry of the multiplication table
 p, # the prime of G
 F, # ground field
 L, # the Lie algebra to be constructed
 B, # Basis of L
 vv, x, # elements of L
 comp, # homogeneous component
 grading, # list of homogeneous components
 pcgps, # list of pc groups, isom to the elts of `grades'.
 hom_pcg, # list of isomomorphisms of `grades[i]' to `pcgps[i]'.
 enum_gens, # List of numbers of elts of `gens' in extrep.
 pp, # Position in a list.
 pimgs, # pth power images
 hm;

 # We do not know the characteristic if `G' is trivial.
 if IsTrivial( G ) then
 Error( "<G> must be a nontrivial p-group" );
 fi;

 # Construct the homogeneous components of `L':

 p:= PrimePGroup( G );
 J:= PCentralSeries( G, p );
 Homs:= List ( [1..Length(J)-1] , x ->
 NaturalHomomorphismByNormalSubgroupNC( J[x], J[x+1] ));
 grades := List ( Homs , Range );
 hom_pcg:= List( grades, IsomorphismSpecialPcGroup );
 pcgps:= List( hom_pcg, Range );
 gens := [];
 enum_gens:= [ ];
 pos := [];
 for i in [1.. Length(grades)] do
 tempgens:= GeneratorsOfGroup( pcgps[i] );
 Append ( gens , tempgens);

 # Record the number that each generator has in extrep.
 Add( enum_gens, List( tempgens, x -> ExtRepOfObj( x )[1] ) );
 Append ( pos , List ( tempgens , x-> i ) );
 od;

 # Construct the field and the multiplication table:

 dim:= Length(gens);
 F:= GF( p );
 T:= EmptySCTable( dim , Zero(F) , "antisymmetric" );
 pimgs := [];
 for i in [1..dim] do
 a:= PreImagesRepresentative( Homs[pos[i]] ,
 PreImagesRepresentative( hom_pcg[pos[i]], gens[i] ) );

 # calculate the p-th power image of `a':

 if pos[i]+1 <= Length(Homs) then
 Add( pimgs, Image( hom_pcg[pos[i]+1],
 Image( Homs[pos[i]+1], a^p) ) );
 else
 Add( pimgs, "zero" );
 fi;

 for j in [i+1.. dim] do
 if pos[i]+pos[j] <= Length( Homs ) then

 # Calculate the commutator [a,b], and map the result into
 # the correct homogeneous component.

 b:= PreImagesRepresentative( Homs[pos[j]],
 PreImagesRepresentative( hom_pcg[pos[j]], gens[j] ));
 c:= Image( hom_pcg[pos[i] + pos[j]],
 Image(Homs[pos[i] + pos[j]], a^-1*b^-1*a*b) );
 e:= ExtRepOfObj(c);
 co:=[];
 for k in [1,3..Length(e)-1] do
 pp:= Position( enum_gens[pos[i]+pos[j]], e[k] );
 t:= Sum( enum_gens{[1..pos[i]+pos[j]-1]}, Length )+pp;
 Add( co, One( F )*e[k+1] );
 Add( co, t );
 od;
 SetEntrySCTable( T, i, j, co );
 fi;

 od;
 od;

 L:= LieAlgebraByStructureConstants( F, T );

 B:= Basis( L );

 # Now we compute the natural grading of `L'.

 grading:= [ ];
 k:= 1;

 for i in [1..Length(enum_gens)] do
 comp:= [ ];
 for j in [1..Length(enum_gens[i])] do
 Add( comp, B[k] );
 k:= k+1;
 od;
 Add( grading, Subspace( L, comp ) );
 od;

 Add( grading, Subspace( L, [ ] ) );

 SetGrading( L, rec( min_degree:= 1,
 max_degree:= Length( grading ) - 1,
 source:= Integers,
 hom_components:= function( d )
 if d in [1..Length(grading)] then
 return grading[d];
 else
 return Last(grading);
 fi;
 end
 )
 );

 vv:= BasisVectors( B );

 # Set the pth-power images of the basis elements of `B':## in <Cite Key="Maas2010"/>. Methods for <Ref Attr="EpimorphismSchurCover"/> are

 ## spin (projective) representation of the symmetric group of degree <A>n</A>
"zero" then
 pimgs[i]:= Zero( L )s# default. If no <A>p</A> is specified, 3 is used by default.

 e:= ExtRepOfObj##
 x:= Zero( L );
 for k in [1,3..Length(e)-1] do
 pp:= Position( enum_gens[pos[i]+1], e[k] );
 t:= Sum( enum_gens{[1..pos[i]]}, Length )+pp;
 x:= x+ One( F )*e[k+1]*vv[t];
 od;
 pimgs[i]:= x;
 fi;
 od;
 SetPthPowerImages( B, pimgs );
 SetIsRestrictedLieAlgebra( L, true );
 SetIsLieNilpotent( L, true );

 hm:= function( g, i )

 local h, e, x, k, pp, f, t;

 if not g in J[i] then
 Error("<g> is not an element of the i-th term of the series used to define <L>");
 fi;

 h:= Image( hom_pcg[i], Image(Homs[i], g ));
 e:= ExtRepOfObj(h);
 x:= Zero(L);
 for k in [1,3..Length(e)-1] do
 pp:= Position( enum_gens[i], e[k] );
 f:= GeneratorsOfGroup( pcgps[i] )[pp];
 t:= Position( gens, f );
 x:= x + e[k+1]*Basis(L)[t];
 od;
 return x;
 end ;

 SetNaturalHomomorphismOfLieAlgebraFromNilpotentGroup( L, hm );

 return L;

end );

Messung V0.5 in Prozent

¤ Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.0.313Bemerkung: ¤

*Bot Zugriff

Wurzel

Suchen

PVS Prover

Isabelle Prover

NIST Cobol Testsuite

Cephes Mathematical Library

Vienna Development Method

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.