/*
This file contains Magma functions for computing generators of the little Weyl group corresponding to the
grading of a simple Lie algebra, defined by a finite order automorphism. This finite order automorphism
is given by its Kac diagram. First we explain how we represent a Kac diagram. 

For inner automorphisms the Kac diagram is the extended Dynkin diagram, with labels 0 or 1. We use the Bourbaki
enumeration of the nodes, where the "extra" nodes comes first (at position 0, so to say). Using that we 
give the Kac diagram as a sequence of 0's and 1's.

Example in E6:
                          
                           0
                           |
                           1
                           |
Kac diagram:         1--0--0--1--0

Sequence:  [ 0, 1, 1, 0, 0, 1, 0 ]


For outer automorphisms we do the same, but in this case the diagram is the twisted
affine Dynkin diagram. For that we use the following enumerations of the nodes of the
diagram:

A_2^(2)       1 =>= 2  (more bonds, but difficult to draw...)

A_2l^(2)      1 =>= 2--3 ... --(l-1) =>= l

                 1
                 |
A_{2l-1}^(2)  2--3--4--....--(l-1) =<= l

D_l^(2)       1 =<= 2--3 ... --(l-1) =>= l

E_6^(2)       1--2--3=<=4--5

D_4^(3)       1--2=<=3 (again, more bonds needed...)


The main functions of this file are:

**  LittleWeylGroup( type, rank, diag ) **

Here diag is a Kac diagram of an *inner* automorphism of the simple Lie algebra of type and rank
given by the first two inputs. The output is a record with the following fields

  theta      : the automorphism with Kac diagram diag, 
  L          : the Lie algebra (i.e., the domain of theta), 
  c          : a sequence containing the basis vectors of the Cartan subspace,
  grading    : a sequence of sequences containing the basis vectors of the homogeneous components of L,
               corresponding to the grading induced by theta,
  generators : a sequence of complex reflections generating the little Weyl group,
  F          : the field (extension of Q) over which everything is defined,
  Wc         : the little Weyl group, as subgroup of GL( #c, F ).


Example:

> time r:= LittleWeylGroup( "E", 6, [0,0,0,0,1,0,0] );
Time: 50.760
> Wc:= r`Wc;
> Wc;
MatrixGroup(3, Cyclotomic Field of order 3 and degree 2)
Generators:
    [-w - 1      0      0]
    [     0      1      0]
    [     0      0      1]

    [  1/3*(-w + 1)    1/3*(w + 2) 1/3*(-2*w - 1)]
    [   1/3*(w + 2)   1/3*(-w + 1)  1/3*(2*w + 1)]
    [   1/3*(w - 1)   1/3*(-w + 1)   1/3*(-w + 1)]

    [1 0 0]
    [0 w 0]
    [0 0 1]
> R:= InvariantRing( Wc );
> FundamentalInvariants(R);
[
    x1^6 + 10*x1^3*x2^3 + 10*x1^3*x3^3 + x2^6 - 10*x2^3*x3^3 + x3^6,
    x1^6*x2^3 - x1^6*x3^3 + x1^3*x2^6 - x1^3*x3^6 + x2^6*x3^3 - x2^3*x3^6,
    x1^12 - 58*x1^9*x2^3 - 58*x1^9*x3^3 - 210*x1^6*x2^6 + 336*x1^6*x2^3*x3^3 - 210*x1^6*x3^6 - 58*x1^3*x2^9 - 
        336*x1^3*x2^6*x3^3 - 336*x1^3*x2^3*x3^6 - 58*x1^3*x3^9 + x2^12 + 58*x2^9*x3^3 - 210*x2^6*x3^6 + 58*x2^3*x3^9 + x3^12
]


**  ProbableLittleWeylGroup( type, rank, diag ) **

This is the same as the previous function, with two differences. It gives a group which contains
the little Weyl group, and in many cases (all cases we tried) is equal to it. (Also, if the Lie 
algebra is E8, F4, or G2 then the function does the same as the previous one.) Secondly, it also
works for outer automorphisms.


*/

BilinearFormMat:= function( C )

   rank:= NumberOfRows( C );
   
   while not IsSymmetric(C) do 

   d:= [ 1 : i in [1..rank ] ];
   for i in [1..rank] do
       for j in [i+1..rank] do
           if C[i][j] ne C[j][i] then
              if C[i][j] eq -1 then
                 d[j]:= -C[j][i];
              else
                 d[i]:= -C[i][j];
              end if;
           end if;
       end for;
   end for;

   B:= Transpose( C );
   for i in [1..rank] do
       B[i] *:= d[i];
   end for;

   C:= Transpose(B);
   end while;

   return C;

end function;

ip:= function( B, v, w )


    n:= NumberOfRows( B );
    sum:= 0;
    for i in [1..n] do
        for j in [1..n] do
            sum +:= v[i]*w[j]*B[i][j];
        end for;
    end for;

    return sum;

end function;

