########################################################
##  rootsys.g     W.A.de Graaf      (2011)
########################################################
##  This is a GAP4 script 
##  
##  It comes as auxiliary material to the paper  "Angular 
##  momentum non conserving symmetries in bosonic models"  
##  by L.Fortunato, W.A. de Graaf, J. Phys. A: Math. Theor. 
##  44 (2011) 145206 (12pp) 
##
########################################################

rootsystem:=    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.
          facs0, num, fam, f, b, c, r; 

    # 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 DeterminantMat( KillingMatrix( Basis( L ) ) ) = Zero( F ) 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

        if Length(j) = 1 then
           Add( newB, j ); 
        else
           V:= Basis( VectorSpace( F, j, "basis" ), j );
           M:= List( j, x -> Coefficients( V, i*x ) );
           f:= CharacteristicPolynomial( F, M );
           facs:= Set(Factors( f ));
           num:= IndeterminateNumberOfUnivariateLaurentPolynomial(f);
           fam:= FamilyObj( f );

           facs0:= [ ];
           for l in facs do
               if Degree(l) = 1 then
                  Add( facs0, l );
               elif Degree(l) = 2 then # we just take square roots...
                  cf:= CoefficientsOfUnivariatePolynomial(l);
                  b:= cf[2]; c:= cf[1];
                  r:= (-b+Sqrt(b^2-4*c))/2;
                  Add( facs0, PolynomialByExtRep( fam, [ [], -r, [num,1], 1] ) );
                  r:= (-b-Sqrt(b^2-4*c))/2;
                  Add( facs0, PolynomialByExtRep( fam, [ [], -r, [num,1], 1] ) );
               else
                  return fail;
               fi;
           od;

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

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