#############################################################################
##
#W ispoly.gi               POLENTA package                     Bjoern Assmann
##
## Methods for testing if a rational matrix group is polycyclic
##
#Y 2005
##

#############################################################################
##
#F  POL_Logarithm( x )
##
## IN: x ..... unipotent matrix
##
##  OUT: log(x)
##
POL_Logarithm := function( x )
    local n,mat,matPow,log,i;
    n := Length( x[1] );
    mat := x - x^0;
    matPow := mat;
    log := mat;
    for i in [2..n-1] do
        matPow := matPow*mat;
        log := log + (-1)^(i-1)*(i)^-1*matPow;
    od;
    return log;
end;

#############################################################################
##
#F POL_Exponential( v )
##
## IN: v ..... matrix which is conjugated to an
##             upper or lower triangular matrix with 0's on the diagonal
##
## OUT: exp(x)
##
POL_Exponential := function( v )
    local n,exp, vPow, fac, i;
    n := Length( v[1] );
    #exp := IdentityMat(n) + v;
    exp := v^0 + v;
    vPow := v;
    #fac := 1;
    fac := v[1][1]^0;
    for i in [2..n-1] do
        vPow := vPow*v;
        fac := fac*i;
        exp := exp + (fac^-1)*vPow;
    od;
    return exp;
end;

#############################################################################
##
#F POL_CloseMatrixSpaceUnderLieBracket( V )
##
## IN:  V ....  vector space < M_nxn( Q )
##
## OUT: smallest Lie algebra containing V
##
POL_CloseMatrixSpaceUnderLieBracket := function( V )
    local basis,n,com,i,j,basis_ext, V_ext;

    basis := Basis( V );
    n := Length( basis );

    for i in [1..n] do
        for j in [i+1..n] do
            com := basis[i]*basis[j]-basis[j]*basis[i];
            if not com in V then
               basis_ext := POL_CopyVectorList( basis );
               Add( basis_ext, com );
               V_ext := VectorSpace( Rationals, basis_ext, "basis" );
               return POL_CloseMatrixSpaceUnderLieBracket( V_ext );
            fi;
        od;
    od;
    return V;
end;

#############################################################################
##
#F
##
## IN: unipo_gens ..... matrices that generate a unipotent matrix group
##
## OUT: The Lie algebra L( <unipo_gens> )
##
POL_LieAlgebra := function( unipo_gens )
    local logs,V,V_ext;

    # get logs
    logs := List( unipo_gens, POL_Logarithm );

    # compute Q-span
    V := VectorSpace( Rationals, logs );

    # close under Lie bracket
    V_ext := POL_CloseMatrixSpaceUnderLieBracket( V );

    return V_ext;
end;

#############################################################################
##
#F POL_CloseLieAlgebraUnderGrpAction( gensG, L )
##
## IN: gensG ... generators of rational matrix group G of degree n
##     L  ...... nilpotent Lie Algebra generated by matrices of degree n
##
## OUT: L^G, where G acts on L by l^g = log( exp(l)^g );
##
POL_CloseLieAlgebraUnderGrpAction := function( gensG, L )
    local basis,b,g,exp,u,l,basis_ext,V_ext,L_ext;
    basis := Basis( L );
    for b in basis do
        for g in gensG do
            exp := POL_Exponential( b );
            u := exp^g;
            l := POL_Logarithm( u );
            if not l in L then
               basis_ext := POL_CopyVectorList( basis );
               Add( basis_ext, l );
               V_ext := VectorSpace( Rationals, basis_ext, "basis" );
               L_ext := POL_CloseMatrixSpaceUnderLieBracket( V_ext );
               return POL_CloseLieAlgebraUnderGrpAction( gensG, L_ext );
            fi;
        od;
    od;
    return L;
end;

#############################################################################
##
#F
##
## IN:  L ......... Lie algebra
##      g ......... group element
##
## OUT: Induced action of g on the Lie algebra, i.e. the linear mapping
##      which maps l to log( exp(l)^g )
##
POL_InducedActionToLieAlgebra := function( g, L )
    local basis,n,inducedAction,i,exp,u,log,coeff;

    basis := Basis( L );
    n := Length( basis );
    inducedAction := [];

    for i in [1..n] do
        exp := POL_Exponential( basis[i] );
        u := exp^g;
        log := POL_Logarithm( u );
        coeff := Coefficients( basis, log );
        Add( inducedAction, coeff );
    od;
    return inducedAction;
end;

