/*
This file contains an implementation of a MAGMA function 

       UnitsOfAbelianGroupRing( o )

where o is a sequence of integers. If o = [ d1, d2, ..., dk ] then the unit group of 
the group ring over Z of the group C_{d1} x C_{d2} x ... x C_{dk} is constructed.
It is assumed that di divides d_{i+1}.

The function returns four values: G, A, u, ord

G   : is the abelian group, constructed as a polycyclic group
A   : is its group algebra over Q
u   : generating set for the unit group (ZG)*
ord : list of integers, giving the multiplicative orders of the elements of u, a 0 means that
      the order is infinity.

Example:

> Attach("cycunits.m");
> Attach("grpalgunits.m");
> G,A,u,ord:= UnitsOfAbelianGroupRing( [2,10] );
> u;
[ G.1 * G.2 * G.3^2, -G.1, -G.1 * G.2 * G.3^4, 6*Id(G) + G.3 - 4*G.3^2 - 4*G.3^3 + G.3^4 - 4*G.2 + 2*G.2 * G.3 + 4*G.2 * 
G.3^2 + 2*G.2 * G.3^3 - 4*G.2 * G.3^4 + 5*G.1 + G.1 * G.3 - 3*G.1 * G.3^2 - 3*G.1 * G.3^3 + G.1 * G.3^4 - 4*G.1 * G.2 + 2*G.1
* G.2 * G.3 + 4*G.1 * G.2 * G.3^2 + 2*G.1 * G.2 * G.3^3 - 4*G.1 * G.2 * G.3^4, Id(G) - G.3 - G.3^2 + G.3^3 + G.3^4 + 2*G.2 * 
G.3 - G.2 * G.3^3 - G.2 * G.3^4, 4*Id(G) - G.3 - G.3^2 - G.3^3 - G.3^4 + 2*G.2 - 2*G.2 * G.3 - 2*G.2 * G.3^3 + 2*G.2 * G.3^4 
+ 2*G.1 * G.3 - 2*G.1 * G.3^2 - 2*G.1 * G.3^3 + 2*G.1 * G.3^4 + G.1 * G.2 + G.1 * G.2 * G.3 - 3*G.1 * G.2 * G.3^2 + G.1 * G.2
* G.3^3 + G.1 * G.2 * G.3^4, -25*Id(G) + 11*G.3 + 29*G.3^2 + 11*G.3^3 - 25*G.3^4 + 8*G.2 - 24*G.2 * G.3 - 24*G.2 * G.3^2 + 
8*G.2 * G.3^3 + 32*G.2 * G.3^4 + 27*G.1 - 9*G.1 * G.3 - 36*G.1 * G.3^2 - 9*G.1 * G.3^3 + 27*G.1 * G.3^4 - 12*G.1 * G.2 + 
28*G.1 * G.2 * G.3 + 28*G.1 * G.2 * G.3^2 - 12*G.1 * G.2 * G.3^3 - 32*G.1 * G.2 * G.3^4 ]
> ord;
[ 2, 2, 10, 0, 0, 0, 0 ]


*/

