Quellcode-Bibliothek lierep.gi Sprache: unbekannt

#############################################################################
##
##  This file is part of GAP, a system for computational discrete algebra.
##  This file's authors include Willem de Graaf, and Craig A. Struble.
##
##  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 modules over Lie algebras.
##

###########################################################################
##
#R  IsZeroCochainRep( <c> )
##
DeclareRepresentation( "IsZeroCochainRep", IsPackedElementDefaultRep, [1] );

##############################################################################
##
#M  Cochain( <V>, <s>, <list> )
##
##
InstallMethod( Cochain,
        "for a module over a Lie algebra, an integer and an object",
        true, [ IsAlgebraModule, IsInt, IsObject ], 0,
        function( V, s, obj )

    local fam,type;

         if IsLeftAlgebraModuleElementCollection( V ) then
             if IsRightAlgebraModuleElementCollection( V ) then
                Error("cochains are note defined for bi-modules");
             else
                 if not IsLieAlgebra( LeftActingAlgebra( V ) ) then
                     TryNextMethod();
                 fi;
             fi;
         else
             if not IsLieAlgebra( RightActingAlgebra( V ) ) then
                 TryNextMethod();
             fi;
         fi;

    # Every s-cochain has the same type, so we store the types in the
    # module. The family of an s-cochain knows about its order (s), and
    # about the underlying module. 0 is not a position in a list, so we store
    # the type of the 0-cochains elsewhere.

         if not IsBound( V!.cochainTypes ) then
            V!.cochainTypes:= [ ];
         fi;
         if s = 0 then
           if not IsBound( V!.zeroCochainType ) then
             fam:= NewFamily( "CochainFamily", IsCochain );
             fam!.order:= s;
             fam!.module:= V;
             type:= NewType( fam, IsZeroCochainRep );
             V!.zeroCochainType:= type;
           else
             type:= V!.zeroCochainType;
           fi;
           return Objectify( type, [ obj ] );
         fi;

         if not IsBound( V!.cochainTypes[s] ) then
            fam:= NewFamily( "CochainFamily", IsCochain );
            fam!.order:= s;
            fam!.module:= V;
            type:= NewType( fam, IsPackedElementDefaultRep );
            V!.cochainTypes[s]:= type;
         else
            type:= V!.cochainTypes[s];
         fi;
         return Objectify( type, [ Immutable( obj ) ] );

end );

##############################################################################
##
#M  ExtRepOfObj( <coch> ) . . . . . . . . . . . . . . . for a cochain
##
InstallMethod( ExtRepOfObj,
        "for a cochain",
        true, [ IsCochain and IsPackedElementDefaultRep ], 0,
        c -> c![1] );

##############################################################################
##
#M  PrintObj( <coch> ) . . . . . . . . . . . . . . . for cochains
##
##
InstallMethod( PrintObj,
       "for a cochain",
       true, [ IsCochain ], 0,
       function( c )

          Print("<",FamilyObj(c)!.order,"-cochain>");
end );

##############################################################################
##
#M  CochainSpace( <V>, <s> ) . . . . . . . for a module over a Lie algebra and
##                                         an integer
##
##
InstallMethod( CochainSpace,
     "for a module over a Lie algebra and an integer",
     true, [ IsAlgebraModule, IS_INT ], 0,
     function( V, s )

       local L,r,n,F,tups,bas,k,t,l;

       L:= ActingAlgebra( V );
       if not IsLieAlgebra( L ) then
         Error("<V> must be a module over a Lie algebra");
       fi;

       r:= Dimension( V );
       F:= LeftActingDomain( L );

       if r = 0 then
         if s = 0 then
           return VectorSpace( F, [], Cochain( V, 0, Zero(V) ), "basis" );
         else
           return VectorSpace( F, [], Cochain( V, s, [] ), "basis" );
         fi;
       fi;

       if s = 0 then
         bas:= List( BasisVectors( Basis( V ) ), x -> Cochain( V, s, x ) );
         return VectorSpace( F, bas, "basis" );
       fi;

       n:= Dimension( L );
       tups:= Combinations( [1..n], s );

    #Every tuple gives rise to `r' basis vectors.

       bas:= [ ];
       for k in [1..r] do
         for t in tups do
           l:= List( [1..r], x -> [] );
           Add( l[k], [ t, One( F ) ] );
           Add( bas, l );
         od;
       od;

       bas:= List( bas, x -> Cochain( V, s, x ) );
       FamilyObj( bas[1] )!.tuples:= tups;
       return VectorSpace( F, bas, "basis" );
end );

##############################################################################
##
#M  \+( <c1>, <c2> ) . . . . . . . . . . . . . . . . . . . for two cochains
#M  AdditiveInverseOp( <c> ) . . . . .  . . . . . . . . . . . . . . . . for a cochain
#M  \*( <scal>, <c> ) . . . . . . . . . . . . . . for a scalar and a cochain
#M  \*( <c>, <scal> ) . . . . . . . . . . . . . . for a chain and a scalar
#M  \<( <c1>, <c2> ) . . . . . . . . . . . . . . . . . . . for two cochains
#M  \=( <c1>, <c2> ) . . . . . . . . . . . . . . . . . . . for two cochains
#M  ZeroOp( <c> ) . . . . . . . . . . . . . . . . . . . .  for a cochain
##
InstallMethod( \+,
    "for two cochains",
    IsIdenticalObj, [ IsCochain and IsPackedElementDefaultRep,
            IsCochain and IsPackedElementDefaultRep ], 0,
    function( c1, c2 )

      local l1,l2,r,l,k,list,i;

      l1:= c1![1]; l2:= c2![1];
      r:= Length( l1 );

  # We `merge the two lists'.

      l:= [ ];
      for k in [1..r] do
        if l1[k] = [] then
          l[k]:= l2[k];
        elif l2[k] = [ ] then
          l[k]:= l1[k];
        else
          list:= List( l1[k], ShallowCopy );
          Append( list, List( l2[k], ShallowCopy ) );
          SortBy( list, t -> t[1] );
          i:= 1;
          while i < Length( list ) do  # take equal things together.
            if list[i][1] = list[i+1][1] then
               list[i][2]:= list[i][2]+list[i+1][2];
               Remove( list, i+1 );
            else
               i:= i+1;
            fi;
          od;
          list:= Filtered( list, x -> x[2]<>0*x[2] );
          l[k]:= list;
        fi;
      od;
      return Objectify( TypeObj( c1 ), [ Immutable( l ) ] );

end );

InstallMethod( \+,
    "for two 0-cochains",
    IsIdenticalObj, [ IsCochain and IsZeroCochainRep,
            IsCochain and IsZeroCochainRep ], 0,
    function( c1, c2 )

      return Objectify( TypeObj( c1 ), [ c1![1] + c2![1] ] );
end );

InstallMethod( AdditiveInverseOp,
     "for a cochain",
     true, [ IsCochain and IsPackedElementDefaultRep ], 0,
     function( c )

       local l,lc,k,i;

       l:= [ ];
       lc:= c![1];
       for k in [1..Length(lc)] do
         l[k]:= List( lc[k], ShallowCopy );
         for i in [1..Length(l[k])] do
           l[k][i][2]:= -l[k][i][2];
         od;
       od;
       return Objectify( TypeObj( c ), [ Immutable( l ) ] );
end );

InstallMethod( AdditiveInverseOp,
     "for a 0-cochain",
     true, [ IsCochain and IsZeroCochainRep ], 0,
     function( c )

      return Objectify( TypeObj( c ), [ -c![1] ] );
end );

InstallMethod( \*,
     "for scalar and cochain",
     true, [ IsScalar, IsCochain and IsPackedElementDefaultRep ], 0,
     function( scal, c )

       local l,lc,k,i;

       l:= [ ];
       lc:= c![1];
       for k in [1..Length(lc)] do
         l[k]:= List( lc[k], ShallowCopy );
         for i in [1..Length(l[k])] do
           l[k][i][2]:= scal*l[k][i][2];
         od;
       od;
       return Objectify( TypeObj( c ), [ Immutable( l ) ] );

end );