ExtendedCM:= function( R )

  // Here R is a root datum of crystallographic type

    posR:= PositiveRoots(R);
    hts:= [ &+Eltseq(u) : u in posR ];
    p:= Index( hts, Maximum(hts) );
    sim:= [-posR[p]];
    sim cat:= [ SimpleRoots(R)[i] : i in [1..Rank(R)] ];
    sim:= [ Vector( Rationals(), u ) : u in sim ];

    B:= BilinearFormMat( CartanMatrix(R) );
    C:= ZeroMatrix( Rationals(), #sim, #sim );
    for i in [1..#sim] do
        for j in [1..#sim] do
            C[i][j]:= 2*ip( B, sim[i], sim[j] )/ip( B, sim[j], sim[j] );
        end for;
    end for;

    v:= Eltseq( Kernel(C).1 );
    a:= Lcm( [ Denominator( v[i] ) : i in [1..#v] ] );
    v:= [ a*v[i] : i in [1..#v] ]; 

    return C, v, p;


end function;


make_homs:= function( L, g, imgs )

   // L a Lie algebra generated by g; 
   // imgs is a list of lists of the images of g, we make the 
   // corr homomorphism corr to each list.

   
   bas:= g;
   basim:= imgs;

   V:= RModule( BaseRing(L), Dimension(L) );
   U:= sub< V | [ Eltseq( u ) : u in bas ] >;

   newbas:= bas;
   newbasi:= basim;

   while #newbas gt 0 do

       newb:= [ ];
       newbi:= [ [ ] : i in [1..#imgs] ];
       for i in [1..#g] do
           if Dimension(U) eq Dimension(L) then break; end if;
           for j in [1..#newbas] do
               if Dimension(U) eq Dimension(L) then break; end if;
               u:= g[i]*newbas[j];
               if not V!Eltseq(u) in U then
                  Append( ~bas, u );
                  Append( ~newb, u );

                  for k in [1..#imgs] do
                     
                      w:= imgs[k][i]*newbasi[k][j];
                      Append( ~basim[k], w );
                      Append( ~newbi[k], w );
                  end for;
  
                  U:= sub< V | [ Eltseq(u) : u in bas ] >;
                  
               end if;
           end for;
       end for;
       newbas:= newb; newbasi:= newbi;
   end while;

   U:= RModuleWithBasis( [ V!Eltseq(u) : u in bas ] );
   auts:= [ ];
   for k in [1..#imgs] do
      ii:= [ ];
      for x in Basis(L) do
          cf:= Coordinates( U, V!Eltseq(x) );
          Append( ~ii, &+[ cf[i]*basim[k][i] : i in [1..#cf] ] );
      end for;
      Append( ~auts, hom< L -> L | ii > );
   end for;

   return auts;

end function;

make_hom0:= function( K, L, g, imgs )

   // L a Lie algebra generated by g; 
   // imgs is a list of lists of the images of g, we make the 
   // corr homomorphism corr to each list.

   
   bas:= g;
   basim:= imgs;

   V:= RModule( BaseRing(L), Dimension(L) );
   U:= sub< V | [ Eltseq( u ) : u in bas ] >;

   newbas:= bas;
   newbasi:= basim;

   while #newbas gt 0 do

       newb:= [ ];
       newbi:= [ [ ] : i in [1..#imgs] ];
       for i in [1..#g] do
           if Dimension(U) eq Dimension(L) then break; end if;
           for j in [1..#newbas] do
               if Dimension(U) eq Dimension(L) then break; end if;
               u:= g[i]*newbas[j];
               if not V!Eltseq(u) in U then
                  Append( ~bas, u );
                  Append( ~newb, u );

                  for k in [1..#imgs] do
                     
                      w:= imgs[k][i]*newbasi[k][j];
                      Append( ~basim[k], w );
                      Append( ~newbi[k], w );
                  end for;
  
                  U:= sub< V | [ Eltseq(u) : u in bas ] >;
                  
               end if;
           end for;
       end for;
       newbas:= newb; newbasi:= newbi;
   end while;

   U:= RModuleWithBasis( [ V!Eltseq(u) : u in bas ] );
   auts:= [ ];
   for k in [1..#imgs] do
      ii:= [ ];
      for x in Basis(L) do
          cf:= Coordinates( U, V!Eltseq(x) );
          Append( ~ii, &+[ cf[i]*basim[k][i] : i in [1..#cf] ] );
      end for;
      Append( ~auts, hom< K -> L | ii > );
   end for;

   return auts;

end function;


FiniteOrderInnerAutomorphism:= function( type, diag )


    // finite order auts of the simple Lie algebra, 
    // of order m, that correspond to untwisted diagram diag.

    rank:= #diag-1;
    //L:= SimpleLieAlgebra( type, rank, F );
    R:= RootDatum(type cat IntegerToString(rank));
    C,a,p:= ExtendedCM( R );
    m:= &+[ a[i]*diag[i] : i in [1..rank+1] ];
    m:= Integers()!m;

    F<w>:= CyclotomicField( m );
    L:= LieAlgebra( R, F );

    x,y,h:= ChevalleyBasis(L);
    g:= [ y[p] ];
    g cat:= [ x[i] : i in [1..rank] ];

    ss:= [ diag ];

    g0:= [ [ w^s[i]*g[i] : i in [1..rank+1] ] : s in ss ];
    return make_homs( L, g, g0 )[1],m;


end function;

FiniteOrderOuterAutomorphism:= function( type, rank, d, m, diag )


     // corresponding to the diagram automorphism of 
     // order d.
     // m is the order of the automorphism to be.

     phi:= function( rt )


        if type eq "A"  then
           return Reverse(rt);
        elif type eq "D" and d eq 2 then
           r0:= rt;
           r0[ rank ]:= rt[rank-1];
           r0[rank-1]:= rt[rank];
        elif type eq "D" and d eq 3 then
           r0:= rt;
           r0[1]:= rt[4]; r0[3]:= rt[1]; r0[4]:= rt[3];
        else  //i.e. E6
           r0:= rt;
           r0[1]:= rt[6]; r0[6]:= rt[1];
           r0[3]:=rt[5]; r0[5]:= rt[3];
        end if;
        return r0; 

     end function;

    R:= RootDatum( type cat IntegerToString(rank) );
    if d eq 2 then
       F<w>:= CyclotomicField( m );
       w0:= w;  
       L:= LieAlgebra( R, F );
    else
       F<w>:= CyclotomicField( m );
       w0:= w;
       w3:= w^(Integers()!(m/3));  
       L:= LieAlgebra( R, F );
    end if;
    
    x,y,h:= CanonicalGenerators( L );
    cg0:= [ [], [], [] ];

    idmat:= IdentityMatrix( Integers(), rank );
    sim:= [ Eltseq(idmat[i]) : i in [1..NumberOfRows(idmat)] ];
    for i in [1..#sim] do
        pos:= Index( sim, phi(sim[i]) );
        Append( ~cg0[1], x[pos] );
        Append( ~cg0[2], y[pos] );
        Append( ~cg0[3], h[pos] );     
    end for;

    g0:= x; g0 cat:= y; g0 cat:= h;
    g1:= cg0[1]; g1 cat:= cg0[2]; g1 cat:= cg0[3];
    f:= make_homs( L, g0, [g1] )[1];

    mat:= [ ];
    for i in [1..Dimension(L)] do
        Append( ~mat, Coordinates( L, f(Basis(L)[i]) ) );
    end for;

    mat:= Matrix( F, mat );
    mat0:= mat- IdentityMatrix( F, Dimension(L) );

    sol:= [ Eltseq(u) : u in Basis( Nullspace( mat0 ) ) ];
       
    K:= sub< L | [ L!u : u in sol ] >; 

    xK,yK,hK:= CanonicalGenerators( K );

    if d eq 2 then
       mat0:= mat+IdentityMatrix( F, Dimension(L) );
    else
       mat0:= mat-w3*IdentityMatrix( F, Dimension(L) );
    end if;

    sol:= [ Eltseq(u) : u in Basis( Nullspace( mat0 ) ) ];
    V:= RModule( F, Dimension(L) );
    W:= sub< V | [ V!u : u in sol ] >;
    // get simultaneous kernel of all elements of yK acting on W...
    mat:= [ ];
    for u in yK do
        mm:= [ ];
        for v in Basis(W) do
            Append( ~mm, Coordinates( W, V!(Eltseq(u*L!Eltseq(V!v))) ) );
        end for;
        mm:= Transpose( Matrix( F, mm ) );
        for i in [1..NumberOfRows(mm)] do
            Append( ~mat, Eltseq( mm[i] ) );
        end for;
    end for;

    mat:= Transpose( Matrix( F, mat ) );
    sol:= Nullspace( mat );

    // then x0 will be the "lowest" root vector (see thm in Helgasson).
    x0:= L!Eltseq(V!(W!Eltseq( Basis(sol)[1] )));

    // next we get the "root" corresponding to x0:

    sp:= RModuleWithBasis( [ V!Eltseq(x0) ] );

    rt:= [ ];
    for h in hK do
        Append( ~rt, Coordinates( sp, V!(Eltseq(h*x0) ) )[1] );
    end for;

    // This now is the root wrt the "weight" basis...

    // now get the extended Cartan matrix...

    U:= RModule( F, #rt );
    _,_,_,C:= RootSystem(K);
    sim:= RModuleWithBasis( [ U!Eltseq( C[i] ) : i in [1..NumberOfRows(C)]]);
    rts:= [ Eltseq( Coordinates( sim, U!rt ) ) ];
    idmat:= IdentityMatrix( F, #rt );
    rts cat:= [ Eltseq(idmat[i]) : i in [1..NumberOfRows(idmat)] ];

    B:= BilinearFormMat( C );

    extC:= ZeroMatrix( Integers(), #rts, #rts );
    for i in [1..#rts] do
        for j in [1..#rts] do
            extC[i][j]:= Integers()!( 2*( ip(B,rts[i],rts[j]) )/
                                            ip(B,rts[j],rts[j]) );
        end for;
    end for;

    v:= Basis( Nullspace(Matrix(Rationals(), extC )) )[1];
    vv:= Eltseq(v);
    a:= Lcm( [ Denominator(u) : u in vv ] );
    v:= a*v;

    g:= [ x0 ];
    g cat:= xK;
//print "subalg test....", Dimension( sub< L | g > ) eq Dimension(L);

    g0:= [ [ w0^diag[i]*g[i] : i in [1..#diag] ] ];
    return make_homs( L, g, g0 )[1];

end function;


matrix:= function( h )

    // here h is an automorphism; we compute its matrix
    // (using row convention)

    L:= Domain(h);
    return Matrix( BaseRing(L), [ Coordinates( L, h(L.i) ) : i in 
                                                     [1..Dimension(L)] ] );

end function;


GradingThetaRep:= function( h, m, w )

   // Here h is a finite-order automorphism, of order m;
   // we compute the corr grading, i.e., the eigenspaces corr to
   // w^i, i=0..m-1

   L:= Domain(h);
   mat:= matrix( h );

   sp:= [ ];
   for i in [0..m-1] do
       K:= Kernel( mat-w^i*IdentityMatrix( BaseRing(L), Dimension(L) ) );
       Append( ~sp, [ L!K.i : i in [1..Dimension(K)] ] );
   end for;

   return sp;

end function;

ss_part:= function( mat )


  f:= CharacteristicPolynomial( mat );
  g:= f/Gcd( f, Derivative( f ) );
  g:= Parent(f)!g;

  if f eq g then return mat; end if;

  _,_,w:= ExtendedGreatestCommonDivisor( g, Derivative( g ) );
  w:= w*g;
  B:= mat;
  zero:= 0*B;
  while Evaluate( g, B ) ne 0*B do
    B:= B - Evaluate( w, B );
  end while;

  return B;

end function;

is_nilpotent:= function( K, x )


    b:= [ x*K.i : i in [1..Dimension(K) ] ];
    for j in [1..Dimension(K)] do
        b:= [ x*b[i] : i in [1..Dimension(K)] ];
        if &and[ IsZero(u) : u in b ] then return true; end if;
    end for;

    return &and[ IsZero(u) : u in b ];


end function;

is_nilpotent0:= function( K, x )

  n:= Dimension(K);
  V:= RModule( BaseRing(K), n );
  W:= sub< V | [ V!Eltseq(K.i) : i in [1..n] ] >;
  while true do
     b:= [ x*K!( Eltseq(V!u) ) : u in Basis(W) ];
     W0:= sub< V | [ V!Eltseq(u) : u in b ] >;
     if Dimension(W0) eq Dimension(W) then
        return false;
     elif Dimension(W0) eq 0 then
        return true;
     else
        W:= W0;
     end if;
  end while;

end function;

CartanSubspace:= function( L, gr0, gr1 )


     F:= BaseRing(L);

     s:= L*L;
     r:= Centre(L);

     //make s and r subspaces of L, using vector spaces
     LL:= RModule( BaseRing(L), Dimension(L) );

     ss:= sub< LL | [ LL!Eltseq(L!s.i) : i in [1..Dimension(s)] ] >;

     gg0:= sub< LL | [ LL!Eltseq(u) : u in gr0 ] >;
     gg1:= sub< LL | [ LL!Eltseq(u) : u in gr1 ] >;

     ss1:= ss meet gg1;

     if Dimension(ss1) eq 0 then
        rr:= sub< LL | [ LL!Eltseq( L!u ) : u in Basis(r) ] >;
        rr1:= rr meet gg1;
        return [ L!Eltseq( LL!u ) : u in Basis(rr1) ];
     end if;

     bas1:= [ L!Eltseq( LL!u ) : u in Basis( ss1 ) ];

     Omega:= [0,1];

     dims:= Dimension(s);
     M:= RMatrixSpace( BaseRing(L), dims, dims );
     ad_mats:= [ AdjointMatrix( s, s.i ) : i in [1..dims]];
     basAd:= RMatrixSpaceWithBasis( [ M!u : u in ad_mats ] );



     jord_dec:= function( y ) // y in s.. find Jordan dec...

         if is_nilpotent0( s, y ) then
            return Zero(s);
         end if;
         a:= Matrix( AdjointMatrix( s, y ) );
         a1:= ss_part( Matrix(a) );
         cfs:= Coordinates( basAd, M!a1 );
         return &+[ cfs[i]*s.i : i in [1..#cfs] ];
     end function;

     while true do 

        // first we try to find an "easy" non nilpotent element;
        // if that doesn't exist then we try the random way... 
        found:= false;
/*        for i in [1..#bas1] do
            x:= s!bas1[i];
            sp:= jord_dec(x);
            if not IsZero(sp) then found:= true; break i; end if;
        end for; */

        if not found then
           cfx:= [ Random(Omega) : i in [1..#bas1] ];
           x:= s!(&+[ cfx[i]*bas1[i] : i in [1..#bas1] ]);

           b:= is_nilpotent0( s, x );

           //cfs:= Coefficients(f);
           //if &and[ IsZero( cfs[i] ) : i in [1..#cfs-1] ] then
           if b then
              sp:= Zero(s);
           else
              for i in [1..#bas1] do
                  k:= 0;
                  found:= false;
                  while not found do
                     cfx[i]:= k;
                     x:= s!(&+[ cfx[i]*bas1[i] : i in [1..#bas1] ]);
                     if not is_nilpotent0( s, x ) then
                        found:= true;
                     else
                        k:= k+1;
                     end if;
                  end while;
              end for;
              sp:= jord_dec( x );
           end if;
        end if;

        if not IsZero(sp) then
           
           K:= Centralizer( L, sub< L | [L!sp] > );
           KK:= sub< LL | [ LL!Eltseq(L!u) : u in Basis(K) ] >;
           U:= KK meet gg0;
           gK0:= [ K!( L!Eltseq( LL!u ) ) : u in Basis(U) ];
           U:= KK meet gg1;
           gK1:= [ K!( L!Eltseq( LL!u ) ) : u in Basis(U) ];
     
           cK:= $$( K, gK0, gK1 );
           return [ L!u : u in cK ];
        else

           ss0:= ss meet gg0;
           bas0:= [ L!Eltseq( LL!u ) : u in Basis( ss0 ) ];
           U:= sub< LL | [ LL!Eltseq( u*x ) : u in bas0 ] >;

           if Dimension(U) eq Dimension(ss1) then
              rr:= sub< LL | [ LL!Eltseq( L!u ) : u in Basis(r) ] >;
              rr1:= rr meet gg1;
              return [ L!Eltseq( LL!u ) : u in Basis(rr1) ];
           end if;

        end if;
     end while;


end function;

find_cart:= function( L, gr0, gr1, bound, deg )

    // find a number (bound) of Cartan subspaces, 
    // together with corresponding CSA's and splitting fields.

    cc:= []; HH:=[]; FF:= [];

    K0:= sub< L | gr0 >;
    for i in [1..bound] do

       c:= CartanSubspace( L, gr0, gr1 );
       K:= K0 meet Centralizer(L,sub< L | c >);

       b:= c cat [ L!u : u in Basis( CartanSubalgebra(K) )];
       H:= Centralizer( L, sub< L | b > );

       if Dimension(H) eq Dimension( CartanSubalgebra(L) ) then

          facs:= [  ];
          for x in Basis(H) do
              ad:= AdjointMatrix( L, L!x );
              facs cat:= [ u[1] : u in 
                    Factorization(CharacteristicPolynomial(ad)) ];
          end for;
          f0:= [ ];

          for f in facs do
              if not f in f0 then Append( ~f0, f ); end if;
          end for;
          facs:= f0;

          F:= BaseField(L);
          P:= PolynomialRing( F );
          for f in facs do
              E:= SplittingField( P!f );
              if E ne F then
                 F:= E; P:= PolynomialRing(F);
              end if;
              if Degree(E) ge deg then return [], [], []; end if;
          end for;

          Append( ~cc, c );
          Append( ~HH, H );
          Append( ~FF, F );

       end if;

    end for;

    return cc, HH, FF;

end function;



rootOrder := function( u, v )
  n := #Eltseq(u);
  height := func< v | &+Eltseq(v) >;

  h := height(u);  k := height(v);
  if h ne k then
    if   h gt 0 and k lt 0 then return -1;
    elif h lt 0 and k gt 0 then return 1;
    elif h gt 0 and k gt 0 then return h-k;
    else return k-h;
    end if;
  else
    if h eq 0 then return 0; end if;
    b := exists(i){i : i in [n..1 by -1] | u[i] ne v[i]};
    if u[i] gt 0 or v[i] gt 0 then
      return u[i]-v[i];
    else
      return v[i]-u[i];
    end if;
  end if;
end function;

isPara := function( u, v )
  if exists(i){ i : i in [1..OverDimension(u)] | u[i] ne 0 } then
    return v[i]*u eq u[i]*v;
  else  // u eq 0
    return true;
  end if;
end function;


integerConvert := function( A );
  p := Characteristic( BaseRing( Parent(A) ) );
  if p eq 0 then
    return Vector( Integers(), [ Integers()!a : a in Eltseq(A) ] );
  else
    A := ChangeRing( A, Integers() );
    P := Parent( A );
    return P![ (2*a le p) select a else a-p : a in Eltseq(A) ];
  end if;
end function;

// the action of ad(h) on B
AdOnSubbasis := function( L, h, B )
    F := CoefficientField( L );
    V:= VectorSpace( F, Dimension( L ) );
    W:= VectorSpaceWithBasis([ V |  V!Eltseq(b) :  b in B ]);
    M:= [];
    for j in [1..#B] do
      M cat:= Coordinates( W, V!Eltseq( L!(h*B[j]) ) );
    end for;
    return MatrixAlgebra( F, #B ) ! M;
end function;



function fund_rts_from_pos_rts(posR)
    // A positive root is called fundamental if it is not the sum of two
    // OTHER positive roots. We calculate the set of fundamental roots.
    // removed OTHER from the code to allow nonreduced systems
    
    fundR:= [ ];
    for a in posR do
      issum:= false;
      for i in [1..#posR] do
        for j in [i..#posR] do
          if a eq posR[i]+posR[j] then
            issum:=true; break;
          end if;
        end for;
        if issum then break; end if;
      end for;
      if not issum then
        Append( ~fundR, a );
      end if;
    end for;
    return fundR;
end function;

// 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 )
  s:=0; t:=0;
  rt:=a-b;
  while (rt in R) or (rt eq 0*R[1]) do
    rt:=rt-b;
    s:=s+1;
  end while;
  rt:=a+b;
  while (rt in R) or (rt eq 0*R[1]) do
    rt:=rt+b;
    t:=t+1;
  end while;
  return s-t; //a eq b select s+t else s-t;
end function;

ratrd:= function( L, H )


    // 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.

    F := CoefficientField( L );
    V := VectorSpace( F, Dimension( L ) );

    rvs := [ [ L.i : i in [1..Dimension( L )] ] ];  
    rts := [ [] ];

    for i in [1..Dimension(H)] do

      newrvs:= [ ];  newrts := [ ];
      for j in [1..#rvs] do
        M := AdOnSubbasis( L, H.i, rvs[j] );
        fs := Factorization( MinimalPolynomial( M ) );

        for f in fs do
          B := Basis( NullSpace( Evaluate( f[1], M ) ) );
          Append( ~newrvs, [ &+[ b[l]*rvs[j][l] : l in [1..#rvs[j]] ] : b in B ]);
          Append( ~newrts, rts[j] cat [f[1]] );
        end for;

      end for;
      rvs:= newrvs;  rts := newrts;

    end for;

    // Now we throw away the subspace H.
    _ := exists(i){ i : i in [1..#rvs] | rvs[i][1] in H };
    Remove( ~rvs, i );  Remove( ~rts, i );

    return rvs, rts;
end function;

ratrd0:= function( L, H )

    found:= false;
    while not found do 
        cf:= [ Random([-5..5]) : i in [1..Dimension(H)] ];
        h:= L!(&+[ cf[i]*H.i : i in [1..Dimension(H)] ] );
        M:= AdjointMatrix( L, h ); 
        f:= CharacteristicPolynomial( M );
        as := Roots( f ); 

        if #as eq Dimension(L)-Dimension(H)+1 then
           found:= true;
        end if;
    end while;

/*    for i in [1..Dimension(H)] do
        k:= 0;
        found:= false;
        while not found do
             cf[i]:= k;
             h:= L!(&+[ cf[i]*H.i : i in [1..Dimension(H)] ] );
             M:= AdjointMatrix( L, h ); 
             f:= MinimalPolynomial( M );
             if Degree( f ) eq Dimension(L)-Dimension(H)+1 then
                found:= true;
             else
                k:= k+1;
             end if;
        end while;
    end for;
"found2"; */

    sp:= [ ];
    for a in as do
        B := Basis( Eigenspace( M, a[1] ) );
        Append( ~sp, [ &+[ b[k]*L.k : k in [1..Dimension(L)] ] : b in B ]);
    end for;

    for i in [1..#sp] do 
        if sp[i][1] in H then k:= i; break i; end if;
    end for;
    Remove( ~sp, k );

    return sp;

end function;

rootdec:= function(L,H)

    k := CoefficientField( L );
    ratrvs:= ratrd( L, H );
    n := Dimension( H );

    MultipleFields:= false;
    
    if not MultipleFields then
      K := k;  LK := L;
      HK := sub< LK | [ LK | Vector( L!h ) : h in Basis( H ) ] >;
      V := VectorSpace( K, Dimension(L) );
      rtV := VectorSpace( K, n );
    end if;
    
    allrvs := MultipleFields select [**] else [];  
    allrts := MultipleFields select [**] else [rtV|];
    for idx in [1..#ratrvs] do
    
      rvs := [[ LK | Vector( x ) : x in ratrvs[idx] ]];  rts := [[]];
      for i in [1..Dimension(H)] do
        newrvs := [];  newrts := [];
        for j in [1..#rvs] do
          M := AdOnSubbasis( LK, HK.i, rvs[j] );
          as := Roots( MinimalPolynomial( M ) );
          for a in as do
            B := Basis( Eigenspace( M, a[1] ) );
            Append( ~newrvs, [ &+[ b[l]*rvs[j][l] : l in [1..#rvs[j]] ] : b in B ]);
            Append( ~newrts, rts[j] cat [a[1]] );
          end for;
        end for;
        rvs := newrvs;  rts := newrts;
      end for;
      
      for i in [1..#rvs] do
        Append( ~allrvs, rvs[i] );  Append( ~allrts, rtV!rts[i] );
      end for;
    end for;

    return allrvs, allrts;
end function;


rootsys:= function( L, H )

    k := BaseRing( L );
    B, S := rootdec( L, H );

    rootSp := VectorSpace( k, Dimension(H) );

    Rvecs := [L| B[i][1] : i in [1..#B] ];

    // A set of roots basR is calculated such that the set
    // { [ x_r, x_{-r} ] + Z | r in basR } is a basis of H/Z.
    basH := [  ];
    basR := [  ];
    K := Centre( L );
    i := 1;
    while K ne H and i le #S do
      a:= S[i];
      j:= Position( S, -a );
      h:= B[i][1]*B[j][1];
      if h notin K then
        Append( ~basR, a );  Append( ~basH, h );
        K := sub< H | K, h >;
      end if;
      i:=i+1;
    end while;

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

    posR:= [ ];
    i:=1;
    while #posR  lt #S div 2 do
      a:= S[i];
      if (not a in posR) and (not -a in posR) then
        cf:= k!0;
        j:= 0;
        while cf eq 0 do
          j:= j+1;
          cf:= CartInt( S, a, basR[j] );
        end while;
        if 0 lt cf then
          Append( ~posR, a );
        else
          Append( ~posR, -a );
        end if;
      end if;
      i:=i+1;
    end while;

    // compute fundamental roots
    fundR:= fund_rts_from_pos_rts(posR);

    // Now we calculate the Cartan matrix
    n := #fundR;
    C := Matrix( Integers(), n,n, [ [ CartInt( S, fundR[i], fundR[j] ) : 
      j in [1..n] ] : i in [1..n] ]);
    if not IsCartanMatrix( C ) then
      "This should not happen. Please email magma-bugs@maths.usyd.edu.au";
    end if;
    // We sort the roots in the same order as used in RootDtm.m
    // to prevent inconsistencies between the root datum and the
    // root system

    fundR := [ rootSp | Eltseq( v ) : v in fundR ];
    N := #posR;
    perm := (N ne 0) select Sym(N)!1 else Sym(1)!1;
        
    if IsIndependent( fundR ) then

      W := VectorSpaceWithBasis( fundR );
      posRint := [ integerConvert( Vector( Coordinates( W, v ) ) ) : v in posR ]
;
      Sort( ~posRint, rootOrder, ~perm );
      
    else  // this alternative method is needed over certain characteristics
    
      R := RootDatum( C );
      idxs := [ Position( posR, v ) : v in fundR ];
      posRvecs := Rvecs[[ Position(S,v) : v in posR ]];
      VL := VectorSpace( k, Dimension(L) );
      for r in [n+1..N] do
        e1,e2 := ExtraspecialPair( R, r );
        v := posRvecs[idxs[e1]] * posRvecs[idxs[e2]];
        ok := exists(s){ s : s in [1..N] | isPara(VL!v,VL!posRvecs[s]) };
        if not ok then
          error "Error: This should not happen. Please email magma-bugs@maths.usyd.edu.au";
        end if;
        idxs[r] := s;
      end for;
      perm := Parent(perm)!idxs;
      
    end if;
 
    // apply the permutation to the roots
    sorR := [ posR[ i^perm ] : i in [1..N] ];
    //sorR cat:= [ -sorR[i] : i in [1..#sorR] ];
 
    // We also sort the root vectors.
    PRvecs:= [ Rvecs[ Position( S, sorR[i] ) ] : i in [1..#sorR] ];
    NRvecs:= [ Rvecs[ Position( S, -sorR[i] ) ] : i in [1..#sorR] ];

    // And the fundamental roots
    fundR := sorR[[1..#fundR]];
    
    return sorR, PRvecs, NRvecs, fundR, C;
end function;


cangens:= function( posRv, negRv, rank )

    x:= posRv[[1..rank]];
    y:= negRv[[1..rank]];
    for i in [1..rank] do
        V:= VectorSpaceWithBasis( Matrix( [ x[i] ]  ));
        c:= Coordinates( V, V!( (x[i]*y[i])*x[i] ) )[1];
        y[i]:= y[i]*2/c;
    end for;
    
    return x, y, [ x[i]*y[i] : i in [1..rank] ];
    
end function;

findperm:= function( C, C1 )

    // first we need to find an enumeration of C1 so
    // that Cartan matrices become equal.
    // For the moment we do this in the most inefficient way!

    rank:= NumberOfRows(C);
    C0:= C1;
    Sn:= Sym(rank);
    pp:= [ u : u in Sn ];
    for p in pp do
        for i in [1..rank] do
            for j in [1..rank] do
                C0[i][j]:= C1[i^p][j^p];
            end for;
        end for;

        if C0 eq C then
           return p;
        end if;
    end for;

end function;


irrep:= function( L, type, rank, hw, C, h, x, y )


   // We construct an irrep of L of highest weight hw, 
   // such that the elements of h act diagonally, 
   // here hw is a list of highest weights, we construct the direct sum...
   // L may be defined over some number field, so
   // we make a "standard form" S over the same field,
   // and an isom L -> S. Such that h is mapped to the CSA of S
   //   Then we construct the rep mats of S,
   // and translate to L...
   // C is the Cartan matrix rel to h, x, y

   R:= RootDatum(type cat IntegerToString(rank));
   K:= LieAlgebra( R, Rationals() );
   S:= LieAlgebra( R, BaseField(L) );

   p:= findperm( CartanMatrix(R), C );
   gL:= [ x[i^p] : i in [1..rank] ];
   gL cat:= [ y[i^p] : i in [1..rank] ];
   gL cat:= [ h[i^p] : i in [1..rank] ];

   x1,y1,h1:= CanonicalGenerators(S);
   gS:= x1 cat y1 cat h1;
   f:= make_hom0( L, S, gL, [gS] );
   f:= f[1];

   F:= BaseRing(L);

   if #hw eq 1 then
      V:= HighestWeightModule( K, hw[1] );
      m:= [ ActionMatrix( V, K.i ) : i in [1..Dimension(K)] ];
      m:= [ Matrix( F, u ) : u in m ];
   else
      V1:= HighestWeightModule( K, hw[1] );
      m1:= [ ActionMatrix( V1, K.i ) : i in [1..Dimension(K)] ];
      V2:= HighestWeightModule( K, hw[2] );
      m2:= [ ActionMatrix( V2, K.i ) : i in [1..Dimension(K)] ];
      m:= [ ];
      n1:= Dimension(V1); n2:= Dimension( V2 );
      for k in [1..Dimension(K)] do
          m0:= ZeroMatrix( F, n1+n2, n1+n2 );
          for i in [1..n1] do 
              for j in [1..n1] do
                  m0[i][j]:= m1[k][i][j];
              end for;
          end for;
          for i in [1..n2] do 
              for j in [1..n2] do
                  m0[n1+i][n1+j]:= m2[k][i][j];
              end for;
          end for;
          Append( ~m, m0 );
      end for;
   end if;

   mats:= [ ];
   for i in [1..Dimension(L)] do
       c:= Eltseq(f(L.i));
       Append( ~mats, &+[ c[k]*m[k] : k in [1..#m] ] );
   end for;

   return function(x)
             cf:= Eltseq(x);
             return &+[ cf[i]*mats[i] : i in [1..#mats]];
          end function; 


end function;

weylcenelt:= function( L, x, y, h )

    // x,y,h: can gens


    rank:= #x;
    F:= BaseRing(L);
    R:= RationalFunctionField( F, 2*rank );
    names:= [ ];
    for i in [1..rank] do
        Append( ~names, "x" cat IntegerToString(i) );
    end for;
    for i in [1..rank] do
        Append( ~names, "y" cat IntegerToString(i) );
    end for;

    AssignNames( ~R, names );

    X:= x;
    last:= x;
    rts:= [ ];
    for i in [1..rank] do
        rt:= [ 0 : k in [1..rank] ];
        rt[i]:= 1;
        Append( ~rts, rt );
    end for;
    
    cfx:= [ R.i : i in [1..rank] ];
    cfs:= cfx;
    done:= false;
    while not done do

          newlast:= [ ];
          newrts:= [ ];
          newcfs:= [ ];
          for i in [1..#last] do
              for j in [1..#x] do
                  if not IsZero(x[j]*last[i]) then
                     rt:= rts[i]; rt[j] +:=1;
                     if not rt in newrts then
                        Append( ~newrts, rt );
                        Append( ~newlast, x[j]*last[i] );
                        Append( ~newcfs, R.j*cfs[i] );
                     end if;
                  end if;
              end for;
          end for;
          if #newlast eq 0 then
             done:= true;
          else
             X cat:= newlast;
             cfx cat:= newcfs;
             cfs:= newcfs;
             last:= newlast;
             rts:= newrts;
          end if;

    end while; 

    Y:= y;
    last:= y;
    rts:= [ ];
    for i in [1..rank] do
        rt:= [ 0 : k in [1..rank] ];
        rt[i]:= 1;
        Append( ~rts, rt );
    end for;

    //cfy:= [ R.i : i in [rank+1..2*rank] ];
    cfy:= [ R.i^-1 : i in [1..rank] ];
    cfs:= cfy;
    done:= false;
    while not done do

          newlast:= [ ];
          newrts:= [ ];
          newcfs:= [ ];
          for i in [1..#last] do
              for j in [1..#y] do
                  if not IsZero(y[j]*last[i]) then
                     rt:= rts[i]; rt[j] +:=1;
                     if not rt in newrts then
                        Append( ~newrts, rt );
                        Append( ~newlast, y[j]*last[i] );
                        Append( ~newcfs, R.j^-1*cfs[i] );
//Append( ~newcfs, R.(rank+j)*cfs[i] );
                     end if;
                  end if;
              end for;
          end for;
          if #newlast eq 0 then
             done:= true;
          else
             Y cat:= newlast;
             cfy cat:= newcfs;
             cfs:= newcfs;
             last:= newlast;
             rts:= newrts;
          end if;

    end while; 

    bas:= h cat X cat Y;
    bas:= [ Vector( R, Eltseq(u) ) : u in bas ];    

//sp:= VectorSpace( R, Dimension(L) );
//Dimension( sub< sp | [ sp!u : u in bas ] > );

    sp:= VectorSpaceWithBasis( bas );
    cfs:= [ One(R) : i in [1..rank]];
    cfs cat:= cfx; cfs cat:= cfy;

    return R, function( elm )

              c:= Coordinates( sp, Vector( R, Eltseq( elm )) );
              c0:= [ c[i]*cfs[i] : i in [1..#c] ];
              return &+[ c0[i]*bas[i] : i in [1..#c0] ];
           end function, 
           [ Numerator(R.i)*Numerator(R.(rank+i))-1 : i in [1..rank] ];

end function;

/*liftaut:= function( L, F, p, f, x, y )

    // here p is a rep of L, eg 27 dim E6 rep
    // f is an automorphism of L, we try to lift it 
    // into the group corr to the image of L in gl(V)

    // returned is a list of matrices (usually has only one elt), a such that
    // rho(f(x)) = a*rho(x)*a^-1, where rho is the corr rep.


    n:= NumberOfRows( p(L.1) );
    g:= x cat y;

    eqns:= [ ];
    for a in g do

        m1:= p(a);
        m2:= p( f(a) );

        for k in [1..n] do 
            for l in [1..n] do

                eqn:= [ F!0 : i in [1..n^2] ];
                for i in [1..n] do
                    eqn[i+(k-1)*n]:= eqn[i+(k-1)*n] + m1[i][l];
                    eqn[l+(i-1)*n]:= eqn[l+(i-1)*n] - m2[k][i];
                end for;
                Append( ~eqns, eqn );
            end for;
        end for;

    end for;
    
    sol:= Kernel( Transpose( Matrix( F, eqns ) ) );
    bas:= [ ];
    for a in Basis(sol) do

        A:= ZeroMatrix( F, n, n );
        for k in [1..n] do
            for l in [1..n] do
                A[k][l]:= a[l+(k-1)*n];
            end for;
        end for;
        Append( ~bas, A );
    end for;
    return bas;

end function;*/

liftaut_diag:= function( L, F, p, f, x, y )

    // here p is a rep of L, eg 27 dim E6 rep
    // f is an automorphism of L, we try to lift it 
    // into the group corr to the image of L in gl(V)
     
    // we know that the resulting matrix has to be diagonal...

    // returned is a list of matrices (usually has only one elt), a such that
    // rho(f(x)) = a*rho(x)*a^-1, where rho is the corr rep.


    n:= NumberOfRows( p(L.1) );
    g:= x cat y;

    eqns:= [ ];
    for a in g do

        m1:= p(a);
        m2:= p( f(a) );

        for k in [1..n] do
            for l in [1..n] do

                eqn:= [ F!0 : i in [1..n] ];
                eqn[k]:= eqn[k] + m1[k][l];
                eqn[l]:= eqn[l] - m2[k][l];
                Append( ~eqns, eqn );
            end for;
        end for;

    end for;
    
    sol:= Kernel( Transpose( Matrix( F, eqns ) ) );
    bas:= [ ];
    for a in Basis(sol) do

        A:= ZeroMatrix( F, n, n );
        for k in [1..n] do
            A[k][k]:= a[k];
        end for;
        Append( ~bas, A );
    end for;
    return bas;

end function;


/*checkcan:= function( x, y, h, C )

    rank:= #x;
    if not &and[ x[i]*y[i] eq h[i] : i in [1..rank] ] then
1;
       return false;
    end if;

    if not &and(Flat( [ [ h[i]*x[j] eq C[j,i]*x[j] : i in [1..rank] ] : 
              j in [1..rank]] ) )   then
2;
       return false;
    end if;
 
    if not &and(Flat( [ [ h[i]*y[j] eq -C[j,i]*y[j] : i in [1..rank] ] : 
              j in [1..rank]] ) )   then
3;
       return false;
    end if;
 
    return true;

end function;


autom_H_H:= function( L, x2, y2, h2, C2, rep, hw )


    // two sets of can gens; first is the standard one...
    // return the automorphism mapping the first
    // to the second.
    // C2: Cartan matrix corr to second (ASSUME: simply laced).
    // rep: representation, hw: its highest weight, must make sure that that is
    // preserved.

    x1,y1,h1:= CanonicalGenerators(L);
    rank:= #x1;

    // look for hw vector wrt x2...
    n:= NumberOfRows( rep(L.1) );
    V:= VectorSpace( BaseField(L), n );
    K:= &meet[ sub< V | Basis(Kernel( Transpose(rep(u))) ) > : u in x2 ];

    v0:= Basis(K)[1]; 

    K:= &meet[ sub< V | Basis(Kernel( Transpose(rep(u))) ) > : u in x1 ];

    w0:= Basis(K)[1];

    hw:= [ Coordinates( K, w0*Transpose(rep(h1[i])) )[1] : i in [1..rank] ]; 

    _,_,_,C1:= RootSystem(L);

    // first we need to find an enumeration of C1 so
    // that Cartan matrices become equal.
    // For the moment we do this in the most inefficient way!
    C:= C1;
    Sn:= Sym(rank);
    pp:= [ u : u in Sn ];
    for p in pp do
        for i in [1..rank] do
            for j in [1..rank] do
                C[i][j]:= C1[i^p][j^p];
            end for;
        end for;

        if C eq C2 then
p;
//[v0*Transpose( rep(h2[i^(p^-1)]) ) - hw[i]*v0 : i in [1..rank]];

           if &and[ v0*Transpose( rep(h2[i^(p^-1)]) ) eq hw[i]*v0 : 
                i in [1..rank]] then 
              perm:= p;
              break;
           end if;
        end if;
    end for;
"perm found"; perm;
    x1:= [ x1[i^perm] : i in [1..rank] ];
    y1:= [ y1[i^perm] : i in [1..rank] ];
    h1:= [ h1[i^perm] : i in [1..rank] ];
 

    X1:= x1;
    X2:= x2;
    last1:= x1;
    last2:= x2;
    rts:= [ ];
    for i in [1..rank] do
        rt:= [ 0 : k in [1..rank] ];
        rt[i]:= 1;
        Append( ~rts, rt );
    end for;
    
    done:= false;
    while not done do

          newlast1:= [ ];
          newlast2:= [ ];
          newrts:= [ ];
          for i in [1..#last1] do
              for j in [1..#x1] do
                  if not IsZero(x1[j]*last1[i]) then
                     rt:= rts[i]; rt[j] +:=1;
                     if not rt in newrts then
                        Append( ~newrts, rt );
                        Append( ~newlast1, x1[j]*last1[i] );
                        Append( ~newlast2, x2[j]*last2[i] );
                     end if;
                  end if;
              end for;
          end for;
          if #newlast1 eq 0 then
             done:= true;
          else
             X1 cat:= newlast1;
             X2 cat:= newlast2;
             last1:= newlast1;
             last2:= newlast2;
             rts:= newrts;
          end if;

    end while; 

    Y1:= y1;
    Y2:= y2;
    last1:= y1;
    last2:= y2;
    rts:= [ ];
    for i in [1..rank] do
        rt:= [ 0 : k in [1..rank] ];
        rt[i]:= 1;
        Append( ~rts, rt );
    end for;
    
    done:= false;
    while not done do

          newlast1:= [ ];
          newlast2:= [ ];
          newrts:= [ ];
          for i in [1..#last1] do
              for j in [1..#y1] do
                  if not IsZero(y1[j]*last1[i]) then
                     rt:= rts[i]; rt[j] +:=1;
                     if not rt in newrts then
                        Append( ~newrts, rt );
                        Append( ~newlast1, y1[j]*last1[i] );
                        Append( ~newlast2, y2[j]*last2[i] );
                     end if;
                  end if;
              end for;
          end for;
          if #newlast1 eq 0 then
             done:= true;
          else
             Y1 cat:= newlast1;
             Y2 cat:= newlast2;
             last1:= newlast1;
             last2:= newlast2;
             rts:= newrts;
          end if;

    end while; 

    F:= BaseField(L);

    bas1:= h1 cat X1 cat Y1;
    bas1:= [ Vector( F, Eltseq(u) ) : u in bas1 ];    
    bas2:= h2 cat X2 cat Y2;

    sp:= VectorSpaceWithBasis( bas1 );

    return function( elm )

              c:= Coordinates( sp, Vector( F, Eltseq( elm )) );
              return &+[ c[i]*bas2[i] : i in [1..#c] ];
           end function;


end function;*/


posrtsNF:= function( posR, fund )


   V:= VectorSpaceWithBasis( fund );
   return [ Vector( Coordinates( V, V!u ) ) : u in posR ];


end function;

alphahat:= function( L, h, C, posR, fund )

   // the "hat" of the positive roots; h is the h-part of the
   // can gens, posR, C come out of rootsystem.


   rank:= NumberOfRows(C);
   eqns:= Matrix( Rationals(), 
      [ [ Trace( Matrix(AdjointMatrix( L, h[i] ))*
                 Matrix(AdjointMatrix(L,h[j])) ) : i in [1..rank] ] : 
                                                   j in [1..rank] ] );

   res:= [ ];
   C:= Matrix( Rationals(), C );
   pr:= posrtsNF( posR, fund );
   for i in [1..#pr] do
       sol:= Solution( eqns, Vector(Rationals(), pr[i])*C );
       Append( ~res, &+[ sol[i]*h[i] : i in [1..rank]]);
   end for;

   return res;

end function;

weylperms:= function( posR, fund, C )


    G:= Sym( 2*#posR );
    pr:= posrtsNF( posR, fund );
    R:= pr cat [ -u : u in pr ];
    gens:= [ ];
    rank:= #fund;
    for i in [1..#fund] do

        sq:= [ ];
        for j in [1..#R] do
            r:= R[j]- &+[ R[j][k]*C[k][i] : k in [1..rank] ]*pr[i];
            Append( ~sq, Index( R, r ) );
        end for;
        Append( ~gens, G!sq );
    end for;

    return G, sub< G | gens>, gens;

end function;


is_complex_reflection:= function( F, m )


   P<x>:= PolynomialRing( F );
   f:= CharacteristicPolynomial( m );
   f:= P!f;
   d:= Degree(f);

   rts:= Roots(f);

   if d eq 1 then
      return IsRootOfUnity( rts[1][2] );
   end if;

   if #rts ne 2 then return false; end if;
   i:= 1;
   if not rts[1][1] eq 1 then
      i:= 2;
   end if;
   if rts[i][1] ne 1 then return false; end if;
   if rts[i][2] ne d-1 then return false; end if;

   i:= -i+3;
   if not IsRootOfUnity( rts[i][2] ) then 
      return false;
   else
      return true;
   end if;
         

end function;

findgens0:= function( G, ref, k )

    // try to find k gens of G

    ind:= [Integers()!i : i in [1..k] ];
    n:= #ref;
    done:= false;
    while not done do

        g:= [ ref[i] : i in ind ];
        if sub< G | g > eq G then
           return g;
        else
           i:= k;
           while i ge 1 do
               if ind[i] eq n+i-k then 
                  i:= i-1; 
               else
                  break;
               end if;
           end while;

           if i gt 0 then
              ind[i]:= ind[i]+1;
              for j in [i+1..k] do
                  ind[j]:= ind[i]+j-i;
              end for;
           else
              done:= true;
           end if;
        end if;
    end while;

    return [];

end function;

findgens:= function( G, ref, m )

    a:= findgens0( G,ref, m );
    if #a eq 0 then 
       return findgens0( G, ref, m+1 );
    else
       return a;
    end if;

end function;

changering:= function( f, c, H, F )


    L:= Domain(f);
    K:= ChangeRing( L, F );
    imgs:= [ K!Vector( f(Basis(L)[i]) ) : i in [1..Dimension(L)] ];
    h:= hom< K -> K | imgs >;
    H0:= sub< K | [ K!Vector(L!Basis(H)[i]) : i in [1..Dimension(H)] ] >;
    c0:= [ K!Vector(c[i]) : i in [1..#c] ];

    return h, c0, H0;


end function;

/*cartanspdes:= procedure( c )

   for v in c do
       string:= "";
       s:= Eltseq(v);
       first:= true;
       for j in [1..#s] do
           if not IsZero(s[j]) then
              if s[j] eq 1 then
                 if first then
                    print "+"; first:= false;
                 end if;
                 print "v.",j;
              elif s[j] eq -1 then
                 print "-v.",j;
                 first:= false;
              else  
                 if not first then
                    print "+";
                 else
                    first:= false;
                 end if;
                 print "(",s[j],")*v.",j;
              end if;
           end if;
       end for;
    end for;

end procedure;*/

baschange_2comps:= function( p, C, x, y, h, C2, x2, y2, h2, F )

    // p an irrep of L wrt h, y, F base field
    // find a base change such that h2 is diagonal also
    // assume that h, h2 etc are wrt Cartan mat C, C2 
    // assume that the first basis vector, e1, is a hw vec.
  
    // returns a matrix M such that M*p(u)*M^-1 is diagonal for u in h2.

    // this is the version for a module with two direct summands...

    v1:= &meet[ Kernel( Transpose(p(u))) : u in x2 ];
    // so it is dim 2... split into two eigenvectors of h2...
    done:= false;
    while not done do 
       h0:= &+[ Random([-10..10])*u : u in h2 ];
       mat:= [ Coordinates( v1, Basis(v1)[1]*Transpose(p(h0)) ),Coordinates( v1, Basis(v1)[2]*Transpose(p(h0)) ) ];
       mat:= Matrix( F, mat );
       ev:= SetToSequence( Eigenvalues( mat ) );
       if #ev eq 2 then
          done:= true;
          w:= Basis( Eigenspace( mat, ev[1][1] ) )[1];
          w1:= w[1]*Basis(v1)[1]+w[2]*Basis(v1)[2];
          w:= Basis( Eigenspace( mat, ev[2][1] ) )[1];
          w2:= w[1]*Basis(v1)[1]+w[2]*Basis(v1)[2];
       end if;
    end while;

    v1:= &meet[ Kernel( Transpose(p(u))) : u in x ];
    // so it is dim 2... split into two eigenvectors of h...
    done:= false;
    while not done do 
       h0:= &+[ Random([-10..10])*u : u in h ];
       mat:= [ Coordinates( v1, Basis(v1)[1]*Transpose(p(h0)) ),Coordinates( v1, Basis(v1)[2]*Transpose(p(h0)) ) ];
       mat:= Matrix( F, mat );
       ev:= SetToSequence( Eigenvalues( mat ) );
       if #ev eq 2 then
          done:= true;
          w:= Basis( Eigenspace( mat, ev[1][1] ) )[1];
          u1:= w[1]*Basis(v1)[1]+w[2]*Basis(v1)[2];
          w:= Basis( Eigenspace( mat, ev[2][1] ) )[1];
          u2:= w[1]*Basis(v1)[1]+w[2]*Basis(v1)[2];
       end if;
    end while;

    n:= NumberOfRows( p(y[1]) );

    V1:= VectorSpaceWithBasis( [ u1 ] );
    hw1:= [ Coordinates( V1, u1*Transpose( p(u) ) )[1] : u in h ];

    V2:= VectorSpaceWithBasis( [ u2 ] );
    hw2:= [ Coordinates( V2, u2*Transpose( p(u) ) )[1] : u in h ];

    rank:= NumberOfRows(C);
    C0:= C2;
    Sn:= Sym(rank);
    pp:= [ u : u in Sn ];
    V1:= VectorSpaceWithBasis( [ w1 ] );
    V2:= VectorSpaceWithBasis( [ w2 ] );
    for q in pp do
        for i in [1..rank] do
            for j in [1..rank] do
                C0[i][j]:= C2[i^q][j^q];
            end for;
        end for;

        if C0 eq C then
           h1:= [ h2[i^q] : i in [1..#h2] ];
           hhw1:= [ Coordinates( V1, w1*Transpose( p(u) ) )[1] : u in h1 ];
           hhw2:= [ Coordinates( V2, w2*Transpose( p(u) ) )[1] : u in h1 ];
           if hhw1 eq hw1 and hhw2 eq hw2 then 
              pm:= q;
              break q;
           elif hhw1 eq hw2 and hhw2 eq hw1 then
              s:= w1;
              w1:= w2; w2:= s;
              pm:= q; 
              break q;
           end if;
        end if;
    end for;

    h1:= [ h2[i^pm] : i in [1..#h2] ];
    y1:= [ y2[i^pm] : i in [1..#y2] ];


    V:= KModule( F, n );
    bas11:= [ V!u1 ]; bas12:= [ V!u2 ];
    bas21:= [ V!w1 ]; bas22:= [ V!w2 ];
    oldlev11:= bas11; oldlev12:= bas12;
    oldlev21:= bas21; oldlev22:= bas22;
    done:= false; 
  
    Y1:= [ p(z) : z in y ];
    Y2:= [ p(z) : z in y1 ];

    while not done do
        newlev1:= [ ]; newlev2:= [ ];
        U:= sub< V | [ ] >;
        for i in [1..#oldlev11] do
            for j in [1..#y] do
                u:= oldlev11[i]*Transpose( Y1[j] );
                if not u in U then
                   Append( ~newlev1, u );
                   Append( ~newlev2, oldlev21[i]*Transpose( Y2[j] ) );
                   U:= sub< V | newlev1 >;
                end if;
            end for;
        end for;
        if #newlev1 eq 0 then 
           done:= true;
        else
           bas11 cat:= newlev1;
           bas21 cat:= newlev2;
           oldlev11:= newlev1;
           oldlev21:= newlev2; 
        end if;
    end while;

    done:= false;
    while not done do
        newlev1:= [ ]; newlev2:= [ ];
        U:= sub< V | [ ] >;
        for i in [1..#oldlev12] do
            for j in [1..#y] do
                u:= oldlev12[i]*Transpose( Y1[j] );
                if not u in U then
                   Append( ~newlev1, u );
                   Append( ~newlev2, oldlev22[i]*Transpose( Y2[j] ) );
                   U:= sub< V | newlev1 >;
                end if;
            end for;
        end for;
        if #newlev1 eq 0 then 
           done:= true;
        else
           bas12 cat:= newlev1;
           bas22 cat:= newlev2;
           oldlev12:= newlev1;
           oldlev22:= newlev2; 
        end if;
    end while;
    
    A:= Transpose( Matrix( F, [ Eltseq(u) : u in (bas11 cat bas12) ] ) );
    B:= Transpose( Matrix( F, [ Eltseq(u) : u in (bas21 cat bas22) ] ) );

    return A*B^-1;

end function;


baschange:= function( p, C, x, y, h, C2, x2, y2, h2, F )

    // p an irrep of L wrt h, y, F base field
    // find a base change such that h2 is diagonal also
    // assume that h, h2 etc are wrt Cartan mat C, C2 
    // assume that the first basis vector, e1, is a hw vec.
  
    // returns a matrix M such that M*p(u)*M^-1 is diagonal for u in h2.

    v1:= &meet[ Kernel( Transpose(p(u))) : u in x2 ];
    if Dimension(v1) eq 1 then
       v1:= Basis(v1)[1];
       dim1:= true;
    else
       return baschange_2comps( p, C, x, y, h, C2, x2, y2, h2, F );
    end if;

    n:= NumberOfRows( p(y[1]) );
    e1:= [ F!0 : i in [1..n] ];
    e1[1]:= F!1;
    e1:= Vector( F, e1 );

    V:= VectorSpaceWithBasis( [ e1 ] );
    hw:= [ Coordinates( V, e1*Transpose( p(u) ) )[1] : u in h ];

    rank:= NumberOfRows(C);
    C0:= C2;
    Sn:= Sym(rank);
    pp:= [ u : u in Sn ];
    V:= VectorSpaceWithBasis( [ v1 ] );
    for q in pp do
        for i in [1..rank] do
            for j in [1..rank] do
                C0[i][j]:= C2[i^q][j^q];
            end for;
        end for;

        if C0 eq C then
           h1:= [ h2[i^q] : i in [1..#h2] ];
           hw1:= [ Coordinates( V, v1*Transpose( p(u) ) )[1] : u in h1 ];
           if hw1 eq hw then 
              pm:= q;
              break q;
           end if;
        end if;
    end for;

    h1:= [ h2[i^pm] : i in [1..#h2] ];
    y1:= [ y2[i^pm] : i in [1..#y2] ];


    V:= KModule( F, n );
    e1:= V!Vector( e1 );
    v1:= V!v1;
    bas1:= [ e1 ];
    bas2:= [ v1 ];
    oldlev1:= [ e1 ];
    oldlev2:= [ v1 ];
    done:= false;
  
    Y1:= [ p(z) : z in y ];
    Y2:= [ p(z) : z in y1 ];

    while not done do
        newlev1:= [ ]; newlev2:= [ ];
        U:= sub< V | [ ] >;
        for i in [1..#oldlev1] do
            for j in [1..#y] do
                u:= oldlev1[i]*Transpose( Y1[j] );
                if not u in U then
                   Append( ~newlev1, u );
                   Append( ~newlev2, oldlev2[i]*Transpose( Y2[j] ) );
                   U:= sub< V | newlev1 >;
                end if;
            end for;
        end for;
        if #newlev1 eq 0 then 
           done:= true;
        else
           bas1 cat:= newlev1;
           bas2 cat:= newlev2;
           oldlev1:= newlev1;
           oldlev2:= newlev2; 
        end if;
    end while;
    
    A:= Transpose( Matrix( F, [ Eltseq(u) : u in bas1 ] ) );
    B:= Transpose( Matrix( F, [ Eltseq(u) : u in bas2 ] ) );

    return A*B^-1;

end function;


weylautom:= function( p, r )

    // p: permutation in W, r : output of previous fct, 
    // produce corr autom of L, this time as autom!

   
    rank:= # r`fund;
    N:= #r`posR;
    x:= [ ]; y:= [ ];

    for i in [1..rank] do
        j:= i^(p);
        if j le N then
           Append( ~x, r`posRv[j] );
           Append( ~y, r`negRv[j] );
        else
           Append( ~x, r`negRv[j-N] );
           Append( ~y, r`posRv[j-N] );
        end if;
    end for;

    for i in [1..rank] do
        V:= VectorSpaceWithBasis( Matrix( [ x[i] ]  ));
        c:= Coordinates( V, V!( (x[i]*y[i])*x[i] ) )[1];
        y[i]:= y[i]*2/c;
    end for;

    x0,y0,h0:= cangens( r`posRv, r`negRv, rank );

    pr:= posrtsNF( r`posR, r`fund );

    // compute the images: first the images of the pos rt vecs,
    // then those of the neg rt vecs, then those of h.

   imgsx:= [ ];
   posRv:= r`posRv;
   for i in [1..rank] do
       V:= VectorSpaceWithBasis( [ Vector( Eltseq( posRv[i] ) ) ] );
       c:= Coordinates( V, Vector( Eltseq( x0[i] ) ) )[1];
       Append( ~imgsx, x[i]/c );
   end for;

   imgsy:= [ ];
   negRv:= r`negRv;
   for i in [1..rank] do
       V:= VectorSpaceWithBasis( [ Vector( Eltseq( negRv[i] ) ) ] );
       c:= Coordinates( V, Vector( Eltseq( y0[i] ) ) )[1];
       Append( ~imgsy, y[i]/c );
   end for;

   for i in [rank+1..#pr] do
       
       ind:= 0;
       for k in [1..rank] do
           pos:= Index( pr, pr[i]-pr[k] );
           if pos gt 0 then
              ind:= k;
              break;
           end if;
       end for;

       V:= VectorSpaceWithBasis( [ Vector( Eltseq( posRv[i] ) ) ] );
       c:= Coordinates( V, Vector( Eltseq( posRv[ind]*posRv[pos] ) ) )[1];
   
       Append( ~imgsx, imgsx[ind]*imgsx[pos]/c );

       V:= VectorSpaceWithBasis( [ Vector( Eltseq( negRv[i] ) ) ] );
       c:= Coordinates( V, Vector( Eltseq( negRv[ind]*negRv[pos] ) ) )[1];
   
       Append( ~imgsy, imgsy[ind]*imgsy[pos]/c );

   end for;

   imbas:= imgsx cat imgsy cat [ x[i]*y[i] : i in [1..rank] ];
   bas:= posRv cat negRv cat h0;
   bas:= [ Vector( Eltseq(u) ) : u in bas ];    

   sp:= VectorSpaceWithBasis( bas );

   return function( elm )

              c:= Coordinates( sp, Vector( Eltseq( elm )) );
              return &+[ c[i]*imbas[i] : i in [1..#c] ];
           end function;


end function;


RF := recformat< theta, hom: Map, L,H : AlgLie, 
c,ref,posR,posRv,negRv,fund,gG,extraeqns,simref,liftrefs : SeqEnum,
permW : GrpPerm, C, U, thetalift : AlgMatElt, p : UserProgram >;


ThetaCen:= function( f, c, H, hw )

   // hw: highest weight of the rep

   L:= Domain(f);
   F:= BaseField(L);
   posR, posRv, negRv, fund, C:= rootsys( L, H );
   x,y,h:= cangens( posRv, negRv, #fund );
   a:= alphahat( L, h, C, posR, fund );      
   a:= a cat [ -u : u in a ];
   p:= [ Index( a, f(u) ) : u in a ];
   SN, W, simref:= weylperms( posR, fund, C );
   Z:= Centralizer( W, SN!p );

   rank:= #fund;
   bas:= a[[1..rank]];
   sp:= VectorSpaceWithBasis( [ Vector( Eltseq(u) ) : u in bas ] );

   gens:= [ Matrix( F, [ Coordinates( sp, Vector(Eltseq(a[i^(Z.j)])) ) : 
               i in [1..rank] ]) : j in [1..Ngens(Z)] ];

   //now wrt c!

   spc:= VectorSpaceWithBasis( [ Vector( Eltseq(u) ) : u in c ] );

   gensc:= [ ];
   for x in gens do
       
       m:= [ ];
       for v in c do
           cf:= Vector( F, Coordinates( sp, Vector(Eltseq( v )) ))*x;
           u:= &+[ cf[i]*bas[i] : i in [1..rank] ];
           Append( ~m, Coordinates( spc, Vector( u ) ) );
       end for;
       Append( ~gensc, Matrix(m) );
   end for;

//return Z, bas, gens, gensc;


   gl:= GL( #c, F );
   G:= sub< gl | gensc >;

   ref:= [ ];
   elms:= [ u : u in G ];
   for u in elms do
       if is_complex_reflection(F,u) then 
          Append( ~ref, u );
       end if;
   end for;

   gG:= findgens( G, ref, #c );

//return gG;

//"gens done; start irrep";

   p:= irrep( L, SemisimpleType(L)[1], rank, hw, C, h, x, y );

//"irrep done";
   
//tt0:= Cputime();

   R, b, extraeqns:= weylcenelt( L, x, y, h );
   v:= liftaut_diag( L, R, p, b, x, y );
   v:= v[1];
   dens:= [ ];
   for i in [1..NumberOfColumns(v)] do
       for j in [1..NumberOfRows(v)] do
           Append( ~dens, Denominator( v[i][j] ) );
       end for;
   end for;
   v:= v*Lcm(dens);
//Cputime()-tt0;
//"Weyl cen elt done";
//tt0:= Cputime();
   x2,y2,h2:= CanonicalGenerators(L);
   C2:= CartanMatrix( RootDatum(L) );
   M:= baschange( p, C, x, y, h, C2, x2, y2, h2, F );
   p0:= function(x) return M*p(x)*M^-1; end function;
   th:= liftaut_diag( L, F, p0, f, x2, y2 ); th:= th[1];
   th:= M^-1*th*M;
//Cputime()-tt0;

   r:= rec< RF | theta:= f, thetalift:= th,
     hom:= hom< Z -> G | [ <Z.i,G!gensc[i]> : i in [1..#gensc] ] >,
     L:= L, H:= H,
     c:= c, ref:= ref, posR:= posR, posRv:= posRv, negRv:= negRv,
     fund:= fund, gG:= gG,permW:= W, C:=C, 
     U:= v, simref:= simref,
     p:=p, extraeqns:= extraeqns  >;
//tt0:= Cputime();     
   liftrefs:= [ ];
   for pm in simref do
      w0:= weylautom( pm, r );
      mat:= Matrix( [ Coordinates( L, w0(L.i) ) : i in [1..Dimension(L)]] );
      b0:= Basis( Kernel( mat-IdentityMatrix(F,Dimension(L)) ) );
      K:= sub< L | [ L!u : u in b0 ] >;
      HK:= CartanSubalgebra( Centralizer( K, K meet H ) );

      pR, pRv, nRv, fnd, Cpm:= rootsys( L, HK );

      if #fnd lt #simref then // non splitting CSA... old method
//"grrr non split CSA";

          facs:= [  ];
          for x in Basis(HK) do
              ad:= AdjointMatrix( L, L!x );
              facs cat:= [ u[1] : u in 
                    Factorization(CharacteristicPolynomial(ad)) ];
          end for;
          f0:= [ ];

          for fc in facs do
              if not fc in f0 then Append( ~f0, fc ); end if;
          end for;
          facs:= f0;

          F0:= BaseField(L);
          P:= PolynomialRing( F0 );
          for fc in facs do
              E:= SplittingField( P!fc );
              if E ne F0 then
                 F0:= E; P:= PolynomialRing(F0);
              end if;
          end for;

          L0:= ChangeRing( L, F0 );
          p1:= function(x)
               cf:= Eltseq(x);
               u:= Matrix( E, p(Zero(L)) );
               for i in [1..#cf] do
                   v:= Matrix( E, p(L.i) );
                   u +:= cf[i]*v;
               end for;
               return u;
          end function;

          w1:= function(x)
               cf:= Eltseq(x);
               u:= Zero(L0);
               for i in [1..#cf] do
                   v:= L0!Eltseq( w0(L.i) );
                   u +:= cf[i]*v;
               end for;
               return u;
          end function;

          HL0:= sub< L0 | [ L0!Eltseq( L!u ) : u in Basis(HK) ] >;
          pR, pRv, nRv, fnd, Cpm:= rootsys( L0, HL0 );
          xp,yp,hp:= cangens( pRv, nRv, #fnd );
          M:= baschange( p1, C, [ L0!Eltseq(u) : u in x], [ L0!Eltseq(u) : u in y], [ L0!Eltseq(u) : u in h], 
                                                       Cpm, xp, yp, hp, F0 );
          p0:= function(x) return M*p1(x)*M^-1; end function;
          ww:= liftaut_diag( L, E, p0, w1, xp, yp ); ww:= ww[1];
          ww:= M^-1*ww*M;

          // PROBLEM: ww lies in GL(F,m), up to a scalar factor in F0...
          // look for an element in ww not in F, divide by it, and hope that gets the
          // whole matrix over F...
          for i in [1..NumberOfRows(ww)] do 
              for j in [1..NumberOfRows(ww)] do 
                  if not ww[i][j] in F then
                     ww:= ww*( ww[i][j]^-1 );
                     break i;
                  end if;
              end for;
           end for;

          ww:= Matrix( F, ww );
          Append( ~liftrefs, ww );

//F0;
//BaseField(L);
//print ww;
//print liftaut( L, F, p, w0, x, y );
         //ww:= liftaut( L, F, p, w0, x, y );
         //Append( ~liftrefs, ww[1] );
      else

         xp,yp,hp:= cangens( pRv, nRv, #fnd );
         M:= baschange( p, C, x, y, h, Cpm, xp, yp, hp, F );
         p0:= function(x) return M*p(x)*M^-1; end function;
         ww:= liftaut_diag( L, F, p0, w0, xp, yp ); ww:= ww[1];
         ww:= M^-1*ww*M;
         Append( ~liftrefs, ww );
      end if;
   end for;
//print "liftrefs time.... ",Cputime()-tt0;
   r`liftrefs:= liftrefs;
   return r;

end function;


ThetaCenSC:= function( f, c, H )

   L:= Domain(f);
   F:= BaseField(L);
   posR, posRv, negRv, fund, C:= rootsys( L, H );
   x,y,h:= cangens( posRv, negRv, #fund );
   a:= alphahat( L, h, C, posR, fund );      
   a:= a cat [ -u : u in a ];
   p:= [ Index( a, f(u) ) : u in a ];

   SN, W:= weylperms( posR, fund, C );

//   return W,SN!p, C;

   Z:= Centralizer( W, SN!p );

   rank:= #fund;
   bas:= a[[1..rank]];
   sp:= VectorSpaceWithBasis( [ Vector( Eltseq(u) ) : u in bas ] );

   gens:= [ Matrix( F, [ Coordinates( sp, Vector(Eltseq(a[i^(Z.j)])) ) : 
               i in [1..rank] ]) : j in [1..Ngens(Z)] ];

   //now wrt c!

   spc:= VectorSpaceWithBasis( [ Vector( Eltseq(u) ) : u in c ] );

   gensc:= [ ];
   for x in gens do
       
       m:= [ ];
       for v in c do
           cf:= Vector( F, Coordinates( sp, Vector(Eltseq( v )) ))*x;
           u:= &+[ cf[i]*bas[i] : i in [1..rank] ];
           Append( ~m, Coordinates( spc, Vector( u ) ) );
       end for;
       Append( ~gensc, Matrix(m) );
   end for;

//   return Z, bas, gens, gensc;

   gl:= GL( #c, F );
   G:= sub< gl | gensc >;

   ref:= [ ];
   elms:= [ u : u in G ];
   for u in elms do
       if is_complex_reflection(F,u) then 
          Append( ~ref, u );
       end if;
   end for;

   gG:= findgens( G, ref, #c );

//return G,gG,ref;

   return rec< RF | theta:= f, 
     hom:= hom< Z -> G | [ <Z.i,G!gensc[i]> : i in [1..#gensc] ] >,
     L:= L, H:= H, 
     c:= c, ref:= ref, posR:= posR, posRv:= posRv, negRv:= negRv,
     fund:= fund, gG:= gG,permW:= W  >;
     
   

end function;



weylautomSC:= function( p, r )

    // p: permutation in W, r : output of previous fct, 
    // produce corr autom of L

   
    rank:= # r`fund;
    N:= #r`posR;
    x:= [ ]; y:= [ ];

    for i in [1..rank] do
        j:= i^(p);
        if j le N then
           Append( ~x, r`posRv[j] );
           Append( ~y, r`negRv[j] );
        else
           Append( ~x, r`negRv[j-N] );
           Append( ~y, r`posRv[j-N] );
        end if;
    end for;

    for i in [1..rank] do
        V:= VectorSpaceWithBasis( Matrix( [ x[i] ]  ));
        c:= Coordinates( V, V!( (x[i]*y[i])*x[i] ) )[1];
        y[i]:= y[i]*2/c;
    end for;

    x0,y0,h0:= cangens( r`posRv, r`negRv, rank );

    pr:= posrtsNF( r`posR, r`fund );

    // compute the images: first the images of the pos rt vecs,
    // then those of the neg rt vecs, then those of h.

   imgsx:= [ ];
   posRv:= r`posRv;
   for i in [1..rank] do
       V:= VectorSpaceWithBasis( [ Vector( Eltseq( posRv[i] ) ) ] );
       c:= Coordinates( V, Vector( Eltseq( x0[i] ) ) )[1];
       Append( ~imgsx, x[i]/c );
   end for;

   imgsy:= [ ];
   negRv:= r`negRv;
   for i in [1..rank] do
       V:= VectorSpaceWithBasis( [ Vector( Eltseq( negRv[i] ) ) ] );
       c:= Coordinates( V, Vector( Eltseq( y0[i] ) ) )[1];
       Append( ~imgsy, y[i]/c );
   end for;

   for i in [rank+1..#pr] do
       
       ind:= 0;
       for k in [1..rank] do
           pos:= Index( pr, pr[i]-pr[k] );
           if pos gt 0 then
              ind:= k;
              break;
           end if;
       end for;

       V:= VectorSpaceWithBasis( [ Vector( Eltseq( posRv[i] ) ) ] );
       c:= Coordinates( V, Vector( Eltseq( posRv[ind]*posRv[pos] ) ) )[1];
   
       Append( ~imgsx, imgsx[ind]*imgsx[pos]/c );

       V:= VectorSpaceWithBasis( [ Vector( Eltseq( negRv[i] ) ) ] );
       c:= Coordinates( V, Vector( Eltseq( negRv[ind]*negRv[pos] ) ) )[1];
   
       Append( ~imgsy, imgsy[ind]*imgsy[pos]/c );

   end for;

   return imgsx cat imgsy cat [ x[i]*y[i] : i in [1..rank] ], h0;

end function;



poleqnsSC:= function( p, r )

   // p: permutation of Weyl group


   L:= r`L;
   theta:= r`theta;
   rank:= #r`fund;

   F:= BaseField(L);
   names:= [ ];
   for i in [1..rank] do
       Append( ~names, "x" cat IntegerToString(i) );
   end for;

   for i in [1..rank] do
       Append( ~names, "y" cat IntegerToString(i) );
   end for;

   Q:= PolynomialRing( F, 2*rank, "grevlex" );
   AssignNames( ~Q, names );

   x:= [ Q.i : i in [1..rank] ];
   y:= [ Q.i : i in [rank+1..2*rank] ];

   bas:= r`posRv;
   n:= #bas;
   bas:= bas cat r`negRv;

   imgs, h0:= weylautomSC( p, r );

   bas cat:= h0;
   B:= VectorSpaceWithBasis( Matrix( bas ) ); 

   pr:=  posrtsNF( r`posR, r`fund );
   pr:= [ Vector( Integers(), u ) : u in pr ];
   eqns:= [ ];

   posn:= [1..rank] cat [n+1..n+rank];
   for pos in posn do

       u:= Coordinates( B, Vector(imgs[pos]) );
       if pos le rank then
          u:= x[pos]*Vector( Q, u );
       else
          u:= y[pos-n]*Vector( Q, u );
       end if;

       v:= Zero(Q)*u;
       for j in [1..Dimension(L)] do
           if not IsZero(u[j]) then
              v:= v + u[j]*Vector( Q, Coordinates(B,Vector(theta(bas[j]))) );
           end if;
       end for;

       u:= Coordinates( B, Vector(theta(bas[pos])) );
       u:= Vector( Q, u );
       for j in [1..Dimension(L)] do
           if j le n then
              u[j]:= u[j]*(&*[ x[k]^pr[j][k] : k in [1..rank] ] );
           elif j le 2*n then
              u[j]:= u[j]*(&*[ y[k]^pr[j-n][k] : k in [1..rank] ] );
           end if;
       end for;     

       w:= Zero(Q)*u;
       for j in [1..Dimension(L)] do
           if not IsZero(u[j]) then
              w:= w + u[j]*Vector( Q, Coordinates(B,Vector(imgs[j])) );
           end if;
       end for;

       eqns cat:= Eltseq( v-w );

   end for;       

   eqns cat:= [ x[i]*y[i]-One(Q) : i in [1..rank] ];
   eqns:= SequenceToSet( eqns );
   eqns:= SetToSequence( eqns );

   return eqns;

end function;


/*check:= function( p, r )

   bas:= r`posRv;
   n:= #bas;
   bas:= bas cat r`negRv;

   imgs, h0:= weylautom( p, r );

   bas cat:= h0;
   B:= VectorSpaceWithBasis( Matrix( bas ) ); 
   N:= #bas;
   for i in [1..N] do
       for j in [i+1..N] do

           cf:= Coordinates( B, Vector( bas[i]*bas[j] ) );
           u:= &+[ cf[k]*imgs[k] : k in [1..N] ];
           if u ne imgs[i]*imgs[j] then
u-imgs[i]*imgs[j];
              return [i, j]; 
           end if;
       end for;
   end for;

   return true;


end function; */

/*
> pms:= r`simref;
> s:= InverseWordMap(W);
> pm:= Random( W );
> FG:= FreeGroup(6);
> fct:= pm@s;
> a:= Eltseq( Evaluate( fct, FG ) );
> a;
[ 3, 2, 4, 5, 1, 2, 3, 4, 2, 1, 5, 4, 2, 3, 6, 3, 2, 1, 4, 2, 3 ]
> pm1:= &*[ pms[i] : i in a ];
> pm1 eq pm;
true
*/


poleqns:= function( perm, r )

   // perm:= permutation in the Weyl group

   L:= r`L;
   F:= BaseField(L);

//tt0:= Cputime();

/*   w0:= weylautom( perm, r );
   w1:= liftaut( L, F, r`p, w0, r`x, r`y );
   w1:=w1[1];       */

   lfts:= r`liftrefs;
   FG:= FreeGroup( #lfts );
   s:= InverseWordMap(r`permW);
   fct:= perm@s;
   a:= Eltseq( Evaluate( fct, FG ) );
   for i in [1..#a] do if a[i] lt 0 then a[i]:= -a[i]; end if; end for;
   w1:= &*[ lfts[i] : i in a ];

   rank:= #r`fund;

/*   names:= [ ];
   for i in [1..rank] do
       Append( ~names, "x" cat IntegerToString(i) );
   end for;
   P:= PolynomialRing( F, rank, "grevlex" );
   AssignNames( ~P, names ); */

   U:= r`U;
   th:= r`thetalift;
   //th:= liftaut( L, F, r`p, r`theta ); th:= th[1];
//"liftaut done for theta...";
//   Q:=w1*U-U*w1;

   Q:= th*w1*U -w1*U*th;

   eqns:= Eltseq(Q);           
   eqns:= SequenceToSet( eqns);
   eqns:= SetToSequence(eqns); 

   eqns:= [ Numerator(u) : u in eqns ] cat r`extraeqns;
   P:= ChangeOrder( Parent( eqns[1] ), "grevlex" );
   return [ P!p : p in eqns ];

end function;

OP := recformat< 
theta: Map, 
L : AlgLie, 
c,grading, generators : SeqEnum,
F : Fld,
Wc : Grp >;


LittleW:= function( type, rank, f, order, scon )

    // size: the size of the little Weyl group, which we assume known...

   if scon then
      sc:= true;
   elif type eq "E" then 
      if rank eq 6 then
         hw:= [[1,0,0,0,0,0]]; sc:= false;
      elif rank eq 7 then
         hw:= [[0,0,0,0,0,0,1]]; sc:= false;
      elif rank eq 8 then
         sc:= true;
      end if;
   elif type eq "D" and rank mod 2 eq 1 then
      hw:= [ 0 : i in [1..rank] ];
      hw[rank]:= 1;
      hw:= [ hw ];
      sc:= false;
   elif type eq "D" and rank mod 2 eq 0 then
      hw:= [ 0 : i in [1..rank] ];
      hw[1]:= 1;
      hw:= [ hw ];
      hw1:= [ 0 : i in [1..rank] ];
      hw1[rank]:= 1;
      Append( ~hw, hw1 );
      sc:= false;
   elif type eq "A" then
      hw:= [ 0 : i in [1..rank] ];
      hw[1]:= 1;
      hw:= [ hw ];
      sc:= false;
   elif type  eq "B" then
      hw:= [ 0 : i in [1..rank] ];
      hw[rank]:= 1;
      hw:= [ hw ];
      sc:= false;
   elif type eq "C" then
      hw:= [ 0 : i in [1..rank] ];
      hw[1]:= 1;
      hw:= [ hw ];
      sc:= false;
   end if;

   //f,order:= FiniteOrderInnerAutomorphism( type, diag );
   L:= Domain(f);
   F:= BaseField(L);
   w:= F.1;
   gr:= GradingThetaRep( f, order, w );
   cc,HH,FF:= find_cart( L, gr[1], gr[2], 1, Infinity() );

   if #cc[1] eq 0 then return [], [], [], [], []; end if;

   k:= 1;
   deg:= Degree( FF[1] );
   while deg gt Degree(F) and k le 9 do
      k:= k+1;
      cc0,HH0,FF0:= find_cart( L, gr[1], gr[2], 1, deg );
      if #FF0 gt 0 then
         if Degree(FF0[1]) lt deg then
            cc:= cc0; HH:= HH0; FF:= FF0; deg:= Degree(FF[1]);
         end if;
      end if;
   end while;

   if deg gt Degree(F) then
      f, c, H:= changering( f, cc[1], HH[1], FF[1] );
      F:= FF[1];
   else
      c:= cc[1]; H:= HH[1];
   end if;
   if sc then
      r:= ThetaCenSC( f, c, H );
   else
      r:= ThetaCen( f, c, H, hw );
   end if;

   g:= r`gG;
   //gl:= GL( #c, F );
   //G0:= sub< gl | [ gl!a : a in g ] >;

   //if #G0 eq size then return g, c, F; end if;

   hom:= r`hom;
   K:= [ u : u in Kernel( hom )];
//print "size kernel... ",#K;

//   GBc:= 0;
//   GBt:= 0;

   gens:= [ ];
   for u in g do
       k:= 1;
       found:= false;
       pre:= u@@hom;
       while k le #K and not found do
          p:= pre*K[k]; 
          if sc then
             qq:= poleqnsSC( p, r );
          else
             qq:= poleqns( p, r );
          end if;
//print "start GB... ",k;
//t:= Cputime();
          G:= GroebnerBasis(qq);
//GBt +:= Cputime()-t;
//GBc +:= 1;
//print "end GB...";
          if G[1] ne 1 then
             found:= true;
             Append( ~gens, u );
          end if;
          k:= k+1;
       end while;
   end for;

   if #gens lt #g then
      gl:= GL( #c, F );
      G0:= sub< gl | [ gl!a : a in g ] >;
      ref:= r`ref;
      ref0:= [ ];
      for u in ref do 
          if not u in g then
             Append( ~ref0, u );
          end if;
      end for;
//print "length ref0...",#ref0;      
      for u in ref0 do
          k:= 1;
          found:= false;
          pre:= u@@hom;
          while k le #K and not found do
             p:= pre*K[k];
             if sc then
                qq:= poleqnsSC( p, r );
             else    
                qq:= poleqns( p, r );
             end if;
//print "start GB... ",k;
//t:= Cputime();
             G:= GroebnerBasis(qq);
//GBt +:= Cputime()-t;
//GBc +:= 1;
//print "end GB...";
             if G[1] ne 1 then
                found:= true;
                Append( ~gens, u );
                S0:= sub< gl | [ gl!a : a in gens ] >;
                p:= Integers()!( #G0/#S0 );
                if IsPrime(p) then break u; end if;
             end if;
             k:= k+1;
          end while;
      end for;
   else
//print "no GB, GB time", GBc, GBt;
      return gens, c, F, f, gr;
   end if;    

   gl:= GL( #c, F );
   G:= sub< gl | gens >;

   ref:= [ ];
   elms:= [ u : u in G ];
   for u in elms do
       if is_complex_reflection(F,u) then 
          Append( ~ref, u );
       end if;
   end for;

   gG:= findgens( G, ref, #c );

//print "size little Weyl group... ",#G;
//print "no GB, GB time", GBc, GBt;
   return gG, c, F, f, gr;

end function;



intrinsic LittleWeylGroup(type::MonStgElt, rank::RngIntElt , diag::SeqEnum ) -> Rec
{The little Weyl group of the grading defined by the Kac diagram diag.}

   // NEED #diag = rank+1

   require #diag eq rank+1: "Diagram does not define an inner automorphism.";

   f,order:= FiniteOrderInnerAutomorphism( type, diag );
   g,c,F,f,gr:= LittleW( type, rank, f, order, false );

   require #g gt 0: "The rank is zero.";

   gl:= GL( #c, F );
   r:= rec< OP | theta:= f, c:= c, generators:= g, F:= F, Wc:= sub< gl | g >, L:= Domain(f), grading:= gr >;

   return r;

end intrinsic;

intrinsic ProbableLittleWeylGroup(type::MonStgElt, rank::RngIntElt , diag::SeqEnum ) -> Rec
{A group containing the little Weyl group of the grading defined by the Kac diagram diag, and in many cases equal to it.}

   if #diag eq rank+1 then
      f,order:= FiniteOrderInnerAutomorphism( type, diag );
   else
      if type eq "A" then
         if rank mod 2 eq 0 then
            lab:= [ 2 : i in [1..#diag] ];
            lab[ 1 ]:= 1;
            d:= 2;
         else
            lab:= [ 2 : i in [1..#diag] ];
            lab[ #lab ]:= 1;
            lab[1]:= 1; lab[2]:= 1;
            d:= 2;
         end if;
      elif type eq "D" then
         if rank eq 4 and #diag eq 3 then
            lab:= [1,2,1];
            d:= 3;
         else
            lab:= [ 1 : i in [1..#diag] ];
            d:= 2;
         end if;
      elif type eq "E" and rank eq 6 then
         //lab:= [1,2,3,2,1];
         lab:= [1,1,2,3,2];
         d0:= [ diag[i] : i in [1,5,4,3,2] ];
         diag:= d0;
         d:= 2;
      end if;
      order:= &+[ lab[i]*diag[i] : i in [1..#diag] ];
      order:= d*order;
      f:= FiniteOrderOuterAutomorphism( type, rank, d, order, diag );
   end if;

   g,c,F, f, gr:= LittleW( type, rank, f, order, true );

   require #g gt 0: "The rank is zero.";

   gl:= GL( #c, F );
   r:= rec< OP | theta:= f, c:= c, generators:= g, F:= F, Wc:= sub< gl | g >, L:= Domain(f), grading:= gr >;

   return r;

end intrinsic;