idpots:= function( A, G )

    ct:= CharacterTable( G );
    ct:= [ Eltseq(ct[i]) : i in [1..#ct] ];
    cg:= ConjugacyClasses(G);
    eG:= [ A!u[3] : u in cg ];
    
    K:= Parent( ct[1][1] );
    n:= Conductor(K);
    F<w>:= CyclotomicField( n );
    _,f:= IsIsomorphic( K, F );
    ct:= [ [ ct[i][j]@f : j in [1..#ct[i]] ] : i in [1..#ct] ];

    // now we need to group the elts of ct in Galois orbits...

    gal:= [ ];
    for i in [1..n] do
        if Gcd(i,n) eq 1 then
           Append( ~gal, hom< F -> F | [w^i] > );
        end if;
    end for;
    
    orbs:= [ ];
    dones:= [ ];
    for k in [1..#ct] do
        if not ct[k] in dones then
           orb:= [ ];
           for f in gal do
               u:= [ z@f : z in ct[k] ];
               if not u in orb then
                  Append( ~orb, u );
               end if;
           end for;
           Append( ~orbs, orb );
           dones cat:= orb;
        end if;
    end for;
        
    osums:= [ &+[ Vector( F, u ) : u in orbs[i] ] : i in [1..#orbs] ];
    ids:= [ ];
    for i in [1..#osums] do     
        e:= Zero(A);
        for k in [1..#G] do
            pos:= Position( eG, eG[k]^-1 );
            e +:= osums[i][pos]*eG[k];
        end for;
        Append( ~ids, e/#G );
    end for;

    return ids;

end function;

splitting_elt:= function( T )

  n:= Dimension(T);

  for i in [1..n] do
      A:= sub< T | [ Basis(T)[i] ] >;    
      if Dimension(A) eq Dimension(T) then
         return Basis(T)[i];
      end if;
  end for;
 
  a:= Basis(T)[1];
  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;



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;


grpalg_data:= function( G )


    A:= GroupAlgebra( Rationals(), G );
    e:= idpots( A, G );
    I:= [ ideal< A | [u] > : u in e ];

    eG:= [ g : g in G ];
    iscyc:= false;
    for g in eG do
       if Order(g) eq #G then
          iscyc:= true; c0:= A!g; break;
       end if;
    end for;

    if not iscyc then       
       c0:= splitting_elt( A );
    end if;
    mc0:= Matrix( [ Coordinates( A, c0*A.i ) : i in [1..Dimension(A)] ] );

    phi:= [ EulerPhi( k ) : k in [1..#G] ];

    F:= [ ];  // the cyclotomic fields,
    u:= [ ];  // the generators: u[k] is the pre-image of w in I[k].

    for k in [1..#I] do

        B:= I[k];
        if Dimension(B) eq 1 then
           Append( ~F, 1 );
           Append( ~u, One(B) );
        else
           found:= false;
           if iscyc then
              c:= c0*e[k];
           else
              c:= splitting_elt( B );
           end if;
           mc:= Transpose( Matrix( [ Coordinates( B, c*B.i ) : i in [1..Dimension(B)] ] ) );
           f:= MinimalPolynomial( mc );

           K:= NumberField( f );
           d:= Degree(K);
           s:= 0;
           for i in [1..#phi] do 
               F0<w0>:= CyclotomicField( i );
               bl,f0:= IsIsomorphic( F0, K );
               if bl then
                  s:= i;
                  w:= Eltseq(f0(w0));
                  break;
               end if;
           end for;

           if s eq 0 then "ERROR...no cyclotomic field found!"; end if;

           Append( ~F, s );
           q:= Zero(B);
           for i in [1..#w] do 
               if i eq 1 then
                  q +:= w[i]*One(B);
               else
                  q +:= w[i]*c^(i-1);
               end if;
           end for;
           Append( ~u, q );
        end if;  
    end for;

    found:= false;

    func:= function(a) return a; end function;

    units:= [ ];
    orders:= [ ];
    projs:= [ ];
    mps:= [ ];
  
    for i in [1..#F] do
        if F[i] eq 1 then

  	   a:= -A!One( I[i] );
           Append( ~units, [a] );
           Append(~orders, [2] );

           ffi:= function( elm )
          
              b:= e[i]*elm;
              if b eq -One(I[i]) then
                 return [1];
              else
                 return [0];
              end if;
           end function;
           Append( ~projs, ffi );

           mp:= function( a )
             return Coordinates( I[i], e[i]*a )[1]/Coordinates( I[i], e[i] )[1];
           end function;
 
           mpi:= function( b )
             return b*e[i];
           end function;
           Append( ~mps, [ mp, mpi ] );

        else
         
           K,v:= UnitsOfCyclotomicField( F[i] );
           Grp,map:= MultiplicativeGroup(v);
           Append( ~orders, [ Order(Grp.i) : i in [1..Ngens(Grp)] ] );
           v:= [ map( Grp.i ) : i in [1..Ngens(Grp)] ];
           uu:= [ ];
           for j in [1..#v] do
               cf:= Eltseq(v[j]);
               a:= Zero(A);
               for k in [1..#cf] do
                   if k eq 1 then
                      a:= a + cf[k]*One(I[i]);
                   else
                      a:= a + cf[k]*u[i]^(k-1);
                   end if;
               end for;
               Append( ~uu, a );
           end for;
           
           Append( ~units, uu );

           elms:= [ Coordinates( A, One(I[i]) ) ] cat [ Coordinates( A, u[i]^(k-1) ) : k in [2..Degree(K)] ]; 
           V:= KModule( Rationals(), Dimension(A) );
           V0:= KModuleWithBasis( [ V!s : s in elms ] );
           
           ffi:= function( a ) // a in A project it onto K, write as a product of units...
    
               cf:= Coordinates( V0, V!Coordinates( A, e[i]*a ) );
               b:= K!cf;
               vv:= [b] cat v;
               H,m:= MultiplicativeGroup( vv );
               FA:= FreeAbelianGroup( #vv );
               q:= hom< FA -> H | [ x@@m : x in vv ] >;
               Nuc:= Kernel(q);
               cf:= Eltseq( FA!(Nuc.1) );
               return [ -cf[i] : i in [2..#cf] ];
           end function;
           Append( ~projs, ffi );

           mp:= function( a ) // a in A, project to K

              cf:= Coordinates( V0, V!Coordinates( A, e[i]*a ) );
              return K!cf;

           end function;

           mpi:= function( b ) // b in K, map to A (ie to I[i]\subset A)

              cf:= Eltseq(b);
              a:= Zero(A);
              for k in [1..#cf] do
                  if k eq 1 then
                     a:= a + cf[k]*One(I[i]);
                  else
                     a:= a + cf[k]*u[i]^(k-1);
                  end if;
              end for;
              return a;
          end function;

          Append( ~mps, [ mp, mpi ] );

        end if;
    end for;

    return A,units,orders,I,e,func,projs,F,u,mps;


end function;

grpalg_func:= function( G )

    A:= GroupAlgebra( Rationals(), G );
    c:= splitting_elt( A );
    mc:= Matrix( [ Coordinates( A, c*A.i ) : i in [1..Dimension(A)] ] );
    N:= Exponent( G );
    K<w>:= CyclotomicField( N );
    P<x>:= PolynomialRing( K );
    h:= P!CharacteristicPolynomial(mc);

    r:= Roots(h);
    r:= [ u[1] : u in r ];
    X:= evecs( K, mc, r );
    X:= Matrix( K, [ Eltseq(u) : u in X ] );
    Xi:= X^-1;

    func:= function( a )  // a in A

       ma:= Matrix( [ Coordinates( A, a*A.i ) : i in [1..Dimension(A)] ] );
       ma:= X*ma*Xi;
       return [ ma[i][i] : i in [1..NumberOfColumns(ma)] ];

    end function;

    return func;


end function;


multrel1:= function( b )
  
   if &and[ IsOne(u) : u in b ] then
      return IdentityMatrix(Integers(),#b);
   end if;
   if &and[ IsZero(u) : u in b ] then  // strange case!!!! (happens if idpot is not the id mat!)
      return IdentityMatrix(Integers(),#b);
   end if;

   c:= [ ];
   for u in b do
       if not IsOne(u)  and not u in c then Append( ~c, u ); end if;
   end for;
   if #c eq 0 then Append( ~c, b[1] ); end if;

   //take 400 of them at a time...
   offset:= 400;
   k:= 1;
   c0:= [ ];
   while k+offset-1 le #c do
       d:= [ c[i] : i in [k..k+offset-1] ];
       H,m:=  MultiplicativeGroup(d);
       c0 cat:= [ m(H.j) : j in [1..Ngens(H)] ];
       k:= k+offset;
   end while;

   c0 cat:= [ c[i] : i in [k..#c] ];

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

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

end function;


multrels:= function( a, f ) // list of elements of the group alg, with the function...


   m:= [ f(u) : u in a ];
   n:= #m[1];

   Lam:= Lattice(multrel1( [ u[1] : u in m ] ));
   mat:= HermiteForm(Matrix( Basis(Lam) ));

   for i in [2..n] do
       sat:= true;
       for k in [1..NumberOfRows(mat)] do
           if not IsOne( &*[ m[j][i]^mat[k][j] : j in [1..NumberOfColumns(mat)] ] ) then
              sat:= false;
              break k;
           end if;
       end for;
       if not sat then
          Lam:= Lam meet Lattice(multrel1( [ u[i] : u in m ] ));
          mat:= HermiteForm(Matrix( Basis(Lam) ));
       end if;
   end for;

   return mat;

end function;

index:= function( small, big, f )

    mat:= multrels( small cat big, f );
    mat:= HermiteForm(mat);
    
    m:= #small;
    n:= #big;
    if NumberOfRows(mat) lt m then
       print "small is not contained in big";
       return 0;
    end if; 
    for i in [1..m] do
        if mat[i][i] ne 1 then
           print "generator ",small[i], " not contained in big group...";
        end if;
    end for;
    
    H:= [ ];
    for i in [1..NumberOfRows(mat)] do
        Append( ~H, [ mat[i][j] : j in [m+1..m+n] ] );
    end for;

    H:= HermiteForm( Matrix(H) );
    return H, mat;

end function;

normal_gens:= function( gens, f )

    if #gens eq 0 then return gens,[]; end if;

    m:= multrels( gens, f );
    S,P,Q:= SmithForm( m );

    k:= 1;
    while S[k][k] eq 1 do 
        k:= k+1; 
        if k gt NumberOfRows(S) then 
           return [], [];
        end if;
    end while;

    ord:= [ ];
    for i in [k..NumberOfRows(S)] do Append( ~ord, S[i][i] ); end for;
    ord cat:= [ 0 : i in [NumberOfRows(S)+1..NumberOfRows(Q)] ];

    Qi:= Q^-1;
    return [ &*[ gens[i]^Qi[j][i] : i in [1..#gens] ] : j in [k..NumberOfRows(Qi)] ], ord;

end function;



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;

//Input:  Matrix group G, matrix g in Generic(G)
//Output: boolean (g in G?) and, if true, then an SLP w in G`SLPGroup
//        such that g = Evaluate(w,UserGenerators(G));
//        If G`SLPGroup and G`UserGenerators are not defined, then
//        it will be set w.r.t. the generators of G.
Element2Word := function(G, g)
    if not Type(G) eq GrpMat or not g in Generic(G) then 
        error "input has to be matrix group G and elt g in Generic(G)"; 
    end if;
    if not assigned G`UserGenerators then
        G`UserGenerators := [G.i : i in [1..Ngens(G)]];
    end if;
    if not assigned G`SLPGroup or 
       not #UserGenerators(G) eq #UserGenerators(G`SLPGroup) then
        G`SLPGroup := SLPGroup(#UserGenerators(G));
    end if;
    
  //copy original gens
    ug    := [Generic(G)!u : u in UserGenerators(G)];
    Gc    := sub<Generic(G) | [G.i : i in [1..Ngens(G)]]>;
    T     := CompositionTree(Gc);
    f,w   := CompositionTreeElementToWord(Gc,g);
    if not f then
        error "element is not in group (?)";
        return false,_;
    end if;

  //write w.r.t. original user generators
    pos   := [Position(ug,Gc.i):i in [1..Ngens(Gc)]];
    map   := CompositionTreeNiceToUser(Gc);
    w     := Image(map,w);
    w     := Evaluate(w,[G`SLPGroup.i : i in pos]);

  //test
    assert g eq Evaluate(w, UserGenerators(G));
    
    return true,w;
end function;





//Input:  Matrix group G, matrix g in Generic(G)
//Output: boolean (g in G?) and, if true, then a list of integers k
//        such that 
//
//           g = &*[UserGenerators(G)[i]^k[i] : i in [1..#UserGenerators(G)];
//
//        If G`UserGenerators is not defined, then
//        it will be set w.r.t. the generators of G.
Element2WordAbelian := function(G, g)
    if not IsAbelian(G) then error "group has to be abelian"; end if;
    f, w := Element2Word(G,g);
    if not f then return false,_; end if;
    n    := #UserGenerators(G);
    ug   := UserGenerators(FreeAbelianGroup(n));
    w    := Evaluate(w,ug);
    exp  := Eltseq(w);

  //test
    assert g eq &*[UserGenerators(G)[i]^(exp[i]) : i in [1..n]];

    return true,exp;
end function;


normalgens_mats:= function( G, mats )

   // here mats is a list of commuting matrices over a finite field, lying in G = GL(d,F);
   // we compute normal generators.

   gg:= [ mats[1] ];
   H:= sub< G | [ mats[1] ] >;
   orders:= [  Order( mats[1] ) ];
   for k in [2..#mats] do
       // determine its smallest power in G...
       h:= mats[k];
       f:= 1;
       while not h in H do
          h *:= mats[k];
          f+:= 1;
       end while;
       if f gt 1 then
          len:= #gg+1;
          mt:= [ ];
          for i in [1..#gg] do
              v:= [ 0 : j in [1..len] ];
              v[i]:= orders[i];
              Append( ~mt, v );
          end for;

          bool, exp:= Element2WordAbelian(H,h);

          Append ( ~mt, exp cat [-f] );
          mt:= Matrix( Integers(), mt );
          S,P,Q:= SmithForm( mt );

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

          Append( ~orders, Order(mats[k]) );
          Append( ~gg, mats[k] );

          neword:= [ ];         
          for i in [1..NumberOfColumns(S)] do
              if S[i][i] ne 1 then
                 a:= One(G);
                 for j in [1..NumberOfColumns(Q)] do
                     s:= Qi[i][j] mod orders[j];
                     a *:= gg[j]^s;
                 end for;
                 Append( ~newgg, a );
                 Append( ~neword, S[i][i] );
              end if;
          end for;

          gg:= newgg;
          orders:= neword;
          H:= sub< G | gg >;
          H`UserGenerators:= gg;
       end if;
   end for;

   return gg, orders;

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, &+[ 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; 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;


ext_cyc:= function( A, e1, e2, I1, I2, u1, u2, ord1, ord2, z0, n0, mps, i0, ee  )

   // assume I2 is an order in a number field; here we assume that it is the 
   // cyclotomic field of order n0, and that z0 in I2 is a primitive n0-th root of unity.

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

    // Now O is a subset of J, Oi are subsets of Ii, Aei are their preimages in A...
    // (How about Ae1[i] = A!O1[i]?).

    //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] ];

    b:= b1 cat b2;

    F<w0>:= CyclotomicField( n0 );
    cycO:= RingOfIntegers( F );

    //make basis of I2 consisting of powers of z0...
    bas2:= [ One(I2) ] cat [ z0^k : k in [1..Dimension(I2)-1] ];
    V2:= KModule( Rationals(), Dimension(I2) );
    U2:= KModuleWithBasis( [ V2!Coordinates( I2, v ) : v in bas2 ] );

    phi:= function( a ) // a in I2, image in F...

         return F!Coordinates( U2, V2!Coordinates( I2, a ) );

    end function;

    // get the image of e2O \cap O (ie ker alpha_2) in cycO...

    I:= ideal< cycO | [ cycO!(phi(I2!u)) : u in b2 ] >;

    Q,map:= quo< cycO | I >;
    G,mG:= UnitGroup( Q );

    // now make maps Ok --> Q, for I2 this is easy...

    alpha2:= function( a ) //a in O2

         return (cycO!phi(a))@map;

    end function;

    M1:= KModule( Rationals(), Dimension(I1) );
    VO1:= KModuleWithBasis( [ M1!Coordinates( I1, u ) : u in O1 ] );

    alpha1:= function( a ) // a in O1

        cf:= Coordinates( VO1, M1!Coordinates( I1, a ) );
        v:= &+[ cf[i]*Ae1[i] : i in [1..#cf] ];
        v:= id*v;  // now in O...         

        return (cycO!(phi(I2!(e2*v))))@map;

    end function;

    g1:= [ alpha1(u)@@mG : u in u1 ];
    H1:= sub< G | g1 >;
    g2:= [ alpha2(u)@@mG : u in u2 ];
    H2:= sub< G | g2 >;

    H:= H1 meet H2;
    
    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 >;
    f1:= hom< Ab1 -> H1 | g1 >;
    f2:= hom< Ab2 -> H2 | g2 >;

    a1:= [ Ab1!u : u in SetToSequence(Generators( Kernel( f1 ) )) ];

    for g in Generators(H) do
        Append( ~a1, preim( f1, H1!g ) );  
    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 Generators(H) do
        Append( ~a2, preim( f2, H2!g ) );
    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 > --> H.

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

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

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

    for j in [1..Ngens(H)] do
        N:= Order( H.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:= [ ];
    for v in u1 do
        Append( ~inv1, &+[ mps[j][2]( (mps[j][1](v))^-1 ) : j in [1..i0-1] ] );
    end for;

    inv2:= [ ];
    for v in u2 do 
        Append( ~inv2, mps[i0][2]( (mps[i0][1](v))^-1 ) );
    end for;

    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;


intrinsic UnitsOfAbelianGroupRing( o::SeqEnum ) -> GrpPC, AlgGrp, SeqEnum, SeqEnum
{Generators of the unit group of an abelian group ring}


    G:= AbelianGroup( GrpPC, o );
    A,u,ord,I,e,f,projs,fields,prims,mps:= grpalg_data( G );

    id:= e[1]; J:= I[1]; v:= u[1]; ov:= ord[1]; 
    for i in [2..#e] do
        w0:= [ ];

        if Dimension( I[i] ) gt 1 and #w0 eq 0 then
           w0:= ext_cyc( A, id, e[i], J, I[i], v, u[i], ov, ord[i], prims[i], fields[i], mps, i, e );
        end if;
        
        if #w0 eq 0 then
           w0:= extendunits( A, id, e[i], J, I[i], v, u[i], ov, ord[i] );
        end if;

        if i lt #e then
           ip:= &+[ e[k] : k in [i+1..#e] ];
        else
           ip:= Zero(A);
        end if;

        FA:= FreeAbelianGroup( &+[ #ord[k] : k in [1..i] ] );
        mods:= Flat( [ ord[k] : k in [1..i] ]);
        U,map:= quo< FA | [ mods[k]*FA.k : k in [1..#mods] ] >;

        g1:= [ ];
        for w in w0 do
            vw:= [ ];
            for k in [1..i] do
                vw cat:= projs[k]( w );
            end for;
            Append( ~g1, (FA!vw)@map );
        end for;

        mat:= [ Eltseq(u) : u in g1 ];
        for k in [1..Ngens(U)] do
            if Order(U.k) gt 0 then
               cf:= [ 0 : j in [1..Ngens(U)] ];
               cf[k]:= Order(U.k);
               Append( ~mat, cf );
            end if;
        end for;

        ker:= Kernel( Matrix( Integers(), mat ) );
        rels:= [ ];
        for k in [1..Dimension(ker)] do
            cf:= [ (ker.k)[j] : j in [1..#w0] ];
            Append( ~rels, cf );
        end for;

        S,P,Q:= SmithForm( Matrix( Integers(), rels ) );

        k:= 1;
        while S[k][k] eq 1 do 
           k:= k+1; 
        end while;

        ov:= [ ];
        for j in [k..NumberOfRows(S)] do Append( ~ov, S[j][j] ); end for;
        ov cat:= [ 0 : j in [NumberOfRows(S)+1..NumberOfRows(Q)] ];

        Qi:= Q^-1;

        v:= [ ];
        id:= id+e[i];
        J:= ideal< A | [id] >; 
        invs:= [ inv( id, J, u ) : u in w0 ];
        for j in [k..NumberOfRows(Qi)] do
            
            cf:= Eltseq( Qi[j] );
            a:= id;
            for l in [1..#cf] do
                if cf[l] gt 0 then
                   a *:= w0[l]^cf[l];
                end if;
                if cf[l] lt 0 then
                   a *:= invs[l]^(-cf[l]);
                end if;
            end for;
            Append( ~v, a ); 
        end for;

    end for;

    return G, A, v, ov;

end intrinsic;