InstallMethod( \*,
     "for scalar and cochain",
     true, [ IsScalar and IsZero, IsCochain and IsPackedElementDefaultRep ], 0,
     function( scal, c )

       return Zero( c );
end );

InstallMethod( \*,
     "for scalar and 0-cochain",
     true, [ IsScalar, IsCochain and IsZeroCochainRep ], 0,
     function( scal, c )

       return Objectify( TypeObj( c ), [ scal*c![1] ] );
end );

InstallMethod( \*,
     "for cochain and scalar",
     true, [ IsCochain and IsPackedElementDefaultRep, IsScalar ], 0,
     function( c, scal )

       local l,lc,k,i;

       l:= [ ];
       lc:= c![1];
       for k in [1..Length(lc)] do
         l[k]:= List( lc[k], ShallowCopy );
         for i in [1..Length(l[k])] do
           l[k][i][2]:= scal*l[k][i][2];
         od;
       od;
       return Objectify( TypeObj( c ), [ Immutable( l ) ] );

end );

InstallMethod( \*,
     "for cochain and scalar",
     true, [ IsCochain and IsPackedElementDefaultRep, IsScalar and IsZero ], 0,
     function( c, scal )

       return Zero( c );
end );

InstallMethod( \*,
     "for 0-cochain and scalar",
     true, [ IsCochain and IsZeroCochainRep, IsScalar ], 0,
     function( c, scal )

        return Objectify( TypeObj( c ), [ scal*c![1] ] );
end );

InstallMethod( \<,
     "for two cochains",
     true, [ IsCochain and IsPackedElementDefaultRep,
             IsCochain and IsPackedElementDefaultRep ],0,
     function( c1, c2 )
        return c1![1]<c2![1];
end );

InstallMethod( \=,
     "for two cochains",
     true, [ IsCochain and IsPackedElementDefaultRep,
             IsCochain and IsPackedElementDefaultRep ],0,
     function( c1, c2 )
        return c1![1]=c2![1];
end );

InstallMethod( ZeroOp,
     "for a cochain",
     true, [ IsCochain and IsPackedElementDefaultRep ], 0,
     function( c )

        local list;

        list:= List( c![1], x -> [] );
        return Objectify( TypeObj( c ), [ Immutable( list ) ] );
end );

InstallMethod( ZeroOp,
     "for a 0-cochain",
     true, [ IsCochain and IsZeroCochainRep ], 0,
     function( c )

       return Objectify( TypeObj( c ), [ Zero( c![1] ) ] );
end );

#############################################################################
##
#M  NiceFreeLeftModuleInfo( <C> ) . . . . . . . . . for a module of cochains
#M  NiceVector ( <C>, <c> ) . . . . .for a module of cochains and a cochain
#M  UglyVector( <C>, <v> ) . . . . . for a module of cochains and a row vector
##
InstallHandlingByNiceBasis( "IsCochainsSpace", rec(
    detect := function( R, gens, V, zero )
      return IsCochainCollection( V );
      end,

    NiceFreeLeftModuleInfo := function( C )

        local G,tups,g,l,k,i;

  # We collect together the tuples occurring in the generators of `C'
  # and store them in `C'. If the dimension of `C' is small with respect
  # to the number of possible tuples, then this leads to smaller nice
  # vectors.

        if ElementsFamily( FamilyObj( C ) )!.order = 0 then
          return true;
        fi;

        G:= GeneratorsOfLeftModule( C );
        tups:= [ ];
        for g in G do
          l:= g![1];
          for k in [1..Length(l)] do
            for i in [1..Length(l[k])] do
              AddSet( tups, l[k][i][1] );
            od;
          od;
        od;
        return tups;
      end,

    NiceVector := function( C, c )
      local tt,l,v,k,i,p;

      if IsZeroCochainRep( c ) then
        return Coefficients( Basis( FamilyObj( c )!.module ), c![1] );
      elif not IsPackedElementDefaultRep( c ) then
        TryNextMethod();
      fi;
      tt:= NiceFreeLeftModuleInfo( C );
      l:= c![1];

   # Every tuple gives rise to dim V entries in the nice Vector
   # (where V is the Lie algebra module).

      v:= ListWithIdenticalEntries( Length(l)*Length(tt),
                                     Zero( LeftActingDomain( C ) ) );
      if v = [ ] then v:= [  Zero( LeftActingDomain( C ) ) ]; fi;

      for k in [1..Length(l)] do
        for i in [1..Length(l[k])] do
          p:= Position( tt, l[k][i][1] );
          if p = fail then return fail; fi;
          v[(k-1)*Length(tt)+p]:= l[k][i][2];
        od;
      od;
      return v;
      end,

    UglyVector := function( C, vec )
      local l,tt,k,j,i,fam;

  # We do the inverse of `NiceVector'.

      fam:= ElementsFamily( FamilyObj( C ) );
      if fam!.order = 0 then

        return Objectify( fam!.module!.zeroCochainType, [
                 LinearCombination( Basis( fam!.module ), vec ) ]  );
      fi;

      l:= [ ];
      tt:= NiceFreeLeftModuleInfo( C );
      k:= 1;
      j:=0;
      while j <> Length( vec ) do
        l[k]:= [ ];
        for i in [j+1..j+Length(tt)] do
          if vec[i] <> 0*vec[i] then
            Add( l[k], [ tt[i-j], vec[i] ] );
          fi;
        od;
        k:= k+1;
        j:= j+ Length(tt);
      od;

      return Objectify( fam!.module!.cochainTypes[ fam!.order ],
                       [ Immutable(l) ] );
      end ) );

##############################################################################
##
#F   ValueCochain( <c>, <y1>, ... ,<ys> )
##
##
InstallGlobalFunction( ValueCochain,
       function( arg )

         local c,ys,V,L,cfs,le,k,cfs1,i,j,cf,val,vs,ind,
               sign, # sign of a permutation.
               p,ec;

     # We also allow for lists as argument of the function.
     # Such a list must then consist of the listed arguments.

         if IsList( arg[1] ) then arg:= arg[1]; fi;

         c:= arg[1];
         if not IsCochain( c ) then
           Error( "first arggument must be a cochain" );
         fi;

         if FamilyObj( c )!.order = 0 then
           return c![1];
         fi;

         ys:= arg{[2..Length(arg)]};
         if Length( ys ) <> FamilyObj( c )!.order then
           Error( "number of arguments is not equal to the order of <c>" );
         fi;

         V:= FamilyObj( c )!.module;
         L:= ActingAlgebra( V );
         cfs:= [ List( ys, x -> Coefficients( Basis(L), x ) ) ];
         le:= Length( ys );
         k:= 1;

   # We expand the list of coefficients to a list of elements of the form
   #
   #         [ [2/3,2], [7,1], [1/3,3] ]
   #
   # meaning that there we have to evaluate 2/3*7*1/3*c( x_2, x_1, x_3 ).

         while k <= le do

           cfs1:= [ ];
           for i in [1..Length(cfs)] do
             for j in [1..Length(cfs[i][k])]  do
               if cfs[i][k][j] <> 0*cfs[i][k][j] then
                  cf:= ShallowCopy( cfs[i] );
                  cf[k] := [ cf[k][j], j ];
                  Add( cfs1, cf );
               fi;
             od;
           od;
           cfs:= cfs1;
           k:= k+1;
         od;

   # We loop over the expanded list, and add the values that we get.

         ec:= c![1];
         val:= Zero( V );
         vs:= BasisVectors( Basis( V ) );
         for i in [1..Length( cfs )] do
           cf:= Product( List( cfs[i], x -> x[1] ) );
           ind:= List( cfs[i], x -> x[2] );
           sign:= SignPerm( Sortex( ind ) );
           for k in [1..Length(ec)] do
             p:= PositionProperty( ec[k], x -> x[1] = ind );
             if p <> fail then
               val:= val +  ec[k][p][2]*sign*cf*vs[k];
             fi;
           od;
         od;

         return val;

end );

