/*
 This file has implementations of two functions:

* DecompositionOfAlgebraicLieAlgebra( L )

   Here L is a sequence of matrices, over a field of characteristic 0, forming a basis of an algebraic
   Lie algebra. This function returns three sequences: S, T, N. S spans a Levi subalgebra, T a torus, 
   and N the nilpotent radical (maximal ideal consisting of nilpotent matrices). All three subalgebras
   are algebraic. T+N is the solvable radical of L.

* GeneratorsOfArithmeticGroup( L )

  Here L is a sequence of matrices, over Q, spanning a solvable algebraic Lie algebra. Let G be the
  corresponding connected algebraic group. This function returns a list of generators of a finite index
  subgroup of G(Z). If L is a torus, or if L consists of nilpotent matrices, then these actually generate
  G(Z). 

Example:

> Attach("arithgrp.m"); 
> P<x>:= PolynomialRing( Rationals() );
> A:= CompanionMatrix( x^4+2*x^3-x+1 );
> L:= [ A^0, A, A^2+2/3*A^3];
> GeneratorsOfArithmeticGroup( L );
Dirichlet time.... 0.000
Full torus time... 0.000
Splitting field time... 0.010
dimension splitting field... 8
Ge time... 0.060
Ge time... 0.240
Group structure... [ 2, 0 ]
[
    [-1  0  0  0]
    [ 0 -1  0  0]
    [ 0  0 -1  0]
    [ 0  0  0 -1],

    [ 0  0  1  1]
    [-1  1  0 -1]
    [ 1 -2  1  2]
    [-2  3 -2 -3]
]


Here L is an algebraic diagonalisable Lie algebra. Currently the function prints some
information about what it is doing, that may ony be understandable to the author of 
the code.

*/

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//  Decomposition of an algebraic Lie algebra: decomposition(L);
//
//

fitting_1:= function( L, K )

   // K is a nilpotent Lie subalgebra of L
   // we compute the Fitting one component L_1(\ad K)
   // both L and K are given as matrix Lie algebras

   n:= NumberOfRows( Basis(L)[1] );
   M:= RMatrixSpace( Rationals(), n, n );
   comp:= sub< M | [ M!u : u in Basis(L) ] >;
   done:= false;
   while not done do
      b:= Basis( comp );
      e:= [ ];
      for i in [1..Dimension(K)] do
          x:= M!K.i;
          for y in b do
              Append( ~e, x*y-y*x );
          end for;
      end for;
      comp0:= sub< M | [ M!u : u in e ] >;
      if Dimension(comp0) eq Dimension(comp) then
         done:= true;
      else
         comp:= comp0;
      end if;
   end while;

   return [ L!u : u in Basis( comp )];

end function;