#############################################################################
##
#F POL_IsIntegerList( list )
##
##
POL_IsIntegerList := function( list )
    local z;
    for z in list do
        if not z in Integers then
            return false;
        fi;
    od;
    return true;
end;

#############################################################################
##
#F POL_IsIntegralActionOnLieAlgebra( gens, L )
##
## IN: gens ....... list of matrices
##     L ........... Lie algebra on which gens acts
##                   via l^g = Log( Exp(l)^g )
##
## OUT: returns true if for all g in gens, the minimal polynomial of
##      the induced action of  g,g^-1 to L is integral.
##      false otherwise.
##
POL_IsIntegralActionOnLieAlgebra := function( gens, L )
    local g,ind,pol,coeffs,constTerm,bool;
    Info( InfoPolenta, 3, "Testing whether action on Lie alg. is integral" );
    for g in gens do
        Info( InfoPolenta, 3, "Generator that will be induced:" );
        Info( InfoPolenta, 3, g );
        ind := POL_InducedActionToLieAlgebra( g, L );
        Info( InfoPolenta, 3, "Induced action to Lie algebra:" );
        Info(InfoPolenta, 3, ind );
        if not Trace( ind ) in Integers then
            Info( InfoPolenta, 3, "Trace not integral\n" );
            return false;
        fi;
        pol := CharacteristicPolynomial( Rationals, Rationals, ind );
        # pol := MinimalPolynomial( Rationals, ind );
        Info( InfoPolenta, 3, "Characteristic polynomial:" );
        Info(InfoPolenta, 3, pol );
        coeffs := CoefficientsOfLaurentPolynomial( pol );

        # test if pol_g in Z[x]
        if not POL_IsIntegerList( coeffs[1] ) then
            return false;
        fi;

        # test if pol_g^-1 in Z[X]
        constTerm := Value( pol, 0 );
        bool := (constTerm = 1 ) or ( constTerm = -1);
        if not bool then
           return false;
        fi;
    od;
    return true;
end;

#############################################################################
##
#F POL_IsIntegralActionOnLieAlgebra_Beals( gens, L )
##
## IN: gens ....... list of matrices
##     L ........... Lie algebra on which gens acts
##                   via l^g = Log( Exp(l)^g )
##
## OUT: returns true if for all g in gens, the minimal polynomial of
##      the induced action of  g,g^-1 to L is integral.
##      false otherwise.
##
POL_IsIntegralActionOnLieAlgebra_Beals := function( gens, L )
    local gens_ind, G_ind, lat;

    gens_ind :=  List( gens, x-> POL_InducedActionToLieAlgebra( x, L ));
    G_ind := GroupByGenerators( gens_ind );

    lat := InvariantLattice( G_ind );
    if lat = fail then
        return false;
    else
        return true;
    fi;

end;

#############################################################################
##
#F
##
## IN: gensU_p ...... normal subgroup generators for U_p
##     gensG ........ the generators of the parent group G
##     gensGU ......  representatives of pc sequence of GU
##
POL_IsFinitelgeneratedU_p := function( gensU_p, gensG, gensGU )
    local L,L_ext,isIntegral,gensU_p_mod,i,method;

    # catch trivial case G/U = 1
    if Length( gensGU ) = 0 then
       return true;
    fi;

    # catch trivial case U_p = 1
    gensU_p_mod := [];
    for i in [1..Length( gensU_p )] do
        if not gensU_p[i]=gensU_p[i]^0 then
            Add( gensU_p_mod,gensU_p[i] );
        fi;
    od;
    if Length( gensU_p_mod ) = 0 then
       return true;
    fi;

    # get Lie algebra corresponding to <gensU_p>
    L := POL_LieAlgebra( gensU_p_mod );
    Info( InfoPolenta, 3, "Dimension L: ", Dimension( L ) );

    # close it under action of G
    L_ext := POL_CloseLieAlgebraUnderGrpAction( gensG, L );
    Info( InfoPolenta, 3, "Dimension L_ext: ", Dimension( L_ext ) );
    Info( InfoPolenta, 3, "Basis of L_ext: ", Basis( L ) );

    # test whether the action of gensGU on the Lie algebra is integral
    method := "beals";
    if method = "char" then
        isIntegral := POL_IsIntegralActionOnLieAlgebra_Beals( gensG, L_ext );
    else
        isIntegral := POL_IsIntegralActionOnLieAlgebra( gensGU, L_ext );
    fi;

    return isIntegral;
end;

#############################################################################
##
#E