#############################################################################
##
#V  LieCoboundaryOperator
##
##  Takes an s-cochain, and returns an (s+1)-cochain.
##
InstallGlobalFunction( LieCoboundaryOperator,

     function( c )

       local s,V,L,bL,n,fam,tups,type,list,t,val,cfs,k,q,r,elts,z,inp,sn,F;

       s:= FamilyObj( c )!.order;
       V:= FamilyObj( c )!.module;
       L:= ActingAlgebra( V );
       bL := BasisVectors( Basis( L ) );
       n:= Dimension( L );
       F:= LeftActingDomain( V );

   # We get the type of the (s+1)-cochains, and store the tuples we need
   # in the family (so that in the next call of `LieCoboundaryOperator' we
   # do not need to recompute them).

       if IsBound( V!.cochainTypes[s+1] ) then
          fam:= FamilyType( V!.cochainTypes[s+1] );
          if IsBound( fam!.tuples ) then
            tups:= fam!.tuples;
          else
            tups:= Combinations( [1..n], s+1 );
            fam!.tuples:= tups;
          fi;
       else
          tups:= Combinations( [1..n], s+1 );
          fam:= NewFamily( "CochainFamily", IsCochain );
          fam!.order:= s+1;
          fam!.module:= V;
          fam!.tuples:= tups;
          type:= NewType( fam, IsPackedElementDefaultRep );
          V!.cochainTypes[s+1]:= type;
       fi;

       list:= List( [1..Dimension(V)], x -> [] );
       for t in tups do

   # We calculate \delta(c)(x_{i_1},...,x_{i_s+1}) (where \delta denotes
   # the coboundary operator). We use the definition of \delta as given in
   # Jacobson, Lie Algebras, Dover 1979, p. 94. There he writes about right
   # modules. We cater for left and right modules; for left modules we have
   # to add a - when acting.

         val:= Zero( V );
         sn:= (-1)^s;
         for q in [1..s+1] do
           elts:= bL{t};
           z:= elts[q];
           Remove( elts, q );
           inp:= [c]; Append( inp, elts );
           if IsLeftAlgebraModuleElementCollection( V ) then
             val:= val - sn*( z^ValueCochain( inp ) );
           else
             val:= val + sn*( ValueCochain( inp )^z );
           fi;
           sn:= -sn;

           for r in [q+1..s+1] do
             elts:= bL{t};
             z:= elts[q]*elts[r];
             Unbind( elts[q] ); Unbind( elts[r] );
             elts:= Compacted( elts );
             inp:= [ c ]; Append( inp, elts ); Add( inp, z );
             val:= val+(-1)^(q+r)*ValueCochain( inp );
           od;
         od;

         cfs:= Coefficients( Basis(V), val );
         for k in [1..Length(cfs)] do
           if cfs[k] <> 0*cfs[k] then
             Add( list[k], [ t, cfs[k] ] );
           fi;
         od;

       od;

       return Cochain( V, s+1, list );

end );

##############################################################################
##
#M  Coboundaries( <V>, <s> ) . . . . . . . . . for alg module and integer
##
##
InstallMethod( Coboundaries,
    "for module over a Lie algebra and an integer",
    true, [ IsAlgebraModule, IS_INT ], 0,
    function( V, s )

      local Csm1,gens;

   # if s=0, then the space is zero.

      if s = 0 then
          return VectorSpace( LeftActingDomain(V),
                         [ ], Cochain( V, 0, Zero(V) ), "basis" );
      fi;

   # The s-coboundaries are the images of the (s-1)-cochains under
   # the coboundary operator.

      Csm1:= CochainSpace( V, s-1 );
      gens:= List( GeneratorsOfLeftModule( Csm1 ), x ->
                                       LieCoboundaryOperator(x) );
      if Length(gens) = 0 then
          return VectorSpace( LeftActingDomain(V),
                         [ ], Cochain( V, s, [] ), "basis" );
      fi;
      return VectorSpace( LeftActingDomain(V), gens );

end );

InstallMethod( Cocycles,
    "for module over a Lie algebra and an integer",
    true, [ IsAlgebraModule, IS_INT ], 0,
    function( V, s )

      local Cs,gens,Bsp1,B,eqmat,sol;

  # The set of s-cocycles is the kernel of the coboundary operator,
  # when restricted to the space of s-cochains.

      Cs:= CochainSpace( V, s );
      if IsTrivial(Cs) then return Cs; fi;
      gens:= List( GeneratorsOfLeftModule( Cs ), x ->
                                       LieCoboundaryOperator(x) );

      Bsp1:= VectorSpace( LeftActingDomain(V), gens );
      B:= Basis( Bsp1 );

      if Dimension( Bsp1 ) > 0 then
          eqmat:= List( gens, x -> Coefficients( B, x ) );
          sol:= NullspaceMat( eqmat );
          sol:= List( sol, x -> LinearCombination(
                        GeneratorsOfLeftModule(Cs),x));
          return Subspace( Cs, sol, "basis" );
      else
          # in this case the every cochain is a cocycle.
          return Cs;
      fi;

end );

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

          local   C,  refl,  rank,  i,  m,  G,  RM,  j;

          C:= CartanMatrix( R );

      # We calculate a list of simple reflections that generate the
      # Weyl group. The reflections are given by the matrices of their
      # action on the fundamental weights. Let r_i denote the i-th
      # simple reflection, and \lambda_i the i-th fundamental weight.
      # Then r_i(\lambda_j) = \lambda_j -\delta_{ij} \alpha_i, where
      # \alpha_i is the i-th simple root. Furthermore, in the basis of
      # fundamental weights the coefficients of the simple root \alpha_i
      # are the i-th row of the Cartan matrix C. So the matrix of the
      # i-th reflection is the identity matrix, with C[i] subtracted
      # from the i-th row. So the action of a reflection with matrix m
      # on a weight \mu = [ n_1,.., n_l] (list of integers) is given by
      # \mu*m.

          refl:= [ ];
          rank:= Length( C);
          for i in [1..rank] do
              m:= IdentityMat( rank, rank );
              m[i]:=m[i]-C[i];
              Add( refl, m );
          od;
          G:= Group( refl );
          SetIsWeylGroup( G, true );
          RM:=[];
          for i in [1..rank] do
              RM[i]:= [ ];
              for j in [1..rank ] do
                  if C[i][j] <> 0 then
                      Add( RM[i], [j,C[i][j]] );
                  fi;
              od;
          od;
          SetSparseCartanMatrix( G, RM );
          SetRootSystem( G, R );
          return G;
end );

#############################################################################
##
#M  ApplySimpleReflection( <SC>, <i>, <w> )
##
##
InstallMethod( ApplySimpleReflection,
   "for a sparse Cartan matrix, index and weight",
   true, [ IsList, IS_INT, IsList ], 0,

function( SC, i, w )

          local   p, ni;

          ni:= w[i];
          if ni = 0 then return; fi;
          for p in SC[i] do
              w[p[1]]:= w[p[1]]-ni*p[2];
          od;

end );

############################################################################
##
#M  LongestWeylWordPerm( <W> ) . . . . . . . . . . . . . . . for a Weyl group
##
##
InstallMethod(LongestWeylWordPerm,
              "for Weyl group",
              true, [ IsWeylGroup ], 0,
          function( W )

          local   M,  rho,  p;

          M:= SparseCartanMatrix( W );

     # rho will be the Weyl vector (in the basis of fundamental weights).

          rho:= List( [1..Length(M)], x -> -x );
          p:= 1;

          while p <> fail do
              ApplySimpleReflection( M, p, rho );
              p:= PositionProperty( rho, x -> x < 0 );
          od;

          return PermList( List( [1..Length(M)], x -> Position( rho, x ) ) );

end );