intrinsic DecompositionOfAlgebraicLieAlgebra( L::SeqEnum ) -> SeqEnum, SeqEnum, SeqEnum
{L is a list of basis elements of an algebraic Lie algebra in char 0; this returns three lists: Levi subalgebra, torus and nilpotent radical}

   // here L is an algebraic matrix Lie algebra
   // this returns three lists of matrices: the first
   // spans a Levi subalgebra, the second a maximal torus
   // in the solvable radical, the third spans the
   // nilpotent radical. L is the direct sum of these.

   n:= NumberOfRows( L[1] );

   V:= MatrixLieAlgebra( Rationals(), n );
   LL:= sub< V | [ V!u : u in L ] >;
   _,K:= HasLeviSubalgebra( LL );

   S:= SolvableRadical( LL );
   bS:= Basis(S);
   t:= [ ];
   for i in [1..Dimension(S)] do
       for j in [1..Dimension(S)] do
           cf:= Coordinates( S, bS[i]*bS[j] );
           for k in [1..#cf] do
               if cf[k] ne 0 then
                  Append( ~t, <i,j,k,cf[k]> );
               end if;
           end for;
       end for;
   end for;
   ss:= LieAlgebra< Rationals(), Dimension(S) | t >;
   hh:= CartanSubalgebra( ss );
   H:= [ ];
   for u in Basis(hh) do
       cf:= Eltseq( ss!u );
       Append( ~H, &+[ cf[k]*bS[k] : k in [1..#cf] ] );
   end for;
   H:= sub< S | [ S!u : u in H] >;

   //H:= CartanSubalgebra( S );

   M:= RMatrixSpace( Rationals(), n, n );
   N:= fitting_1( S, H );

   T:= [ ];
   for u in Basis(H) do
       sem,nil:= JordanDecomposition( Matrix(u) );
       Append( ~T, sem ); 
       Append( ~N, nil );
   end for;

   T:= sub< M | [ M!u : u in T ] >;
   T:= Basis(T);
   N:= sub< M | [ M!u : u in N ] >;
   N:= Basis(N);
   K:= sub< M | [ M!u : u in Basis(K) ] >;
   K:= Basis(K);

   return K,T,N;

end intrinsic;

////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//  Extension function from the grp alg files...
//
mult:= function( T, m, a, b )

    c:= [ ];
    for k in [1..#m] do
        s:= 0;
        for i in [1..#m] do 
            for j in [1..#m] do 
                s+:= T[i][j][k]*a[i]*b[j];
            end for;
        end for;
        c[k]:= s mod m[k];
    end for;

    return c;

end function;

enum:= function( A, m, one, gen ) // enumerate the group...

    elms:= [ one ];
    for g in gen do
        if not g in elms then

           old:= elms;
           while #old gt 0 do
              new:= [ ];
              for u in old do
                  a:= Eltseq( (A!u)*(A!g) ); for k in [1..#m] do a[k]:= a[k] mod m[k]; end for;
                  if not a in elms then
                     Append( ~new, a );
                  end if;
              end for;
              elms cat:= new;
              old:= new;
           end while;
        end if;
    end for;

    return elms;

end function;


enum_withexps:= function( A, m, one, ord, gen ) // enumerate the group...

    elms:= [ one ];
    exps:= [ [ ] ];
    for j in [1..#gen] do
   
        g:= gen[j];
        for i in [1..#exps] do exps[i] cat:= [0]; end for;
        new:= [ ];
        newexps:= [ ];
        for k in [1..#elms] do
            a:= elms[k];
            exp:= exps[k];
            for i in [1..ord[j]-1] do
                a:= Eltseq( (A!a)*(A!g) );
                for l in [1..#m] do a[l]:= a[l] mod m[l]; end for;
                Append( ~new, a );
                exp0:= exp;
                exp0[j]:= i;
                Append( ~newexps, exp0 );
            end for;
        end for;

        elms cat:= new;
        exps cat:= newexps;

    end for; 

    ParallelSort( ~elms, ~exps );
    return elms, exps;

end function;

pow:= function( A, m, one, g, n ) // compute g^n by repeated squaring... 

   if n eq 0 then return one; end if;
   b:= [ ];
   k:= n;
   while k ne 0 do
       s:= k mod 2;
       Append( ~b, s );
       k:= Integers()!((k-s)/2);
   end while;

   r:= #b-1;
   s:= [ g ];
   a:= g;
   for i in [1..r] do
       a:= Eltseq( (A!a)^2 );
       for k in [1..#m] do a[k]:= a[k] mod m[k]; end for;
       Append( ~s, a );
   end for;

   c:= one;
   for i in [1..r+1] do
       if b[i] eq 1 then
          c:= Eltseq( (A!c)*(A!s[i]) );
          for k in [1..#m] do c[k]:= c[k] mod m[k]; end for;
       end if;
   end for;

   return c;

end function;

enum_withexps0:= function( A, m, one, gen ) // enuerate the group...

    elms:= [ one ];
    exps:= [ [ ] ];
    for g in gen do
        if not g in elms then
           
           for i in [1..#exps] do exps[i] cat:= [0]; end for;
           old:= elms;
           oldex:= exps;
           while #old gt 0 do
              new:= [ ];
              newex:= [ ];
              for l in [1..#old] do

                  a:= Eltseq( (A!old[l])*(A!g) );
                  for k in [1..#m] do a[k]:= a[k] mod m[k]; end for;

                  if not a in elms then
                     Append( ~new, a );
                     ex:= oldex[l]; ex[#ex]+:= 1;
                     Append( ~newex, ex ); 
                  end if;
              end for;
              elms cat:= new;
              exps cat:= newex;
              ParallelSort( ~elms, ~exps );
              old:= new;
              oldex:= newex;
           end while;
        else
           for i in [1..#exps] do exps[i] cat:= [0]; end for;
        end if;
    end for;

    return elms, exps;

end function;



gens:= function( A, m, one, G )

   // G is a group given as set of elts, we find a generating set

   elms:= [ one ];
   gens:= [ ];

   while  #elms lt #G do
       // find an element outside elms...

       g:= G[1]; k:= 1;
       while g in elms do 
            k:= k+1; g:= G[k];
       end while;

       Append( ~gens, g );
       old:= elms;
       while #old gt 0 do
           new:= [ ];
           for u in old do
               a:= Eltseq( (A!u)*(A!g) ); for k in [1..#m] do a[k]:= a[k] mod m[k]; end for;
               if not a in elms then
                  Append( ~new, a );
               end if;
           end for;
           elms cat:= new;
           old:= new;
       end while;
   end while;
       
   return gens;

end function;

normalgens:= function( A, m, one, gens )

   G:= [ one ];
   gg:= [ ];
   exps:= [ [ ] ];
   orders:= [ ];
   for k in [1..#gens] do
       if not gens[k] in G then
          // determine its smallest power in G...
          h:= gens[k];
          f:= 1;
          while not h in G do
              h:= Eltseq( (A!h)*(A!gens[k]) ); for i in [1..#m] do h[i]:= h[i] mod m[i]; end for;
              f+:= 1;
          end while;
          p:= Index( G, h );
          len:= #gg+1;
          mat:= [ ];
          for i in [1..#gg] do
              v:= [ 0 : j in [1..len] ];
              v[i]:= orders[i];
              Append( ~mat, v );
          end for;
          Append ( ~mat, exps[p] cat [-f] );
          mat:= Matrix( Integers(), mat );
          S,P,Q:= SmithForm( mat );

          newgg:= [ ];
          Qi:= Q^-1;

          s:= f;
          while h ne one do
               h:= Eltseq( (A!h)*(A!gens[k]) ); for o in [1..#m] do h[o]:= h[o] mod m[o]; end for;
               s +:= 1;
          end while;

          Append( ~orders, s );
          Append( ~gg, gens[k] );

          neword:= [ ];         
          for i in [1..NumberOfColumns(S)] do
              if S[i][i] ne 1 then
                 a:= one;
                 for j in [1..NumberOfColumns(Q)] do
                     s:= Qi[i][j] mod orders[j];
                     p:= pow( A, m, one, gg[j], s );
                     a:= Eltseq( (A!a)*(A!p) ); for o in [1..#m] do a[o]:= a[o] mod m[o]; end for;
                 end for;
                 Append( ~newgg, a );
                 Append( ~neword, S[i][i] );
              end if;
          end for;

          gg:= newgg;
          orders:= neword;

          G, exps:= enum_withexps( A, m, one, orders, gg );
       end if;
   end for;

   return gg, orders, G, exps;

end function;


order:= function( A, J, e )

    // basis of the Z-order of J = e*A, spanned by all e*a, a in a basis of A...

    m:= [ Coordinates( J, Basis(A)[i]*e ) : i in [1..Dimension(A)] ];
    d:= 1;
    for i in [1..#m] do
        for j in m[i] do
            d:= Lcm( d, Denominator( j ) );
        end for;
    end for;
    for i in [1..#m] do
        for j in [1..#m[i]] do
            m[i][j]*:= d;
        end for;
    end for;

    m,P:= HermiteForm( Matrix( Integers(), m ) );
    m:= Matrix( Rationals(), m )/d;
    P:= Matrix( Rationals(), P );

    O:= [ ];
    Aelms:= [ ];
    for i in [1..NumberOfRows(m)] do
        a:= &+[ m[i][j]*Basis(J)[j] : j in [1..Dimension(J)] ];
        if not IsZero(a) then 
           Append( ~O, a ); 
           Append( ~Aelms, &+[ Basis(A)[j]*P[i][j] : j in [1..Dimension(A)] ] );
        end if;
    end for;
 
    return O,Aelms;

end function;

preim:= function( f, a )

    // here f is a hom between two fp abelian groups, a in the codomain, we compute a preimage.

    A:= Domain( f );
    B:= Codomain(f);
    rk:= Ngens( B );
    m:= [ Eltseq( A.i@f ) : i in [1..Ngens(A)] ];
    for i in [1..Ngens(B)] do
        d:= Order(B.i);
        if d gt 0 then
           v:= [ 0 : j in [1..rk] ];
           v[i]:= d;
           Append( ~m, v );
        end if;
    end for;
    m:= Matrix(m);
    v:= Vector( Eltseq(a) );

    w:= Solution( m, v );
    w:= [ w[i] : i in [1..Ngens(A)] ];
    return A!Eltseq(w);

end function;

inv:= function( e, I, a )   // I ideal with idpot e, a elt; get b such that ab=e...

    lg:= function(u) 

       if u eq 0 then return u; end if;
       return Floor(Log(Abs(u))/Log(10));
    end function;

    ca:= Coordinates( I, a );
    ce:= Coordinates( I, e );
    n:= Dimension(I);
    eqns:= ZeroMatrix( Rationals(), n, n );
    for i in [1..n] do
        for j in [1..n] do
            cij:= Coordinates( I, Basis(I)[i]*Basis(I)[j] );
            for k in [1..n] do
                eqns[j][k] +:= cij[k]*ca[i];
            end for;
        end for;
    end for;

    w:= Eltseq( Solution( eqns, Vector( Rationals(), ce ) ) );
    return &+[ w[i]*Basis(I)[i] : i in [1..n] ];

end function;

extendunits:= function( A, e1, e2, I1, I2, u1, u2, ord1, ord2 )

    id:= e1+e2;
    J:= ideal< A | [id] >;
    O:= order( A, J, id );
    O1,Ae1:= order( A, I1, e1 );
    O2,Ae2:= order( A, I2, e2 );

    //next compute Oi \cap O:

    M:= KModule( Rationals(), Dimension(J) );
    VO:= KModuleWithBasis( [ M!Coordinates( J, u ) : u in O ] );
    c:= Matrix( Rationals(), [ Coordinates( VO, M!Coordinates(J,u) ) : u in O1 ] );
    c:= Saturation( c );
    b1:= [ &+[ c[i][j]*O[j] : j in [1..#O] ] : i in [1..#O1] ];

    c:= Matrix( Rationals(), [ Coordinates( VO, M!Coordinates(J,u) ) : u in O2 ] );
    c:= Saturation( c );
    b2:= [ &+[ c[i][j]*O[j] : j in [1..#O] ] : i in [1..#O2] ];

    bas:= b1 cat b2;

    m:= Matrix( Integers(), [ Coordinates( VO, M!Coordinates(J,u) ) : u in bas ] );
    S,P,Q:= SmithForm( m );

    offset:= 1;
    while S[offset][offset] eq 1 do 
          offset:= offset+1; 
          if offset gt NumberOfRows(S) then break; end if;
    end while;
    moduli:= [ S[i][i] : i in [offset..NumberOfColumns(S)] ];
    dim:= #moduli;
    T:= [ [] : i in [1..dim] ];

    Qi:= Q^-1;
    VQ:= KModuleWithBasis( [ M!Vector( Rationals(), Qi[i] ) : i in [1..NumberOfRows(Q)] ] );

    for i in [1..dim] do
        for j in [i..dim] do
            a:= &+[ Qi[ i+offset-1 ][k]*O[k] : k in [1..#O] ];
            b:= &+[ Qi[ j+offset-1 ][k]*O[k] : k in [1..#O] ];
            c:= Coordinates( VQ, M!Coordinates( VO, M!Coordinates(J,a*b) ) );
            c:= [ (Integers()!c[k]) mod moduli[k-offset+1] : k in [offset..#c] ];
            T[i][j]:= c;
            if i ne j then T[j][i]:= T[i][j]; end if;
        end for;
    end for;

    one:= Coordinates( VQ, M!Coordinates( VO, M!Coordinates(J,id) ) );
    one:= [ (Integers()!one[k]) mod moduli[k-offset+1] : k in [offset..#one] ];

    Alg:= Algebra< Integers(), #moduli | T >;

    // map gens of O1*, O2*...
    M1:= KModule( Rationals(), Dimension(I1) );
    VO1:= KModuleWithBasis( [ M1!Coordinates( I1, u ) : u in O1 ] );
    g1:= [ ];
    for u in u1 do
        cf:= Coordinates( VO1, M1!Coordinates( I1, u ) );
        a:= &+[ cf[i]*Ae1[i] : i in [1..#cf] ];
        a:= id*a;  // now in O...
        cf:= Coordinates( VQ, M!Coordinates( VO, M!Coordinates(J,a) ) );
        cf:= [ (Integers()!cf[k]) mod moduli[k-offset+1] : k in [offset..#cf] ];
        Append( ~g1, cf );
    end for;

    M2:= KModule( Rationals(), Dimension(I2) );
    VO2:= KModuleWithBasis( [ M2!Coordinates( I2, u ) : u in O2 ] );
    g2:= [ ];
    for u in u2 do
        cf:= Coordinates( VO2, M2!Coordinates( I2, u ) );
        a:= &+[ cf[i]*Ae2[i] : i in [1..#cf] ];
        a:= id*a;  // now in O...
        cf:= Coordinates( VQ, M!Coordinates( VO, M!Coordinates(J,a) ) );
        cf:= [ (Integers()!cf[k]) mod moduli[k-offset+1] : k in [offset..#cf] ];
        Append( ~g2, cf );
    end for;

    // make first a hom Oi^* --> group, then take the preimages of the intersection...
    // so we make the two image groups, that we enumerate anyway, get normal gens of those,
    // make the homs, get gens of the intersection, take preimages of the intersection
    // (preimages of the gens, plus the kernel?)

    gg1, orders1, grp1, exps1:= normalgens( Alg, moduli, one, g1 );

    if exps1[1] eq [] then exps1[1]:= [ Integers() | ]; end if;
    Ab:= FreeAbelianGroup( #orders1 );
    rr:= [ orders1[i]*Ab.i : i in [1..#orders1] ];
    G1:= quo< Ab | rr >;


    gg2, orders2, grp2, exps2:= normalgens( Alg, moduli, one, g2 );

    if exps2[1] eq [] then exps2[1]:= [ Integers() | ]; end if;
    Ab:= FreeAbelianGroup( #orders2 );
    rr:= [ orders2[i]*Ab.i : i in [1..#orders2] ];
    G2:= quo< Ab | rr >;

    
    if &and[ h in grp2 : h in gg1 ] then
       gg:= gg1; orders:= orders1; grp:= grp1; exps:= exps1;
    else
      grp:= SetToSequence( SequenceToSet( grp1 ) meet SequenceToSet( grp2 ) );
      gg:= gens( Alg, moduli, one, grp );
      gg, orders, grp, exps:= normalgens( Alg, moduli, one, gg );
    end if;

    Ab:= FreeAbelianGroup( #ord1 );
    rr:= [ ];
    for i in [1..#ord1] do
        if ord1[i] gt 0 then
           Append( ~rr, ord1[i]*Ab.i );
        end if;
    end for;
    Ab1:= quo< Ab | rr >;

    Ab:= FreeAbelianGroup( #ord2 );
    rr:= [ ];
    for i in [1..#ord2] do
        if ord2[i] gt 0 then
           Append( ~rr, ord2[i]*Ab.i );
        end if;
    end for;
    Ab2:= quo< Ab | rr >;

    im1:= [];
    for g in g1 do
        p:= Index( grp1, g );
        Append( ~im1, G1!exps1[p] );
    end for;

    f1:= hom< Ab1 -> G1 | im1 >;

    im2:= [];
    for g in g2 do
        p:= Index( grp2, g );
        Append( ~im2, G2!exps2[p] );
    end for;
    f2:= hom< Ab2 -> G2 | im2 >;

    a1:= [ Ab1!u : u in SetToSequence(Generators( Kernel( f1 ) )) ];
    for g in gg do
        p:= Index( grp1, g );
        Append( ~a1, preim( f1, G1!exps1[p] ) );  
    end for;
    A1:= sub< Ab1 | a1 >;
    a1:= [ A1.i : i in [1..Ngens(A1)] ];
    a1:= [ Ab1!u : u in a1 ];

    a2:= [ Ab2!u : u in SetToSequence(Generators( Kernel( f2 ) )) ];
    for g in gg do
        p:= Index( grp2, g );
        Append( ~a2, preim( f2, G2!exps2[p] ) );
    end for;
    A2:= sub< Ab2 | a2 >;
    a2:= [ A2.i : i in [1..Ngens(A2)] ];
    a2:= [ Ab2!u : u in a2 ];

    // now we need to make the maps < ai > --> G1 cap G2.

    Ab:= FreeAbelianGroup( #orders );
    rr:= [ orders[i]*Ab.i : i in [1..#orders] ];
    G:= quo< Ab | rr >;

    im1:= [];
    for g in a1 do
        p:= Index( exps1, Eltseq( g@f1 ) );
        h:= grp1[p];
        p:= Index( grp, h );
        Append( ~im1, G!exps[p] );
    end for;

    f1:= hom< A1 -> G | im1 >;
    ca1:= [ Eltseq( (A1.i)@f1 ) : i in [1..Ngens(A1)] ];

    im2:= [];
    for g in a2 do
        p:= Index( exps2, Eltseq( g@f2 ) );
        h:= grp2[p];
        p:= Index( grp, h );
        Append( ~im2, G!exps[p] );
    end for;

    f2:= hom< A2 -> G | im2 >;
    ca2:= [ Eltseq( (A2.i)@f2 ) : i in [1..Ngens(A2)] ];

    lattices:= [ ];
    len:= Ngens(A1)+Ngens(A2);

    for j in [1..Ngens(G)] do
        N:= Order( G.j );
        R:= Integers(N);
        m:= [ ];
        for i in [1..Ngens(A1)] do
            Append( ~m, [ R!(ca1[i][j]) ] );
        end for;
        for i in [1..Ngens(A2)] do
            Append( ~m, [ R!(-ca2[i][j]) ] );
        end for;
        K:= Kernel( Matrix( R, m ) );
        m:= [ Eltseq( Basis(K)[i] ) : i in [1..Dimension(K)] ];
        m:= Matrix( Integers(), m );
        m:= [ Eltseq(m[i]) : i in [1..NumberOfRows(m)] ];
        for i in [1..len] do
            v:= [ 0 : k in [1..len] ];
            v[i]:= N;
            Append( ~m, v );
        end for;
        m:= HermiteForm( Matrix( Integers(), m ) );
        m0:= [ ];
        for i in [1..NumberOfRows(m)] do
            if not IsZero(m[i]) then Append( ~m0, Eltseq(m[i]) ); end if;
        end for;
        Append( ~lattices, Lattice( Matrix( Integers(), m0 ) ) );
    end for;

    //Now get the intersection of the lattices...

    if #lattices gt 0 then
       L:= Basis( &meet(lattices) );
    else
       L:= IdentityMatrix( Integers(), len );
       L:= [ L[i] : i in [1..len] ];
    end if;

    s:= [ ];
    inv1:= [ inv(e1,I1,v) : v in u1 ];
    inv2:= [ inv(e2,I2,v) : v in u2 ];
    for u in L do
        c:= [ u[i] : i in [1..Ngens(A1)] ];
        a:= &+[ c[i]*a1[i] : i in [1..#c] ];
        c:= Eltseq( a );
        a:= e1;
        for i in [1..#c] do
            if c[i] gt 0 then
               a*:= u1[i]^c[i];
            elif c[i] lt 0 then
               a*:= inv1[i]^(-c[i]);
            end if;
        end for;

        c:= [ u[i] : i in [Ngens(A1)+1..len] ];
        b:= &+[ c[i]*a2[i] : i in [1..#c] ];
        c:= Eltseq( b );
        b:= e2;
        for i in [1..#c] do
            if c[i] gt 0 then
               b*:= u2[i]^c[i];
            elif c[i] lt 0 then
               b*:= inv2[i]^(-c[i]);
            end if;
        end for;

        Append( ~s, a+b );
    end for;

    return s;   
        
end function;


////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//  Generators of the arithmetic group corresponding to a toral algebraic group.
//

evecs:= function( F, mat, evals )

     n:= NumberOfRows(mat);
     ii:= IdentityMatrix( F, n );
     ev:= [ ];
     for x in evals do
        ev cat:= Basis( Kernel( mat-x*ii ) );
     end for;

     return ev;
end function;


splittingfield:= function( F, p )

     K:= F;
     R:= PolynomialRing(K);
p;
     ff:= Factorization( p );
     facs:= [ h[1] : h in ff ];
ff;
     mults:= [ h[2] : h in ff ];
        
     ready:= false;
     while not ready do
         found:= false; 
         i:= 1;

         while not found and i le #facs do
               if Degree( facs[i] ) gt 1 then
                  f:= facs[i];
                  Remove( ~facs, i );
                  m:= mults[i];
                  Remove( ~mults, i );
                  ff:= Factorization( PolynomialRing(K)!f );
ff;
                  if #ff eq 1 and ff[1][2] eq 1 then
                     found:= true;
                  else
                     facs cat:= [ e[1]: e in ff ];
                     mults cat:= [ m*e[2] : e in ff ];
                  end if;
               else
                  i:= i+1;
               end if;
         end while;

         if found then
f;
           K:= ext< K | R!f >;
           R:= PolynomialRing(K);
           facs:= [ R!h : h in facs ];
           ff:= Factorization( R!f );
           facs cat:= [ e[1]: e in ff ];
           mults cat:= [ m*e[2] : e in ff ];
         else
           ready:= true;
         end if;

     end while;

     cfs:= [ K!(-Coefficient(u,0)) : u in facs ]; 

     return K,cfs,mults;

end function;

eltseq:= function( noext, elm )


     c:= Eltseq( elm );
     for i in [1..noext] do
         cnew:= [];
         for x in c do
             cnew cat:= Eltseq( x );
         end for;
         c:= cnew; 
     end for;
     return c;

end function;


splitting_elt:= function( T )


  // first look for basis element:

  n:= Dimension(T);
  mdeg:= 0;
  for i in [1..n] do
      deg:= Degree( MinimalPolynomial( Basis(T)[i] ) );
      if deg eq Dimension(T) then
         return Basis(T)[i];
      else
  	 if deg gt mdeg then mdeg:= deg; k:=i; end if;
      end if;
  end for;

  a:= Basis(T)[k];
  A:= sub< T | [a] >;

  while true do
     if Dimension(A) eq Dimension(T) then
        return a;
     end if;
     for i in [1..n] do
	 if not Basis(T)[i] in A then    
	    b:= Basis(T)[i]; break;
         end if;
     end for;

     k:= 1;
     B:= sub< T | [k*a+b] >;
     while Dimension(B) le Dimension(A) do
	   k:= k+1;
           B:= sub< T | [k*a+b] >;
     end while;
     A:= B; a:= k*a+b;
  end while;
 
end function;

GensSimpleTorus:= function( T, ord )

  // here T is a matrix algebra that is a torus isomorphic to a field
  // extension of Q. ord is a basis of an order in T; we return generators
  // of the group of units in ord.

  if Dimension(T) eq 1 then 
     // check whether Basis(T)[1] is the identity...
     if not (Basis(T)[1])^2 eq Basis(T)[1] then
	print "HELP not one!!!!";
     end if;
     return [-Basis(T)[1] ], [2]; 
  end if;

  n:= Dimension(T);
  t0:= splitting_elt( T );

  M:= RMatrixSpace( Rationals(), Degree(T), Degree(T) );
  V:= RMatrixSpaceWithBasis( [ M!(t0^i) : i in [1..n] ] );
//print MinimalPolynomial(t0);
  F:= NumberField( MinimalPolynomial( t0 ) );

  bF:= [ F.1^i : i in [1..n] ];

  cbF:= [ Eltseq( u ) : u in bF ];
  trafo:= Matrix( Rationals(), cbF )^-1;

  cc:= [ Coordinates( V, M!u ) : u in ord ];
  R:= [ &+[ u[j]*bF[j] : j in [1..n] ] : u in cc ];

  O:= Order( R );

tm:= Cputime();
  G,f:= UnitGroup(O);
print "Dirichlet time....",Cputime()-tm;
  g:= [ Vector( Rationals(), Eltseq(F!f(G.i)))*trafo : i in [1..Ngens(G)] ];

  tt:= [ t0^i : i in [1..n] ];


  return [ &+[ u[j]*tt[j] : j in [1..n] ] : u in g ], [ Order(G.i) : i in [1..Ngens(G)] ];



end function;


GensSubTorus:= function( A, latt )

  // here A is a dxd matrix algebra, and a simple torus
  // latt is a basis of a full-dim lattice in a faithful A-submodule
  // of Q^d. 
  
  // get the image of A acting on latt.

  d:= Degree(A);
  V:= RModule( Rationals(), d );
  W:= RModuleWithBasis( [ V!v : v in latt ] );

  B:= [ ];
  for a in Basis(A) do
	  Append( ~B, Transpose( Matrix( Rationals(), 
					 [ Coordinates( W, V!v*Transpose(a) ) : v in latt ] ) ) );
  end for;

  m:= #latt;
  M:= MatrixAlgebra( Rationals(), m );
  BB:= sub< M | [ M!b : b in B ] >;
 
  mat:= [ ];
  for b in B do
      c:= Eltseq( b );
      lcm:= Lcm([ Denominator(r) : r in c ] );
      c:= [ lcm*u : u in c ];
      Append( ~mat, c );
  end for;

  mat:= Saturation( Matrix( Integers(), mat ) );

  ord:= [ Matrix( Rationals(), m, m, Eltseq(mat[i]) ) : i in [1..#B] ];

  g,o:= GensSimpleTorus( BB, ord );

  M:= RMatrixSpace( Rationals(), m, m );
  V:= RMatrixSpaceWithBasis( [ M!b : b in B ] );

  gg:= [ ];
  for x in g do
      c:= Coordinates( V, M!x );
      Append( ~gg, &+[ c[i]*Basis(A)[i] : i in [1..#B] ] );
  end for;

  return gg,o;


end function;

is_intmat:= function( m )

  return &and[ IsIntegral( u ) : u in Eltseq(m) ];

end function;

order:= function(m)  // order of a rational matrix...

    f:= MinimalPolynomial(m);
    fac:= Factorization( f );
    for a in fac do
        if a[2] gt 1 then
           return 0;
        elif not IsCyclotomicPolynomial(a[1]) then
           return 0;
        end if;
    end for;

    mk:= m; k:= 1;
    id:= m^0;
    while mk ne id do 
        mk *:= m;
        k +:= 1;
    end while;
    return k; 

end function;

GensTorus:= function( T )



  // here T is an associative matrix algebra isomorphic 
  // to a torus; we get generators of T^*_L, where L
  // is the lattice Z^n.

  
  e,A:= CentralIdempotents( T );
  deg:= Degree(T);
  gens:= [ ];
  for i in [1..#e] do
      d:= Lcm( [ Denominator(v) : v in Eltseq( e[i] ) ] );
      m:= HermiteForm( Matrix( Integers(), d*Transpose(e[i]) ) );
      latt:= [ ];
      for j in [1..deg] do
          if not IsZero(m[j]) then
             Append( ~latt, Vector( Rationals(), m[j] )/d );
          end if;
      end for;

      g,o:= GensSubTorus( A[i], latt );
      if i eq 1 then
         e1:= e[1]; I1:= A[1];
         u1:= g; ord1:= o;
      else
         e2:= e[i]; I2:= A[i]; u2:= g; ord2:= o;
         w0:= extendunits( T, e1, e2, I1, I2, u1, u2, ord1, ord2 );
         e1:= e1+e2;
         I1:= ideal< T | [e1] >; 
         u1:= w0;
         ord1:= [ order(a) : a in u1 ];     
      end if;
  end for;

  return u1;

end function;


lattice_complement:= function( bV, bW )

  // here bV, bW are bases of spaces in Q^n; we compute bases
  // of their intersection with Z^n, we return a basis of the
  // complement of W in V.

  // first we need integral coordinates:


  n:= #Eltseq( bV[1] );
  d:=  Lcm( [ Denominator(u) : u in Eltseq( Matrix(bV) ) ] );
  mV:= Matrix( Integers(), [ Eltseq( d*u ) : u in bV ] );

  if #bW eq 0 then
     mV:= Saturation( mV );
     return [ mV[i] : i in [1..NumberOfRows(mV)] ];
  end if;

  d:=  Lcm( [ Denominator(u) : u in Eltseq( Matrix(bW) ) ] );
  mW:= Matrix( Integers(), [ Eltseq( d*u ) : u in bW ] );

  mV:= Saturation( mV );
  mW:= Saturation( mW );

  sp:= RSpace( Rationals(), n );

  U:= RSpaceWithBasis( [ sp!mV[i] : i in [1..NumberOfRows(mV)] ] );

  mat:= [ Coordinates( U, sp!mW[i] ) : i in [1..NumberOfRows(mW)] ];

  S,P,Q:= SmithForm( Matrix( mat ) );
  Q:= Q^-1;

  bas:= [ ];
  for i in [NumberOfRows(mW)..NumberOfRows(Q)] do

      Append( ~bas, &+[ Q[i][j]*mV[j] : j in [1..NumberOfColumns(Q)] ] );

  end for;

  return bas;

end function;

minpol:= function( K, noext, a, x )


  vecs:= [ eltseq( noext, One(K) ) ];
  V:= RSpace( Rationals(), #vecs[1] ); 
  U:= sub< V | [ V!vecs[1] ] >;
  vecs:= [ V!vecs[1] ];
  done:= false;
  b:= One(K);
  while not done do
     b:= a*b;
     v:= V!eltseq(noext,b);
     if v in U then
        done:= true;
     else
        Append( ~vecs, v );
        U:= sub< V | vecs >; 
     end if;
  end while;

  U:= RSpaceWithBasis( vecs );
  cf:= Coordinates( U, v );

  f:= x^( #cf );
  for i in [1..#cf] do
      f:= f-cf[i]*x^(i-1);
  end for;
  return f;

end function;

multrels:= function(a)

   if &and[ x eq 1 : x in a ] then
      return IdentityMatrix( Integers(), #a );
   end if;

   if &and[ x in Rationals() : x in a ] then
      F:= CyclotomicField(4);
      a:= [ F!x : x in a];
   end if;

   H,m:= MultiplicativeGroup(a);
   F:= FreeAbelianGroup( #a );
   q:= hom< F -> H | [ x@@m : x in a ] >;
 
   K:= Kernel(q);

   return HermiteForm( 
            Matrix( Integers(), [ Eltseq(F!K.i) : i in [1..Ngens(K)] ] ) );

end function;

GensIntTorus:= function( L )

  // L : Lie algebra given by a basis consisting of matrices;
  // output: generators of the integral torus of the corr alg grp.

  
  deg:= NumberOfRows( L[1] );
  A:= MatrixAlgebra( Rationals(), deg );

  gg:= [ A!u : u in L ];
  Append( ~gg, One(A) );

  T:= sub< A | gg >;
tm0:= Cputime();
  gT:= GensTorus( T );  

print "Full torus time...",Cputime()-tm0; 
  t0:= splitting_elt( T );

tm0:= Cputime();
  //K,evl,mul:= splittingfield( Rationals(),  CharacteristicPolynomial(t0) );
  K,evl:= SplittingField( CharacteristicPolynomial(t0) );
print "Splitting field time...",Cputime()-tm0;
  noext:= 0; 
  dK:= #Basis(K);
  F:= BaseField(K);
  while not (F eq Rationals()) do
      dK *:= #Basis(F);
      F:= BaseField(F);
      noext +:=1;
  end while;
print "dimension splitting field...",dK;
  // L diagonalises over K, get the corr transformation...X

  X:= evecs( K, t0, evl );
  X:= Matrix( K, [ Eltseq(u) : u in X ] );
  Xi:= X^-1;

  bL:= [ X*Matrix( K, u )*Xi : u in L ];

  rr:= [ ];
  for i in [1..deg] do
      row:= [ ];
      for j in [1..#bL] do
	  row cat:= eltseq( noext, bL[j][i][i] );
      end for;
      Append( ~rr, row );
  end for;

  lam_g:= Basis( Kernel( Matrix( Rationals(), rr ) ) );
  if #lam_g eq 0 then
     e:= [ Vector( [ 0 : i in [1..deg] ] ) ];
  else
    // same for the torus
    bT:= [ X*Matrix( K, u )*Xi : u in Basis(T) ];

    rr:= [ ];
    for i in [1..deg] do
      row:= [ ];
      for j in [1..#bT] do
	  row cat:= eltseq( noext, bT[j][i][i] );
      end for;
      Append( ~rr, row );
    end for;

    sp:= Kernel( Matrix( Rationals(), rr ) );
    // we need a basis of the complement of sp in the space spanned by e:

    e:= lattice_complement( lam_g, Basis(sp) );

  end if; 

  // we get the corr character values on the generators computed...


  T_traf:= [ X*u*Xi : u in gT ];

  ch_vals:= [ ];
  for i in [1..#e] do
      ch:= [];
      for u in T_traf do
          Append( ~ch, &*[ u[j][j]^e[i][j] : j in [1..deg] ] );
      end for;
      Append( ~ch_vals, ch );
  end for;
tm0:= Cputime();

  P<var>:= PolynomialRing( Rationals() );

  //ml:= multrels( K, ch_vals[1], [ minpol(K,noext,a,var) : a in ch_vals[1] ] );
  ml:= multrels( ch_vals[1] );

  L:= Lattice( ml );
  for i in [2..#ch_vals] do
	//ml:= multrels( K, ch_vals[i], [ minpol(K,noext,a,var) : 
	//					a in ch_vals[i] ] );
      ml:= multrels( ch_vals[i] );
      L:= L meet Lattice( ml );
  end for;
print "Ge time...",Cputime()-tm0;

  g:= [ ];
  for u in Basis( L ) do

      v:= Eltseq(u);
      Append( ~g, &*[ gT[i]^v[i] : i in [1..#v] ] );
  end for;

  // get the characters that define G, i.e., the intersection of their
  // kernels when restricted to G is 1.

  zm:= [ ];
  for i in [1..deg] do
      v:= [ 0 : k in [1..deg] ];
      v[i]:= 1;  
      Append( ~zm, Vector(v) );
  end for;

  f:= lattice_complement( zm, lam_g );

  T_traf:= [ X*u*Xi : u in g ];

  ch_vals:= [ ];
  for i in [1..#f] do
      ch:= [];
      for u in T_traf do
          Append( ~ch, &*[ u[j][j]^f[i][j] : j in [1..deg] ] );
      end for;
      Append( ~ch_vals, ch );

  end for;

tm0:= Cputime();

  //ml:= multrels( K, ch_vals[1], [ minpol(K,noext,a,var) : a in ch_vals[1] ] );
  ml:= multrels( ch_vals[1] );

  L:= Lattice( ml );
  for i in [2..#ch_vals] do
	//ml:= multrels( K, ch_vals[i], [ minpol(K,noext,a,var) : 
        //                                        a in ch_vals[i] ] ); 
      ml:= multrels( ch_vals[i] );
      L:= L meet Lattice( ml );
  end for;
print "Ge time...",Cputime()-tm0;

  if Dimension(L) eq 0 then 
     m:= Matrix( [ Eltseq(Zero(L)) ] );
  else
     m:= Matrix( Integers(), [ Eltseq(u) : u in Basis(L) ] );
  end if;
  S,P,Q:=SmithForm(m);

  e:= [ ];
  Q:= Q^-1; 
  d:= [ ];
  for i in [1..NumberOfRows(Q)] do
      if i le NumberOfRows(S) then
         if S[i][i] gt 1 then
   	    Append( ~d, S[i][i] );
            Append( ~e, Q[i] );
         end if;
      else
 	 Append( ~d, 0 );
         Append( ~e, Eltseq(Q[i]) );
      end if;
  end for;

  gg:= [ ];
  for v in e do
      u:= Eltseq(v); 
      Append( ~gg, &*[ g[i]^u[i] : i in [1..#u] ] );
  end for;

print "Group structure...",d;
  return gg;

end function;


///////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//  Unipotent groups...
//
//

trans:= function( L, W )

  // here L is a list of nxn matrices; W is a basis of a subspace of F^n;
  // we return the subspace of v\in F^n such that xv in W for all x in L.

  n:= NumberOfColumns( L[1] );
  m:= #L;
  r:= #W;
  F:= BaseRing( L[1] );
  zero:= Zero(F);

  eqns:= [ ];
  for j in [1..m] do
      for k in [1..n] do
          eqn:= [ zero : i in [1..n+r*m] ];
          for i in [1..n] do
              eqn[i]:= L[j][k][i];
          end for;
          for l in [1..r] do
              eqn[n+j+(l-1)*m] := -W[l][k];
          end for;
          Append( ~eqns, eqn );
      end for;
  end for;

  eqns:= Transpose( Matrix( F, eqns ) );

  b:= Basis( Kernel( eqns ) );

  bas:= [ ];
  for u in b do
      w:= Vector( F, [ u[i] : i in [1..n] ] );
      if not IsZero(w) and not w in bas then Append( ~bas, w ); end if;
  end for;

  V:= KModule( F, n );
  bas:= Basis( sub< V | [V!u : u in bas ] >);
  return [V!u : u in bas];

end function;

MaximalFlag:= function( L )

  // Here L is a list of matrices that forms a basis of a
  // Lie algebra consisting entirely of nilpotent matrices.
  // Let V be its natural module; this function returns a 
  // basis b of V with the following property:
  // X*b[i] lies in the span of b[1],..,b[i-1] for all X in K.

  F:= BaseRing( L[1] );

  n:= NumberOfRows( L[1] );
  V:= KModule( F, n );
  bas:= [ ];
  W:= sub< V | [] >;
  sp:= [ Zero(V) ];
  while #bas lt n do
      sp0:= trans( L, sp );
      for u in sp0 do
          if not u in W then 
             Append( ~bas, u );
             W:= sub< V | bas >; 
          end if;
      end for;
      sp:= sp0;
   end while;
    
   return bas;

end function;

Flag:= function( L )

  // Here L is a list of matrices that forms a basis of a
  // Lie algebra consisting entirely of nilpotent matrices.
  // Let V be its natural module; this function returns a 
  // list of bases of subspaces of V such that X*v lies in V_{i-1} for
  // all v in V_i

  F:= BaseRing( L[1] );

  n:= NumberOfRows( L[1] );
  V:= KModule( F, n ); 
  bas:= [ ];
  bb:= [ ];
  W:= sub< V | [ ] >; 
  sp:= [ Zero(V) ];
  while Dimension(W) lt n do
      sp0:= trans( L, sp );
      b0:= [ ];
      for u in sp0 do
          u0:= V!Eltseq(u);
          if not u0 in W then 
             Append( ~bas, u0 );
             Append( ~b0, u0 );
             W:= sub< V | bas >; 
          end if;
      end for;
      Append( ~bb, bas );
      sp:= sp0;
   end while;
    
   return bb;

end function;

saturation:= function( A )


     // we assume that A has lin indep rows, and is an mxn mat
     // over Q; we return a basis of Z^n \cap V, where V is the
     // space spanned by the rows of A
     // we also return the inverse of the matrix Q.
     
     A0:= Matrix( Integers(), Matrix(A)*Lcm( [ Denominator(u) : u in Eltseq(A)] ) );
     S,P,Q:= SmithForm( A0 );
     for i in [1..NumberOfRows(S)] do
         if S[i][i] eq 0 then
            print "ZERO IN SATURATION\n";
         end if;
         S[i][i]:= 1;
     end for;
     Qi:= Q^-1;
     return P^-1*S*Qi, Qi;


end function;


intersection_latt_vs:= function( A, B )

   // Here A is an nxn integral matrix whose rows span a full dim lattice L
   // B is an mxn rational matrix whose rows span a subspace W of Q^n;
   // we return an nxn integral matrix whose rows span L, and whose
   // first m rows span W cap L.

   n:= NumberOfRows(A);
   V:= KModule( Rationals(), n );
   BL:= KModuleWithBasis( [ V!A[k] : k in [1..n] ] );

   B0:= Matrix( Rationals(), [ Coordinates( BL, B[k] ) : k in [1..NumberOfRows(B)] ] );
   _,Qi:= saturation( B0 );
   mat:= Qi*A;

   M:= HermiteForm( Matrix( Integers(), [ Eltseq(mat[k]) : k in [1..NumberOfRows(B)] ] ) );
   M:= [ M[k] : k in [1..NumberOfRows(M)] ];
   M cat:= [ mat[k] : k in [NumberOfRows(B)+1..NumberOfRows(mat)] ];
   return Matrix( Integers(), [ Eltseq(u) : u in M ] );

end function;

intersection_latt_vs_idmat:= function( B )

   // Same as previous function; but assuming A is the identity...


   _,mat:= saturation( B );
   M:= HermiteForm( Matrix( Integers(), [ Eltseq(mat[k]) : k in [1..NumberOfRows(B)] ] ) );
   M:= [ M[k] : k in [1..NumberOfRows(M)] ];
   M cat:= [ mat[k] : k in [NumberOfRows(B)+1..NumberOfRows(mat)] ];
   return Matrix( Integers(), [ Eltseq(u) : u in M ] );

end function;

Lcomplements:= function( A, V1, V2 )
   
   if IsIdentity(A) then
      W:= intersection_latt_vs_idmat( Matrix(V2) );
   else
      W:= intersection_latt_vs( A, Matrix(V2) );
   end if;

   t:= #V2;
   w:= [ W[i] : i in [1..t] ];
   n:= NumberOfRows(A);

   V:= KModule( Rationals(), n );
   bas:= KModuleWithBasis( [ V!Eltseq(u) : u in w ] ); 
   A0:= Matrix( Rationals(), [ Coordinates( bas, V!Eltseq(x) ) : x in V1 ] );

   B:= intersection_latt_vs_idmat( A0 );
   C:= B*Matrix( Integers(), [ Eltseq(u) : u in w ] );
   C:= [ Eltseq(C[i]) : i in [1..NumberOfRows(C)] ];
   C cat:= [ Eltseq(W[i]) : i in [t+1..NumberOfRows(W)] ];
   return Matrix( Integers(), C );

end function;

integral_rels:= function( A )

   // here A is an mxn matrix with rational coefficients
   // we output a basis of the lattice of (e1,...,em)\in \Z^m
   // with e1A[1] + ... + emA[m]\in \Z^n 

   n:= NumberOfColumns( A );
   A0:= [ Eltseq(A[i]) : i in [1..NumberOfRows(A)]];
   I0:= IdentityMatrix( Rationals(), n );
   A0 cat:= [ Eltseq(I0[i]) : i in [1..n] ]; 
   sol:= Matrix( Rationals(), [ Eltseq(u) : u in Basis(Kernel( Matrix( Rationals(), A0 ) )) ] );
   B:= saturation( sol );
   m:= NumberOfRows(A);
   B:= [ [ B[i][j] : j in [1..m] ] : i in [1..NumberOfRows(B)] ];
   return HermiteForm( Matrix( Integers(), B ) );
  
end function;



GensIntUnipot:= function( g, L )


    // Here g is a list of nilpotent matrices, spanning a Lie algebra;
    // with unipotent group G.
    // L is a basis of a lattice on which g acts. 
    // We return a list of generators of the corr arithmetic group;
    // they form a polycyclic series.

    exp:= function( mat )

        // mat assumed nilpotent...
        // Note: this can be done much more efficiently...
        
        mat:= Matrix( Rationals(), mat );
        return &+[ mat^i/Factorial(i) : i in [0..NumberOfRows(mat)] ];

    end function;

    log:= function( mat )

       // mat is assumed unipotent...
       
       A:= mat-mat^0;
       return &+[ (-1)^(i+1)*A^i/i : i in [1..NumberOfRows(A)] ];

    end function;

    if #g eq 0 then
       // just the identity
       return [ ];
    end if;

print "DIM  ",NumberOfRows(L);

    m:= NumberOfRows(L);

    fl:= Flag( g );
    n:= #fl;
//print "Len Flag  ",n,"\n";
//fl[1]; fl[n-1];
    Lc:= Lcomplements( L, fl[1], fl[n-1] );

//print Lc;
    k:= #fl[1];
    l:= #fl[n-1];

    W_1:= [ Lc[i] : i in [l+1..m] ];
    L1:= [ Lc[i] : i in [1..l] ];
    
    V:= KModule( Rationals(), #Eltseq(L[1]) );
    B1:= KModuleWithBasis( [ V!Eltseq(u) : u in L1 ] );

    L2:= [ Lc[i] : i in [k+1..m] ];
    B2:= KModuleWithBasis( [ V!Eltseq(Lc[i]) : i in [1..NumberOfRows(Lc)] ] );

    g0:= [ ];
    n1:= #L1;
    n2:= #L2;

    for x in g do
        m:= ZeroMatrix( Rationals(), n1+n2, n1+n2 );
        for i in [1..n1] do

            cf:= Coordinates( B1, V!( Eltseq( (V!Eltseq(L1[i]))*Transpose(x) ) ) );
            for j in [1..#cf] do
                m[j][i]:= cf[j];
            end for;
        end for;

        for i in [1..n2] do
            cf:= Coordinates( B2, V!Eltseq( (V!Eltseq(L2[i]))*Transpose(x) ) );
            cf:= [ cf[j] : j in [k+1..#cf] ];
            for j in [1..#cf] do
                m[n1+j][n1+i]:= cf[j];
            end for;
        end for;
      
        Append( ~g0, m );

    end for;

    s:= NumberOfRows( g[1] );
    Ms:= KMatrixSpace( Rationals(), s, s );
    K:= KMatrixSpaceWithBasis( [ Ms!u : u in g ] );

    f:= function( mt ) 
         cf:= Coordinates( K, Ms!mt );
         return &+[ cf[i]*g0[i] : i in [1..#cf] ];
    end function;
    
    M12:= KMatrixSpace( Rationals(), n1+n2, n1+n2 );
    sp:= sub< M12 | [ ] >;
    b12:= [M12| ];
    inds:= [ ];
    for i in [1..#g0] do
        if not M12!g0[i] in sp then
           Append( ~inds, i );
           Append( ~b12, M12!g0[i] );
           sp:= sub< M12 | b12 >;
        end if;
    end for;

    sp:= KMatrixSpaceWithBasis( b12 );
    fi:= function( mt )

         cf:= Coordinates( sp, M12!mt );
         return &+[ cf[i]*g[inds[i]] : i in [1..#cf] ];
    end function;
           
    g00:= b12;

//print "G00  ",g00;

    eqns:= [ ];

    for x in g do
        eqn:= [ ];
        for j in [1..n1] do
            eqn cat:= Eltseq( V!Eltseq(L1[j])*Transpose(x) );
        end for;
        Append( ~eqns, eqn );
    end for;
    sol1:= Kernel( Matrix(eqns) );
    
    eqns:= [ ];
    for x in g do
        eqn:=[ ];
        for j in [1..n2] do
            cf:= Coordinates( B2, V!Eltseq( V!Eltseq(L2[j])*Transpose(x) ) );
            cf:= [ cf[r] : r in [k+1..#cf] ];
            eqn cat:= cf;  
        end for;
        Append( ~eqns, eqn );
    end for;
    sol2:= Kernel( Matrix(eqns) );
    if Dimension(sol1) eq 0 then
       sol:= Basis(sol2);
    elif Dimension(sol2) eq 0 then
       sol:= Basis(sol1);
    else
       sol:= Basis( sol1 meet sol2 );
    end if;

    n0:= [ &+[ x[i]*g[i] : i in [1..#g] ] : x in sol ];


    // NB: as we took coefficients of the Lie alg elms wrt a lattice basis;
    // the image lattice is just the standard one...
//print g00,"\n\n";
    r1:= $$( g00, IdentityMatrix( Rationals(), n1+n2 ) );
//print "r1 ",#r1,"\n\n";

     s0:= k;
     t:= #W_1;  

    lam:= function( mat )

        m:= ZeroMatrix( Rationals(), s0, t );
        for i in [1..t] do
            cf:= Coordinates( B2, V!Eltseq(V!Eltseq(W_1[i])*Transpose(mat)) );
            cf:= [ cf[l] : l in [1..s0] ];
            for j in [1..s0] do
                m[j][i]:= cf[j];
            end for;
        end for;
        return m;

    end function;

    lamN:= [ lam(x) : x in n0 ]; 

    Ms0t:= KMatrixSpace( Rationals(), s0, t );
    BasLamN:= KMatrixSpaceWithBasis( [ Ms0t | Ms0t!u : u in lamN ] );
    // needed for inverse map...

    Alg:= MatrixAlgebra( Rationals(), NumberOfRows( g[1] ) );
    A0:= sub< Alg | [ Alg!u : u in g ] >;
    
    ass:= sub< Ms0t | [ lam(x) : x in Basis(A0) ] >;
//print "DIM ASS ",Dimension(ass),"\n";
    ba:= Basis(ass);
//ba;
    inds:= [ ];
    for i in [1..s0] do
        for j in [1..t] do
            if not &and[ IsZero(x[i][j]) : x in ba ] then
               Append( ~inds, [i,j] );
            end if;
        end for;
    end for;

//print "LEN inds  ",#inds,"\n";

    
    // we work with the space of s0xt matrices that have
    // nonzero entries only in the places contained in inds.
    
    // make a basis of s0xt mats that is a complement to
    // lamN...

    if #lamN eq 0 then
       basC:= [ ];
       for x in inds do
           i:= x[1]; j:= x[2];
           m:= ZeroMatrix( Rationals(),s0,t );
           m[i][j]:= 1;
           Append( ~basC, m );
       end for;
    else

       lamN0:= [ [ ] : u in lamN ];
       for i in [1..#lamN0] do
           for x in inds do
               Append( ~lamN0[i], lamN[i][x[1]][x[2]] );
           end for;
       end for;

//print "LEN lamN0  ", #lamN0[1],"\n";

       latt:= intersection_latt_vs_idmat( Matrix( lamN0 ) );
       basC:= [ ];
       lamN:= [ ];
       for r in [1..#lamN0] do
           m:= ZeroMatrix( Rationals(), s0, t );
           for l in [1..#inds] do
               i:= inds[l][1]; j:= inds[l][2];
               m[i][j]:= latt[r][ l ];
           end for;
           Append( ~lamN, m );
       end for;
       for r in [#lamN0+1..NumberOfRows(latt)] do
           m:= ZeroMatrix( Rationals(), s0, t );
           for l in [1..#inds] do
               i:= inds[l][1]; j:= inds[l][2];
               m[i][j]:= latt[r][ l ];
           end for;
           Append( ~basC, m );
       end for;
 
    end if;
    
    mC:= #basC;
//print "mC  ",mC,"\n";
    if mC eq 0 or #r1 eq 0 then
       // rather trivial case...
       k_g:= r1; // i.e. everything is in the kernel...
       if mC eq 0 then
          BC:= KMatrixSpaceWithBasis( [ Ms0t!u : u in lamN ] );
       end if;
    else
       basC cat:= lamN;
       BC:= KMatrixSpaceWithBasis( [ Ms0t!u : u in basC ] ); 

       imgs:= [ ]; // will be the images of the gens of G0 under 
                   // Phi_Q

    
       for x in r1 do     

           u:= lam( exp( fi( log(x) ) ) );
           cf:= Coordinates( BC, Ms0t!u );
           Append( ~imgs, [ cf[i] : i in [1..mC] ] );

       end for;

//print "imgs  ",imgs,"\n";

       sol:= integral_rels( Matrix(imgs) );
//print "sol  ",sol,"\n";

       k_g:= [ &*[ r1[i]^sol[j][i] : i in [1..#imgs] ] : j in [1..NumberOfRows(sol)] ];

    end if;

    k_g:= [ exp( fi( log(x) ) ) : x in k_g ];

    // This might not be integral; we find elements of N that
    // make it integral.

    for i in [1..#k_g] do
        cf:= Coordinates( BC, Ms0t!lam( k_g[i] ) );
        cf:= [ cf[j] : j in [mC+1..#cf] ];
        if #cf gt 0 then
           x:= &+[ cf[j]*lamN[j] : j in [1..#cf] ];
           cf:= Coordinates( BasLamN, Ms0t!x );
//if not &and[ x eq 0 : x in cf ] then print "CFT!!  ",cf,"\n"; end if;
           k_g[i]:= k_g[i]*exp( -&+[ cf[j]*n0[j] : j in [1..#cf] ] );
        end if;
    end for;
 
    // for N we take inverse images, as we already found
    // a lattice basis for N...

//print "lamN  ",#lamN,"\n";

    n00:= [ ];
    for x in lamN do
        cf:= Coordinates( BasLamN, Ms0t!x );        
        Append( ~n00, &+[ cf[j]*n0[j] : j in [1..#cf] ] );
    end for;

//print "done\n";
    
    k_g cat:= [ exp(u) : u in n00 ];  
                                     
    return k_g;


end function;

intrinsic GeneratorsOfArithmeticGroup( L::SeqEnum ) -> SeqEnum
{L is a basis of an algebraic matrix Lie algebra over Q; this returns generators of a corresponding arithmetic group if the Lie algebra is solvable}


    s,t,n:= DecompositionOfAlgebraicLieAlgebra( L );
    require #s eq 0 : "Can only do this for solvable Lie algebras";

    if #t gt 0 then
       a:= GensIntTorus( t );
    else
       a:= [ ];
    end if;
    if #n gt 0 then
       b:= GensIntUnipot( n, IdentityMatrix( Integers(), NumberOfRows(L[1]) ) );
    else
       b:= [ ];
    end if;
    return a cat b;
    

end intrinsic;