#############################################################################
##
#M  ConjugateDominantWeight( <W>, <w> )
##
##
InstallMethod( ConjugateDominantWeight,
              "for Weyl group and weight",
              true, [ IsWeylGroup, IsList ], 0,
              function( W, wt )

          local   ww,  M,  p;

          ww:= ShallowCopy( wt );
          M:= SparseCartanMatrix( W );
          p:= PositionProperty( ww, x -> x < 0 );

      # We apply simple reflections until `ww' is dominant.

          while p <> fail do
              ApplySimpleReflection( M, p, ww );
              p:= PositionProperty( ww, x -> x < 0 );
          od;
          return ww;

end);

###########################################################################
##
#M  ConjugateDominantWeightWithWord( <W>, <wt> )
##
##
InstallMethod( ConjugateDominantWeightWithWord,
              "for Weyl group and weight",
              true, [ IsWeylGroup, IsList ], 0,
              function( W, wt )

          local   ww,  M,  p, word;

          ww:= ShallowCopy( wt );
          word:= [ ];
          M:= SparseCartanMatrix( W );
          p:= PositionProperty( ww, x -> x < 0 );
          while p <> fail do
              ApplySimpleReflection( M, p, ww );
              Add( word, p );
              p:= PositionProperty( ww, x -> x < 0 );
          od;
          return [ ww, word ];
end);

#############################################################################
##
#M  WeylOrbitIterator( <w>, <wt> )
##
##  stack is a stack of weights, i.e., a list of elts of the form [ w, ind ]
##  the last elt of this list [w1,i1] is such that the i1-th refl app to
##  w1 gives currentweight. The second to last elt [w2,i2] is such that the
##  i2-th refl app to w2 gives w1 etc.
##
##  the status indicates whether or not to compute a successor
##     status=1 means output current weight w, next one will be g_0(w)
##     status=2 means output g_0(w), where w=current weight, next one will
##              be the successor of w
##     status=3 means output current weight w, next one will be succ(w)
##
##  midLen is the middle length, to where we have to compute
##  permuteMidLen is true if we have to map the weights of length
##  midLen with the longest Weyl element...
##

############################################################################
##
#M  IsDoneIterator( <it> ) . . . . . . . . . . . . for Weyl orbit iterator
##
BindGlobal( "IsDoneIterator_WeylOrbit", it -> it!.isDone );

############################################################################
##
#M  NextIterator( <it> ) . . . . . . . . . . . . for a Weyl orbit iterator
##
##  The algorithm is due to D. M. Snow (`Weyl group orbits',
##  ACM Trans. Math. Software, 16, 1990, 94--108).
##
BindGlobal( "NextIterator_WeylOrbit", function( it )
    local   output,  mu,  rank,  len,  stack,  bound,  foundsucc,
            pos,  i,  nu,  a;

    if it!.isDone then Error("the iterator is exhausted"); fi;

    if it!.status = 1 then
        it!.status:= 2;
        mu:= it!.currentWeight;
        if mu = 0*mu then
            it!.isDone:= true;
        fi;
        return mu;
    fi;

    if it!.status = 2 then
        output:= -Permuted( it!.currentWeight, it!.perm );
    else
        output:= ShallowCopy( it!.currentWeight );
    fi;

    #calculate the successor of curweight

    mu:= ShallowCopy(it!.currentWeight);
    rank:= Length( mu );
    len:= it!.curLen;
    stack:= it!.stack;
    bound:= 1;
    foundsucc:= false;
    while not foundsucc do

        pos:= fail;
        if len <> it!.midLen then
            for i in [bound..rank] do
                if mu[i]>0 then
                    nu:= ShallowCopy(mu);
                    ApplySimpleReflection( it!.RMat, i, nu );
                    if ForAll( nu{[i+1..rank]}, x -> x >= 0 ) then
                        pos:= i; break;
                    fi;
                fi;
            od;
        fi;

        if pos <> fail then
            Add( stack, [ mu, pos ] );
            foundsucc:= true;
        else

            if mu = it!.root then

                # we cannot find a successor of the root: we are done

                it!.isDone:= true;
                nu:= [];
                foundsucc:= true;
            else
                a:= Remove( stack );
                mu:= a[1]; bound:= a[2]+1;
                len:= len-1;
            fi;

        fi;

    od;

    it!.stack:= stack;
    it!.curLen:= len+1;
    it!.currentWeight:= nu;
    if len+1 = it!.midLen and not it!.permuteMidLen then
        it!.status:= 3;
    else
        it!.status:= 1;
    fi;

    return output;

end );

InstallMethod( WeylOrbitIterator,
        "for weights of a W-orbit",
        [ IsWeylGroup, IsList ],

        function( W, wt )

    local   mu,  perm,  nu,  len,  i;

    # The iterator starts at the dominant weight of the orbit.

    mu:= ConjugateDominantWeight( W, wt );

    # We calculate the maximum length occurring in an orbit (the length of
    # an element of the orbit being defined as the minimum number of
    # simple reflections that have to be applied in order to get from the
    # dominant weight to the particular orbit element). This will determine
    # whether we also have to apply the longest Weyl element to the elements
    # of "middle" length.

    perm:= LongestWeylWordPerm(W);
    nu:= -Permuted( mu, perm );
    len:= 0;
    while nu <> mu do
        i:= PositionProperty( nu, x -> x < 0 );
        ApplySimpleReflection( SparseCartanMatrix(W), i, nu );
        len:= len+1;
    od;

    return IteratorByFunctions( rec(
               IsDoneIterator := IsDoneIterator_WeylOrbit,
               NextIterator   := NextIterator_WeylOrbit,
#T no `ShallowCopy'!
               ShallowCopy:= function( iter )
                      return rec( root:= ShallowCopy( iter!.root ),
                        currentWeight:= ShallowCopy( iter!.currentWeight ),
                        stack:= ShallowCopy( iter!.stack ),
                        RMat:= iter!.RMat,
                        perm:= iter!.perm,
                        status:= iter!.status,
                        permuteMidLen:=  iter!.permuteMidLen,
                        midLen:=  iter!.midLen,
                        curLen:= iter!.curLen,
                        maxlen:= iter!.maxlen,
                        noPosR:= iter!.noPosR,
                        isDone:= iter!.isDone );
                     end,
                        root:= mu,
                        currentWeight:= mu,
                        stack:= [ ],
                        RMat:= SparseCartanMatrix(W),
                        perm:= perm,
                        status:= 1,
                        permuteMidLen:=  IsOddInt( len ),
                        midLen:=  EuclideanQuotient( len, 2 ),
                        curLen:= 0,
                        maxlen:= len,
                        noPosR:= Length( PositiveRoots(
                                RootSystem(W) ) ),
                        isDone:= false ) );
end );

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

      local posR,V,lcombs;

      posR:= PositiveRoots( R );
      V:= VectorSpace( Rationals, SimpleSystem( R ) );
      lcombs:= List( posR, r ->
                       Coefficients( Basis( V, SimpleSystem(R) ), r ) );
      return List( lcombs, c -> LinearCombination( CartanMatrix(R), c ) );

end );

#############################################################################
##
#M  DominantWeights( <R>, <maxw> )
##
InstallMethod( DominantWeights,
    "for a root system and a dominant weight",
    true, [ IsRootSystem, IsList ], 0,
    function( R, maxw )

    local n,posR,V,lcombs,dom,ww,newdom,mu,a,levels,heights,pos;

   # First we calculate the list of positive roots, represented in the
   # basis of fundamental weights. `heights' will be the list of heights
   # of the positive roots.

   posR:= PositiveRoots( R );
   V:= VectorSpace( Rationals, SimpleSystem( R ) );
   lcombs:= List( posR, r -> Coefficients( Basis( V, SimpleSystem(R) ), r ) );
   posR:= List( lcombs, c -> LinearCombination( CartanMatrix(R), c ) );

   heights:= List( lcombs, Sum );

   # Now `dom' will be the list of dominant weights; `levels' will be a list
   # (in bijection with `dom') of the levels of the weights in `dom'.

   dom:= [ maxw ];
   levels:= [ 0 ];

   ww:= [ maxw ];

   # `ww' is the list of weights found in the last round. We subtract the
   # positive roots from the elements of `ww'; algorithm as in
   # R. V. Moody and J. Patera, "Fast recursion formula for weight
   # multiplicities", Bull. Amer. math. Soc., 7:237--242.

   while ww <> [] do

     newdom:= [ ];
     for mu in ww do
       for a in posR do
         if ForAll( mu-a, x -> x >= 0 ) and not (mu-a in dom) then
           Add( newdom, mu - a );
           Add( dom, mu-a );
           pos:= Position( mu, dom );
           Add( levels, levels[Position(dom,mu)]+heights[Position(posR,a)] );
         fi;
       od;
     od;
     ww:= newdom;

   od;

   return [dom,levels];

end );

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

     local C, B, roots, i, j;

     C:= CartanMatrix( R );
     B:= NullMat( Length(C), Length(C) );
     roots:= ShallowCopy( PositiveRoots( R ) );
     Append( roots, NegativeRoots( R ) );

     # First we calculate the lengths of the roots. For that we use
     # the following. We have that $\kappa( h_i, h_i ) = \sum_{r\in R}
     # r(h_i)^2$, where $\kappa$ is the Killing form, and the $h_i$
     # are the canonical Cartan generators. Furthermore,
     # $(\alpha_i, \alpha_i) = 4/\kappa(h_i,h_i)$. We note that the roots
     # of R are represented on the basis of the $h_i$, so the $i$-th
     # element of a root $r$, is the value $r(h_i)$.

     for i in [1..Length(C)] do
       B[i][i]:= 4/Sum( List( roots, r -> r[i]^2 ) );
     od;

     # Now we calculate the other entries of the matrix of the bilinear
     # form.

     for i in [1..Length(C)] do
       for j in [i+1..Length(C)] do
         if C[i][j] <> 0 then
           B[i][j]:= C[i][j]*B[j][j]/2;
           B[j][i]:= B[i][j];
         fi;
       od;
     od;

     return B;

end );

#############################################################################
##
#M  DominantCharacter( <R>, <maxw> )
#M  DominantCharacter( <L>, <maxw> )
##
InstallMethod( DominantCharacter,
    "for a root system and a highest weight",
    true, [ IsRootSystem, IsList ], 0,
   function( R, maxw )

   local ww, rank, fundweights, rhs, bilin, i, j, rts, dones, mults,
         lam_rho, clam, WR, refl, grps, orbs, k, mu, zeros, p, O, W, reps,
         sum, a, done_summing, sum1, nu, nu1, mu_rho, gens;

   ww:= DominantWeights( R, maxw );
   rank:= Length( CartanMatrix( R ) );

   # `fundweights' will be a list of the fundamental weights, calculated
   # on the basis of simple roots. `bilin' will be the matrix of the
   # bilinear form of `R', relative to the fundamental weights.
   # We have that $(\lambda_i,\lambda_j) = \zeta_{ji} (\alpha_i,\alpha_i)/2$,
   # where $\zeta_{ji}$ is the $i$-th coefficient in the expression for
   # $\lambda_j$ as a linear combination of simple roots.

   fundweights:= [ ];
   for i in [1..rank] do
     rhs:= ListWithIdenticalEntries( rank, 0 );
     rhs[i]:= 1;
     Add( fundweights, SolutionMat( CartanMatrix(R), rhs ) );
   od;

   bilin:= NullMat( rank, rank );
   for i in [1..rank] do
     for j in [i..rank] do
       bilin[i][j]:= fundweights[j][i]*BilinearFormMat( R )[i][i]/2;
       bilin[j][i]:= bilin[i][j];
     od;
   od;

   # We sort the dominant weights according to level.

   SortParallel( ww[2], ww[1] );

   rts:= ShallowCopy( PositiveRootsAsWeights( R ) );
   Append( rts, -rts );

   # `dones' will be a list of the dominant weights for which we have
   # calculated the multiplicity. `mults' will be a list containing the
   # corresponding multiplicities. `lam_rho' is the weight `maxw+rho',
   # where `rho' is the Weyl vector.

   dones:= [ maxw ];
   mults:= [ 1 ];

   lam_rho:= maxw+List([1..rank], x -> 1 );
   clam:= lam_rho*(bilin*lam_rho);

   WR:= WeylGroup( R );
   refl:= GeneratorsOfGroup( WR );

   # `grps' is a list containing the index lists for the stabilizers of the
   # different weights (i.e., such a stabilizer is generated by the
   # simple reflections corresponding to the indices). `orbs' is a list
   # of orbits of these groups (acting on the roots).

   grps:= [ ]; orbs:= [ ];

   for k in [2..Length(ww[1])] do

       mu:= ww[1][k];

       # We calculate the multiplicity of `mu'. The algorithm is as
       # described in
       # R. V. Moody and J. Patera, "Fast recursion formula for weight
       # multiplicities", Bull. Amer. math. Soc., 7:237--242.
       # First we calculate the orbits of the stabilizer of `mu' (with the
       # additional element -1), acting on the roots.

       zeros:= Filtered([1..rank], x -> mu[x]=0 );
       p:= Position( grps, zeros );
       if p <> fail then
           O:= orbs[p];
       else

           gens:= refl{zeros};
           Add( gens, -IdentityMat(rank) );
           W:= Group( gens );
           O:= Orbits( W, rts );
           Add( grps, zeros );
           Add( orbs, O );
       fi;

       # For each representative of the orbits we calculate the sum occurring
       # in Freudenthal's formula (and multiply by the size of the orbit).

       reps:= List( O, o -> Intersection( o, PositiveRootsAsWeights(R) )[1] );
       sum:= 0;
       for i in [1..Length(reps)] do
           a:= reps[i];
           j:= 1; done_summing:= false;
           sum1:= 0;
           while not done_summing do
               nu:= mu+j*a;
               nu1:= ConjugateDominantWeight( WR, nu );
               if not nu1 in ww[1] then
                   done_summing:= true;
               else

                   p:= Position( dones, nu1 );
                   sum1:= sum1 + mults[p]*(nu*(bilin*a));
                   j:= j+1;
               fi;
           od;
           sum:= sum + Length(O[i])*sum1;
       od;

       mu_rho:= mu+List([1..rank],x->1);

       sum:= sum/( clam - mu_rho*(bilin*mu_rho) );
       Add( dones, mu );
       Add( mults, sum );

   od;

   return [ dones, mults ];

end );

InstallOtherMethod( DominantCharacter,
    "for a semisimple Lie algebra and a highest weight",
    true, [ IsLieAlgebra, IsList ], 0,
   function( L, maxw )
       return DominantCharacter( RootSystem(L), maxw );
end );

###############################################################################
##
#M  DecomposeTensorProduct( <L>, <w1>, <w2> )
##
##
InstallMethod( DecomposeTensorProduct,
     "for a semisimple Lie algebra and two dominant weights",
     true, [ IsLieAlgebra, IsList, IsList ], 0,
    function( L, w1, w2 )

    #W decompose the tensor product of the two irreps of L with hwts
    #w1,w2 respectively. We use Klymik's formula.

    local   R,  W,  ch1,  wts,  mlts,  rho,  i,  it,  ww,  mu,  nu,
            mult,  p;

    R:= RootSystem( L );
    W:= WeylGroup( R );
    ch1:= DominantCharacter( L, w1 );
    wts:= [ ]; mlts:= [ ];
    rho:= ListWithIdenticalEntries( Length( CartanMatrix( R ) ), 1 );

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

       # We loop through all weights of the irrep with highest weight <w1>.
       # We get these by taking the orbits of the dominant ones under the
       # Weyl group.

        it:= WeylOrbitIterator( W, ch1[1][i] );
        while not IsDoneIterator( it ) do

            ww:= NextIterator( it ); #+w2+rho;
            ww:= ww+w2+rho;
            mu:= ConjugateDominantWeightWithWord( W, ww );

            if not ( 0 in mu[1] ) then

              # The stabilizer of `ww' is trivial; so `ww' contributes to the
              # formula. `nu' will be the highest weight of the direct
              # summand gotten from `ww'.

                nu:= mu[1]-rho;
                mult:= ch1[2][i]*( (-1)^Length(mu[2]) );
                p:= PositionSorted( wts, nu );
                if not IsBound( wts[p] ) or wts[p] <> nu then
                    Add( wts, nu, p );
                    Add( mlts, mult, p );
                else
                    mlts[p]:= mlts[p]+mult;
                    if mlts[p] = 0 then
                        Remove( mlts, p );
                        Remove( wts, p );
                    fi;

                fi;
            fi;
        od;
    od;
    return [ wts, mlts ];

end );

###############################################################################
##
#M  DimensionOfHighestWeightModule( <L>, <w> )
##
##
InstallMethod( DimensionOfHighestWeightModule,
        "for a semisimple Lie algebra",
        true, [ IsLieAlgebra, IsList ], 0,
        function( L, w )

    local   R,  l,  B,  M,  p,  r,  cf,  den,  num,  i;

    R:= RootSystem( L );
    l:= Length( CartanMatrix( R ) );
    B:= Basis( VectorSpace( Rationals, SimpleSystem(R) ), SimpleSystem(R) );
    M:= BilinearFormMat( R );
    p:= 1;
    for r in PositiveRoots( R ) do
        cf:= Coefficients( B, r );
        den:= 0;
        num:= 0;
        for i in [1..l] do
            num:= num + cf[i]*(w[i]+1)*M[i][i];
            den:= den + cf[i]*M[i][i];
        od;
        p:= p*(num/den);
    od;

    return p;

end );

############################################################################
##
#M  ObjByExtRep( <fam>, <list> )
#M  ExtRepOfObj( <obj> )
##
InstallMethod( ObjByExtRep,
   "for family of UEALattice elements, and list",
   true, [ IsUEALatticeElementFamily, IsList ], 0,
   function( fam, list )
#+
    return Objectify( fam!.packedUEALatticeElementDefaultType,
                    [ Immutable(list) ] );
end );

InstallMethod( ExtRepOfObj,
   "for an UEALattice element",
   true, [ IsUEALatticeElement ], 0,
   function( obj )
#+
   return obj![1];

end );

###########################################################################
##
#M  PrintObj( <m> ) . . . . . . . . . . . . . . . . for an UEALattice element
##
InstallMethod( PrintObj,
        "for UEALattice element",
        true, [IsUEALatticeElement and IsPackedElementDefaultRep], 0,
        function( x )

    local   lst,  k, i, n;

    # This function prints a UEALattice element; see notes above.

    lst:= x![1];
    n:= FamilyObj( x )!.noPosRoots;
    if lst=[] then
        Print("0");
    else
        for k in [1,3..Length(lst)-1] do
            if lst[k+1] > 0 and k>1 then
                Print("+" );
            fi;
            if not IsOne(lst[k+1]) then
                Print( lst[k+1],"*");
            fi;
            if lst[k] = [] then
                Print("1");
            else

                for i in [1,3..Length(lst[k])-1] do
                    if lst[k][i] <=n then
                        Print("y",lst[k][i]);
                        if lst[k][i+1]>1 then
                            Print("^(",lst[k][i+1],")");
                        fi;
                    elif lst[k][i] <= 2*n then
                        Print("x",lst[k][i]-n);
                        if lst[k][i+1]>1 then
                            Print("^(",lst[k][i+1],")");
                        fi;
                    else
                        Print("( h",lst[k][i],"/",lst[k][i+1]," )");
                    fi;
                    if i <> Length(lst[k])-1 then
                        Print("*");
                    fi;
                od;
            fi;

        od;

    fi;

end );

#############################################################################
##
#M  OneOp( <m> ) . . . . . . . . . . . . . . . . for a UEALattice element
#M  ZeroOp( <m> ) . . . . . . . . . . . . . . .  for a UEALattice element
#M  \<( <m1>, <m2> ) . . . . . . . . . . . . . . for two UEALattice elements
#M  \=( <m1>, <m2> ) . . . . . . . . . . . . . . for two UEALattice elements
#M  \+( <m1>, <m2> ) . . . . . . . . . . . . . . for two UEALattice elements
#M  \AdditiveInverseOp( <m> )     . . . . . . . . . . . . . . for a UEALattice element
##
##
InstallMethod( OneOp,
        "for UEALattice element",
        true, [ IsUEALatticeElement and IsPackedElementDefaultRep ], 0,
        function( x )

    return ObjByExtRep( FamilyObj( x ), [ [], 1 ] );

end );

InstallMethod( ZeroOp,
        "for UEALattice element",
        true, [ IsUEALatticeElement and IsPackedElementDefaultRep ], 0,
        function( x )

    return ObjByExtRep( FamilyObj( x ), [ ] );

end );

InstallMethod( \<,
                "for two UEALattice elements",
        IsIdenticalObj, [ IsUEALatticeElement and IsPackedElementDefaultRep,
                IsUEALatticeElement and IsPackedElementDefaultRep ], 0,
        function( x, y )
    return x![1]< y![1];
end );

InstallMethod( \=,
                "for two UEALattice elements",
        IsIdenticalObj, [ IsUEALatticeElement and IsPackedElementDefaultRep,
                IsUEALatticeElement and IsPackedElementDefaultRep ], 0,
        function( x, y )

    return x![1] = y![1];
end );

InstallMethod( \+,
        "for two UEALattice elements",
        true, [ IsUEALatticeElement and IsPackedElementDefaultRep,
                IsUEALatticeElement and IsPackedElementDefaultRep], 0,
        function( x, y )

    return ObjByExtRep( FamilyObj(x), ZippedSum( x![1], y![1], 0, [\<,\+] ) );
end );

InstallMethod( AdditiveInverseOp,
        "for UEALattice element",
        true, [ IsUEALatticeElement and IsPackedElementDefaultRep ], 0,
        function( x )

    local   ex,  i;

    ex:= ShallowCopy(x![1]);
    for i in [2,4..Length(ex)] do
        ex[i]:= -ex[i];
    od;
    return ObjByExtRep( FamilyObj(x), ex );
end );

#############################################################################
##
#M  \*( <scal>, <m> ) . . . . . . . . .for a scalar and a UEALattice element
#M  \*( <m>, <scal> ) . . . . . . . . .for a scalar and a UEALattice element
##
InstallMethod( \*,
        "for scalar and UEALattice element",
        true, [ IsScalar, IsUEALatticeElement and
                IsPackedElementDefaultRep ], 0,
        function( scal, x )

    local   ex,  i;

    ex:= ShallowCopy( x![1] );
    for i in [2,4..Length(ex)] do
        ex[i]:= scal*ex[i];
    od;
    return ObjByExtRep( FamilyObj(x), ex );
end);

InstallMethod( \*,
        "for UEALattice element and scalar",
        true, [ IsUEALatticeElement and IsPackedElementDefaultRep,
                IsScalar ], 0,
        function( x, scal )

    local   ex,  i;

    ex:= ShallowCopy( x![1] );
    for i in [2,4..Length(ex)] do
        ex[i]:= scal*ex[i];
    od;
    return ObjByExtRep( FamilyObj(x), ex );
end);

#############################################################################
##
#F  CollectUEALatticeElement( <noPosR>, <BH>, <f>, <vars>, <Rvecs>, <RT>,
##                                                          <posR>, <lst> )
##
InstallGlobalFunction( CollectUEALatticeElement,

    function( noPosR, BH, f, vars, Rvecs, RT, posR, lst )

    local   i, j, k, l, p, q, r, s,   # loop variables
            todo,                # list of monomials that still need treatment
            dones,               # list of monomials that don't
            collocc,             # `true' is a collection has occurred
            mon, mon1,           # monomials,
            cf, c1, c2, c3, c4,  # coefficients
            temp,                # for temporary storing
            start, tail,         # beginning and end of a monomial
            h,                   # Cartan element
            rr,                  # list of monomials with coefficients
            type,                # type of a pseudo root system of rank 2
            i1, j1, m, n,        # integers
            a, b,                # roots
            p1, p2, p3, p4,      # positions
            st1,
            has_h,
            mons,
            pol,
            ww,
            mm,
            min,
            WriteAsLCOfBinoms;   # local function.

     WriteAsLCOfBinoms:= function( vars, pol )

        # This function writes the polynomial `pol' in the variables `vars'
        # as a linear combination of polynomials of the form
        # (x_1\choose m_1).....(x_t\choose m_t). (`pol' must tae integral
        # values when evaluated at integral points.)

         local   d,  ind,  e,  fam,  fac,  k,  p,  q,  bin,  cc,  res,
                 mon, dfac;

         if IsConstantRationalFunction( pol ) or vars = [] then
             return [ [], pol ];
         fi;
         d:=  DegreeIndeterminate(pol,vars[1]);
         if d = 0 then
            # The variable `vars[1]' does not occur in `pol', so we can
            # recurse with one variable less.
             return WriteAsLCOfBinoms( vars{[2..Length(vars)]}, pol );
         fi;

         ind:= IndeterminateNumberOfLaurentPolynomial( vars[1] );
         e:= ShallowCopy( ExtRepPolynomialRatFun( pol ) );
         fam:= FamilyObj( pol );

         # `fac' will be contain the monomials of degree `d' in the variable
         # `vars[1]'.
         fac:= [ ];
         for k in [1,3..Length(e)-1] do

             if e[k]<>[] and e[k][1] = ind and e[k][2] = d then
                 Add( fac, e[k] ); Unbind( e[k] );
                 Add( fac, e[k+1] ); Unbind( e[k+1] );
             fi;
         od;
         e:= Compacted( e );
         # `e' now contains the rest of the polynomial.

         p:= PolynomialByExtRepNC( fam, fac )/(vars[1]^d);
         q:= PolynomialByExtRepNC( fam, e );

         # So now we have `pol = vars[1]^d*p+q', where `p' does not contain
         # `vars[1]' and `q' has lower degree in `vars[1]'. We can also
         # write this as (writing x = vars[1])
         #
         #            (x)            (x)
         #    pol = d!(d)p + q - { d!(d) - x^d }p
         #
         # `bin' will be d!* x\choose d.

         bin:= Product( List( [0..d-1], x -> vars[1] - x ) );
         q:= q - (bin-vars[1]^d)*p;
         cc:= WriteAsLCOfBinoms( vars{[2..Length(vars)]}, p );

         # No wwe prepend d!*(x\choose d) to cc.
         dfac := Factorial( d );
         res:=[ ];
         for k in [1,3..Length(cc)-1] do
             mon:=[ vars[1], d ];
             Append( mon, cc[k] );
             Add( res, mon ); Add( res, dfac*cc[k+1] );
         od;
         Append( res, WriteAsLCOfBinoms( vars, q ) );
         for k in [2,4..Length(res)] do
             if res[k] = 0*res[k] then
                 Unbind( res[k-1] ); Unbind( res[k] );
             fi;
         od;

         return Compacted( res );
     end;

    # We collect the UEALattice element represented by the data in `lst'.
    # `lst' represents a UEALattice element in the usual way, except that
    # a Cartan element is now not represented by an index, but by a list
    # of two elements: the element of the Cartan subalgebra, and an integer
    # (meaning `h-k', if the list is [h,k]). The ordering
    # is as follows: first come the `negative' root vectors (in the
    # same order as the roots), then the Cartan elements, and then the
    # `positive' root vectors.

    todo:= ShallowCopy( lst );
    dones:= [ ];

    while todo <> [] do

     # `collocc' will be `true' once a collection has occurred.

        collocc:= false;
        mon:= ShallowCopy(todo[1]);

     # We collect `mon'.

        i:= 1;
        while i <= Length( mon ) - 3 do

            # Collect `mon[i]' and `mon[i+1]'.
            if IsList( mon[i] ) and IsList( mon[i+2] ) then

                # They are both Cartan elements; so we do nothing.
                i:= i+2;
            elif IsList( mon[i] ) and not IsList( mon[i+2] ) then

                #`mon[i]' is a Cartan element, but `mon[i+2]' is not.
                if mon[i+2] > noPosR then

                    # They are in the right order; so we do nothing.
                    i:= i+2;
                else

                    # They are not in the right order, so we swap.
                    # `cf' is the coefficient in [ h, x ] = cf*x,
                    # where h is the Cartan element, x an element from the
                    # root space corresponding to `mon[i+2]'. When swapping
                    # the second element of the list representing the Cartan
                    # element changes.
                    cf:= Coefficients( BH, mon[i][1] )*posR[mon[i+2]];
                    temp:= mon[i];
                    temp[2]:= temp[2] +mon[i+3]*cf;
                    mon[i]:= mon[i+2];
                    mon[i+2]:= temp;

                    # Swap the coefficients.
                    temp:= mon[i+1];
                    mon[i+1]:= mon[i+3];
                    mon[i+3]:= temp;
                    todo[1]:= mon;
                    i:= 1;

                fi;
            elif not IsList( mon[i] ) and IsList( mon[i+2] ) then

                # Here `mon[i]' is no Cartan element, but `mon[i+2]' is. We
                # do the same as above.
                if mon[i] <= noPosR then
                    i:= i+2;
                else
                    cf:= Coefficients( BH, mon[i+2][1] )*posR[mon[i]];
                    temp:= mon[i+2];
                    temp[2]:= temp[2] - mon[i+1]*cf;
                    mon[i+2]:= mon[i];
                    mon[i]:= temp;
                    temp:= mon[i+1];
                    mon[i+1]:= mon[i+3];
                    mon[i+3]:= temp;
                    todo[1]:= mon;
                    i:= 1;
                fi;
            elif mon[i] = mon[i+2] then

                # They are the same; so we take them together. This costs
                # a binomial factor.
                mon[i+1]:= mon[i+1]+mon[i+3];
                todo[2]:= todo[2]*Binomial(mon[i+1],mon[i+3]);

                Remove( mon, i+2 );
                Remove( mon, i+2 );
                todo[1]:= mon;
            elif mon[i] < mon[i+2] then

                # They are in the right order; we do nothing.
                i:=i+2;
            else

                # We swap them. There are two cases: the two roots are
                # each others negatives, or not. In the first case we
                # get extra Cartan elements. In both cases the result of
                # swapping the two elements will be contained in `rr'.
                # To every element of `rr' we then have to prepend
                # `start' and to append `tail'.

                cf:= todo[2];
                Unbind( todo[1] ); Unbind( todo[2] );
                start:= mon{[1..i-1]};
                tail:= mon{[i+4..Length(mon)]};
                if posR[mon[i]] = -posR[mon[i+2]] then
                    i1:= mon[i]; j1:= mon[i+2];
                    m:= mon[i+1]; n:= mon[i+3];
                    h:= Rvecs[i1]*Rvecs[j1];
                    min:= Minimum( m, n );
                    rr:= [ ];
                    for k in [0..min] do
                      mon1:= [ ];
                      if n-k>0 then
                        Append( mon1, [ j1, n-k ] );
                      fi;
                      if k > 0 then
                        Append( mon1, [ [ h, -n-m+2*k ], k ] );
                      fi;
                      if m-k > 0 then
                        Append( mon1, [ i1, m-k ] );
                      fi;
                      Add( rr, mon1 ); Add( rr, 1 );
                    od;

                else

                # In the second case we have to swap two powers of root
                # vectors. According to the form of the root string
                # we distinguish a few cases. In each case we have a
                # different formula for the result.

                    i1:= mon[i]; j1:= mon[i+2];
                    m:= mon[i+1]; n:= mon[i+3];
                    a:= posR[j1]; b:= posR[i1];
                    if a+b in posR then
                       if a+2*b in posR then
                          if a+3*b in posR then
                             type:= "G2+";
                          else
                             if 2*a+b in posR then
                                type:= "G2~";
                             else
                                type := "B2+";
                             fi;
                          fi;
                       elif 2*a+b in posR then
                            if 3*a+b in posR then
                               type:= "G2-";
                            else
                               type:= "B2-";
                            fi;
                       else
                            type:= "A2";
                       fi;
                    else
                       type:= "A1A1";
                    fi;

                    rr:= [ ];
                    if type = "A1A1" then

                       # The elements simply commute.
                       rr:= [ [ j1, n, i1, m ], 1 ];
                    elif type = "A2" then

                       c1:= -RT[j1][i1];
                       c2:= 1;
                       p1:= Position( posR, a+b );
                       for k in [0..Minimum(m,n)] do
                          mon1:= [ ];
                          if n-k > 0 then
                             Append( mon1, [ j1, n-k ] );
                          fi;
                          if m-k > 0 then
                             Append( mon1, [ i1, m-k] );
                          fi;
                          if k>0 then
                             Append( mon1, [ p1, k ] );
                          fi;
                          Add( rr, mon1 );
                          Add( rr, c2 );
                          c2:= c2*c1;
                       od;

                    elif type = "B2+" then

                       c1:= -RT[j1][i1];
                       p1:= Position( posR, a+b );
                       p2:= Position( posR, a+2*b );
                       c2:= -c1*RT[i1][p1]/2;
                       min:= Minimum( m,n );
                       for k in [0..min] do
                          for l in [0..min] do
                             if n-k-l >= 0 and m-k-2*l >= 0 then
                                mon1:= [ ];
                                if n-k-l > 0 then
                                   Append( mon1, [ j1, n-k-l ] );
                                fi;
                                if m-k-2*l > 0 then
                                   Append( mon1, [ i1, m-k-2*l ] );
                                fi;
                                if k > 0 then
                                   Append( mon1, [ p1, k ] );
                                fi;
                                if l > 0 then
                                   Append( mon1, [ p2, l ] );
                                fi;
                                Add( rr, mon1 );
                                Add( rr, c1^k*c2^l );
                             fi;
                          od;
                       od;

                    elif type = "B2-" then

                       c1:= -RT[j1][i1];
                       p1:= Position( posR, a+b );
                       p2:= Position( posR, 2*a+b );
                       c2:= -c1*RT[j1][p1]/2;
                       min:= Minimum( m,n );
                       for k in [0..min] do
                          for l in [0..min] do
                             if n-k-2*l >= 0 and m-k-l >= 0 then
                                mon1:= [ ];
                                if n-k-2*l > 0 then
                                   Append( mon1, [ j1, n-k-2*l ] );
                                fi;
                                if m-k-l > 0 then
                                   Append( mon1, [ i1, m-k-l ] );
                                fi;
                                if k > 0 then
                                   Append( mon1, [ p1, k ] );
                                fi;
                                if l > 0 then
                                   Append( mon1, [ p2, l ] );
                                fi;
                                Add( rr, mon1 );
                                Add( rr, c1^k*c2^l );
                             fi;
                          od;
                       od;

                    elif type = "G2+" then

                       p1:= Position( posR, a+b );
                       p2:= Position( posR, a+2*b );
                       p3:= Position( posR, a+3*b );
                       p4:= Position( posR, 2*a+3*b );
                       c1:= RT[j1][i1];
                       c2:= RT[i1][p1];
                       c3:= RT[p1][p2];
                       c4:= RT[i1][p2]/2;
                       min:= Minimum(m,n);
                       for p in [0..min] do
                          for q in [0..min] do
                             for r in [0..min] do
                                for s in [0..min] do
                                   if n-p-q-r-2*s>=0 and
                                          m-p-2*q-3*r-3*s >=0  then
                                      mon1:= [ ];
                                      if n-p-q-r-2*s > 0 then
                                         Append( mon1, [ j1, n-p-q-r-2*s ] );
                                      fi;
                                      if m-p-2*q-3*r-3*s > 0 then
                                         Append( mon1, [i1,m-p-2*q-3*r-3*s]);
                                      fi;
                                      if p > 0 then
                                         Append( mon1, [ p1, p ] );
                                      fi;
                                      if q > 0 then
                                         Append( mon1, [ p2, q ] );
                                      fi;
                                      if r > 0 then
                                         Append( mon1, [ p3, r ] );
                                      fi;
                                      if s > 0 then
                                         Append( mon1, [ p4, s ] );
                                      fi;
                                      Add( rr, mon1 );
                                      Add( rr, (-1)^(p+r)*(1/3)^(s+r)*(1/2)^q*
                                        c1^(p+q+r+2*s)*c2^(q+r+s)*c3^s*c4^r );
                                   fi;
                                od;
                             od;
                          od;
                       od;
                    elif type = "G2-" then

                       p1:= Position( posR, a+b );
                       p2:= Position( posR, 2*a+b );
                       p3:= Position( posR, 3*a+b );
                       p4:= Position( posR, 3*a+2*b );
                       c1:= RT[j1][i1];
                       c2:= RT[j1][p1]/2;
                       c3:= RT[j1][p2]/3;
                       c4:= (c1*RT[p1][p2]+c3*RT[i1][p3])/2;
                       min:= Minimum(m,n);
                       for p in [0..min] do
                          for q in [0..min] do
                             for r in [0..min] do
                                for s in [0..min] do
                                   if n-p-2*q-3*r-3*s>=0 and
                                                m-p-q-r-2*s >=0 then
                                      mon1:= [ ];
                                      if n-p-2*q-3*r-3*s > 0 then
                                         Append( mon1,[j1, n-p-2*q-3*r-3*s]);
                                      fi;
                                      if m-p-q-r-2*s > 0 then
                                         Append( mon1, [ i1, m-p-q-r-2*s ] );
                                      fi;
                                      if p > 0 then
                                         Append( mon1, [ p1, p ] );
                                      fi;
                                      if q > 0 then
                                         Append( mon1, [ p2, q ] );
                                      fi;
                                      if r > 0 then
                                         Append( mon1, [ p3, r ] );
                                      fi;
                                      if s > 0 then
                                         Append( mon1, [ p4, s ] );
                                      fi;
                                      Add( rr, mon1 );
                                      Add( rr, (-1)^(p+r)*
                                        c1^(p+q+r+s)*c2^(q+r+s)*c3^r*c4^s );
                                   fi;
                                od;
                             od;
                          od;
                       od;
                    elif type = "G2~" then

                       p1:= Position( posR, a+b );
                       p2:= Position( posR, 2*a+b );
                       p3:= Position( posR, a+2*b );
                       c1:= RT[j1][i1];
                       c2:= RT[j1][p1]/2;
                       c3:= RT[i1][p1]/2;
                       min:= Minimum(m,n);
                       for p in [0..min] do
                          for q in [0..min] do
                             for r in [0..min] do
                                if n-p-2*q-r>=0 and m-p-q-2*r >=0 then
                                   mon1:= [ ];
                                   if n-p-2*q-r > 0 then
                                      Append( mon1, [ j1, n-p-2*q-r ] );
                                   fi;
                                   if m-p-q-2*r > 0 then
                                      Append( mon1, [ i1, m-p-q-2*r ] );
                                   fi;
                                   if p > 0 then
                                      Append( mon1, [ p1, p ] );
                                   fi;
                                   if q > 0 then
                                      Append( mon1, [ p2, q ] );
                                   fi;
                                   if r > 0 then
                                      Append( mon1, [ p3, r ] );
                                   fi;
                                   Add( rr, mon1 );
                                   Add( rr, (-1)^(p)*c1^(p+q+r)*c2^q*c3^r);
                                fi;
                             od;
                          od;
                       od;
                    fi;
                fi;  # End of the piece that swapped two elements, and
                     # produced `rr', which we now insert.

                for j in [1,3..Length(rr)-1] do
--> --------------------

--> maximum size reached

--> --------------------

Quellcode-Bibliothek lierep.gi Sprache: unbekannt

[ 0.75Quellennavigators Projekt ]