###########################################################
###########################################################
##
##  2009-2010, 2010-2011 UGA VIGRE Algebra Group
##
##  Last modified: 30 Sept 2011
##   
##  This file contains GAP code to perform calculations necessary
##  for the following papers by the UGA VIGRE Algebra Group:
##
##  [1] First cohomology for finite groups of Lie type:
##      Simple modules with small dominant weights
##  [2] Second cohomology for finite groups of Lie type
## 
###########################################################
###########################################################


###########################################################
###########################################################
##
##  om(j,n)
##
##  Constructs a vector of length n having 1 in position j
##  and zeros elsewhere. Using for constructing fundamental
##  dominant weights.
##
###########################################################
###########################################################

om := function(j,n)
  local out,k;
  out := [];
  for k in [1..n] do
    if k = j then
      Add(out,1);
    else
      Add(out,0);
    fi;
  od;
  return out;
end;


###########################################################
###########################################################
##
##  OneVector(n)
##
##  Constructs a vector of length n having all 1s as entries.
##  Useful for constructing the Weyl weight of a root system.
##
###########################################################
###########################################################

OneVector:=function(n)
  local k, out;
  out:=List([1..n],k->1);
  return out;
end;


###########################################################
###########################################################
##
##  ConstructRootSystemData(type,rank)
##
##  This function produces a GAP record containing various pieces of
##  root system data. It takes as input the Lie type, e.g., "B", and
##  the rank of the desired root system. The output is a record containing
##  the root system, Weyl group, Cartan matrix, inverse of the Cartan
##  matrix, the Weyl weight rho, the rank, and a list of the positive
##  roots (in the weight basis). Within the record, these items are
##  accessed via the following list of names:
##
##  rootsystem, weylgroup, cartanmat, invcartanmat, rho, rank, posroots
##
##  example: ConstructRootSystemData("E",6);
##
###########################################################
###########################################################

ConstructRootSystemData := function(type,n)
  local tmprho, liealgebra, rootsystemrec, tmprootsystem, tmpweylgroup,
        tmpcartanmat, tmpinvcartanmat, tmpposroots, rho, rootystem,
        weylgroup, cartanmat, invcartanmat, posroots;

  tmprho:=OneVector(n);

  liealgebra:=SimpleLieAlgebra(type,n,Rationals);

  tmprootsystem:=RootSystem(liealgebra);
  tmpweylgroup:=WeylGroup(tmprootsystem);
  tmpcartanmat:=CartanMatrix(tmprootsystem);
  tmpinvcartanmat:=tmpcartanmat^(-1);
  tmpposroots:=PositiveRoots(tmprootsystem);

  rootsystemrec:=rec(
    rootsystem:=tmprootsystem,
    weylgroup:=tmpweylgroup,
    cartanmat:=tmpcartanmat,
    invcartanmat:=tmpinvcartanmat,
    rho:=tmprho,
    rank:=n,
    posroots:=tmpposroots,
  );
  
  return rootsystemrec;
end;


###########################################################
###########################################################
##
##  TfqInvariants(weight)
##
##  This function checks for which prime powers q=p^r a given weight is
##  invariant for the finite group T(Fq), i.e., checks for which primes
##  powers q=p^r the weight is equal to (q-1) times a weight. The output
##  is the list of prime powers satisfying this condition.
##
##  example: TfqInvariants([6,14]);
##
###########################################################
###########################################################

TfqInvariants:=function(weight)
  local gcd,factors,i,powers;
  powers:=[];
  gcd:=Gcd(weight);
  if gcd=0 then return [1]; fi;
  factors:=DivisorsInt(gcd);
  for i in [1..Length(factors)] do
    if IsPrimePowerInt(factors[i]+1) then
      Add(powers,factors[i]+1);
    fi;
  od;
  return powers;
end;


###########################################################
###########################################################
##
##  TrInvariants(weight)
##
##  This function checks for which prime powers p^r a given weight is
##  invariant for the Frobenius kernel Tr, i.e., checks for which prime
##  powers q=p^r the weight is equal to p^r times a weight. The output
##  is the list of prime powers satisfying this condition.
##
##  example: TrInvariants([15,35]);
##
###########################################################
###########################################################

TrInvariants:=function(weight)
  local gcd,factors,i,powers;
  powers:=[];
  gcd:=Gcd(weight);
  if gcd=0 then return [1]; fi;
  factors:=DivisorsInt(gcd);
  for i in [1..Length(factors)] do
    if IsPrimePowerInt(factors[i]) then
      Add(powers,factors[i]);
    fi;
  od;
  return powers;
end;


###########################################################
###########################################################
##
##  H2uTfqInvariants(type,rank,lambda)
##
##  This function is designed to check for which prime powers q=p^r the Lie
##  algebra second cohomology group H^2(u,L(lambda)) could have nonzero
##  T(Fq) invariant weights. This function is used for the proof of
##  [2, Lemma 4.2.1], to show for the exceptional Lie types that the only
##  T(Fq)-invariant weight in H^2(u,L(lambda)) is the zero weight.
##
##  The function computes all sums alpha + beta + nu, with alpha and beta
##  ranging over all pairs of nonproportional roots, and nu ranging over
##  all possible dominant weights in L(lambda). The function checks whether
##  each sum is a T(Fq)-invariant weight. The output is the list of prime
##  powers satisfying this condition.
##
##  example: HwuTfqInvariants("B",4,[1,0,0,0]);
##
###########################################################
###########################################################

H2uTfqInvariants:=function(type,rank,lambda)

  local rootsysdata,phi,alpha,beta,nu,mu,powerList;
	
  rootsysdata:=ConstructRootSystemData(type,rank);
  phi:=Concatenation(rootsysdata.posroots,-rootsysdata.posroots);

  powerList:=[];

  # For each dominant weight nu in L(lambda), and each pair of
  # non-proportional roots alpha and beta, compute mu = alpha + beta + nu
  # and determine the list of prime powers p^r for which mu is
  # (p^r-1) times a weight.

  for nu in DominantWeights(rootsysdata.rootsystem,lambda)[1] do
    for alpha in phi do
      for beta in phi do
        if alpha <> beta and alpha <> -beta then
   	  mu:=alpha+beta+nu;
  	  powerList:=Union(powerList,TfqInvariants(mu));
    	fi;
      od;
    od;
  od;
  return powerList;
end;;


###########################################################
###########################################################
##
##  wdotwtpluswtTrInvariants(type,rank,lambda)
##
##  This function computes (w dot lambda) + lambda for all w in the Weyl
##  group. For each Weyl group element, it checks the prime powers p^r
##  for which the sum could be in p^r times the root lattice. The
##  output is the list of prime powers satisfying this condition.
##
##  This function is used for the proof of [2, Theorem 1.2.5], to
##  rule out weights of a certain form in the Lie algebra
##  cohomology group H^1(u,L(lambda)).
##
##  example: wdotwtpluswtTrInvariants("G",2,[0,1]);
##
###########################################################
###########################################################

wdotwtpluswtTrInvariants:=function(type,rank,lambda)

  local powers, rootsysdata, W, rho, Cinv, wlambdaplusrho,
        dotaction, basisChange, j, iteratorcounter;

  powers:=[];
  iteratorcounter:=0;

  rootsysdata:=ConstructRootSystemData(type,rank);
  W:=rootsysdata.weylgroup;
  rho:=rootsysdata.rho;
  Cinv:=rootsysdata.invcartanmat;
  
  for wlambdaplusrho in WeylOrbitIterator(W,lambda+rho) do

    ## For large Weyl groups, print out timing
    ## information every half-million elements.
    iteratorcounter:= iteratorcounter+1;
  	 if IsInt(iteratorcounter/500000) then
    	    Print("element number ",iteratorcounter,"\n");
         fi;

    # compute (w dot lambda) + lambda
    dotaction:= wlambdaplusrho - rho + lambda;

    # convert to root basis
    basisChange:= dotaction*Cinv;

    # check whether basisChange is in the root lattice
    if ForAll(basisChange,j->IsInt(j)) then

      # check prime powers for which basisChange is in root lattice
      powers:=Union(powers,TrInvariants(basisChange));

    fi;

  od;

  return powers;
end;


###########################################################
###########################################################
##
##  HigherWeights(rootsystemdata,tau,listofweights)
##
##  Output is the subset of listofweights that are greater
##  than the given weight tau. The argument rootsystemdata
##  should be a GAP record constructed as the output of the
##  function ConstructRootSystemData defined above.
##
##  example:
##  rootsystemdata:=ConstructRootSystemData("G",2);
##  g2:= [om(0,2),om(1,2),om(2,2)];
##  HigherWeights(rootsystemdata,[1,0],g2);
##
###########################################################
###########################################################

HigherWeights := function(rootsystemdata,tau,listofweights)
  local outList, lambda, tmp;
  outList := [];
  for lambda in listofweights do
    tmp:=(lambda - tau)*rootsystemdata.invcartanmat;
    if ForAll(tmp,j->IsPosInt(j+1)) and lambda <> tau then
       Add(outList,lambda);
    fi;
  od;
  return outList;
end;


###########################################################
###########################################################
##
##  PrimesForLinkage(rootsystemdata,tau,lambda)
##
##  Given a GAP record of root system data, as constructed by
##  the function ConstructRootSystemData, and given a pair of
##  weights tau and lambda, the function gives as output the
##  list of primes for which the two weights are linked under
##  the dot action of the affine Weyl group Wp.
##
##  For large Weyl groups, this function prints timing
##  information every half-million elements.
##
##  example:
##  rootsystemdata:=ConstructRootSystemData("G",2);
##  PrimesForLinkage(rootsystemdata,[1,0],[0,1]);
##
###########################################################
###########################################################

PrimesForLinkage:=function(rootsystemdata,tau,lambda)
  local primes, iteratorcounter, taurho, wtaurho, Rt, dtau,
        dtauminuslambda, basisChange, divisor, tmp;

  primes:=[];
  Print("Checking for linkage between ",tau," and ",lambda,"..\n");

  iteratorcounter:=0;
  taurho := tau + rootsystemdata.rho;

  ## Begin iterating over the Weyl group.
  for wtaurho in WeylOrbitIterator(rootsystemdata.weylgroup,taurho) do

    ## For large Weyl groups, print out timing
    ## information every half-million elements.
    iteratorcounter:= iteratorcounter+1;
  	 if IsInt(iteratorcounter/500000) then
    	    Print("element number ",iteratorcounter,"\n",
                  "tau = ",tau," lambda=",lambda,"\n");

            # Uncomment the following lines to print runtimes.
    	    # Rt:=Runtimes();
    	    # Print(Rt,"\n");
  	 fi;

    ## compute dot action
    dtau := wtaurho - rootsystemdata.rho;

    ## compute dtau - lambda
    dtauminuslambda:= dtau-lambda;

    ## convert dtauminuslambda to root basis
    basisChange:= dtauminuslambda*rootsystemdata.invcartanmat;

    ## check whether basisChange is in the root lattice
    if ForAll(basisChange,j->IsInt(j)) then

      ## if basisChange is in the root lattice, compute
      ## the primes for which basisChange is in p times
      ## the root lattice
      divisor:=Gcd(basisChange);
      tmp:=Set(FactorsInt(divisor));
      primes:=Union(primes,tmp);

    fi;

  od;

  Print(tau," and ",lambda," are linked for the primes ",primes,"\n");
  return primes;
end;


###########################################################
###########################################################
##
##  CheckLinks(type,rank,listofweights,"alllinks","logfile")
##
##  This function checks for what primes p the dominant weights
##  in the list listofweights are linked to zero under the dot
##  action of the affine Weyl group Wp. Here Wp is the affine
##   Weyl group at the prime p for the given Lie type and rank.
##
##  The fourth arguments is optional. If the fourth argument
##  "alllinks" is present, the function will compute the linkage
##  relations among all pairs of weights in listofweights.
##
##  The fifth argument is also optional. If present, the function
##  will log the output of the function to "logfile".
##
##  If only four arguments are present, and if the fourth argument
##  is not "alllinks", the function will assume that the fourth
##  argument is the name for the log file.
##
##  This function is used to check when the weights listed in
##  [1, Table 7.1] are linked to the zero weight. It is also used
##  for the proof of [2, Theorem 1.2.5], to check when the highest
##  long and short roots are linked.
##
##  examples:
##  g2:= [om(0,2),om(1,2),om(2,2)];
##  CheckLinks("G",2,g2,"alllinks","G2data.txt");
##  CheckLinks("G",2,g2,"alllinks");
##  CheckLinks("G",2,g2,"G2data.txt");
##  CheckLinks("G",2,g2);
##
###########################################################
###########################################################

CheckLinks:=function(arg)
  local type, rank, zeroweight, tauList, onlyzero, rootsystemdata,
        log, logging, outList, tau, lambdalist, lambda, primes;

  ## pull data from the list of inputs
  type:=arg[1];
  rank:=arg[2];
  tauList:=arg[3];
  zeroweight:=om(0,rank);

  ## Construct root system data
  rootsystemdata:=ConstructRootSystemData(type,rank);

  ## By default, logging to a file is turned off.
  ## If the optional logfile argument is supplied,
  ## then logging is turned back on.

  logging:="off";

  ## By default, the function only checks linkage
  ## relations with the zero weight.
  onlyzero:=true;

  if Length(arg) > 3 and arg[4] = "alllinks" then
     onlyzero:=false;
  fi;

  if Length(arg) > 3 and arg[4] <> "alllinks" then
     log:=arg[4];
     logging:="on";
  fi;

  if Length(arg) > 4 then
    log:=arg[5];
    logging:="on";
  fi;

  ## If only checking for linkage to zero, add the zero
  ## weight to tauList, if not already present.
  if onlyzero then
     tauList:=Union([zeroweight],arg[3]);
  fi;

  ## Print the status of logging.
  if logging = "on" then
     Print("Logging to file ",log,"\n");
     LogTo(log);
  fi;

  outList:=[];

  ## Loop through values of tau in tauList (i.e., in listofweights).
  for tau in tauList do

    ## If we are only checking linkage relations with zero, then
    ## skip over all nonzero values of tau.
    if onlyzero and tau <> om(0,rank) then
       continue;
    fi;

    ## Determine all weights in tauList that are greater than tau.
    ## In particular, if tau = 0, then tau is only linked to
    ## a dominant weight lambda if lambda > 0 in the root ordering.
    lambdalist:=HigherWeights(rootsystemdata,tau,tauList);

    if lambdalist = [] then
       Print(tau," is maximal..skipping this value of tau\n\n");
       continue;
    fi;    

    for lambda in lambdalist do
      primes:=PrimesForLinkage(rootsystemdata,tau,lambda);
      Add(outList,[primes,tau,lambda]);
    od;

  od;

   # stop logging
   if logging = "on" then
      Print("\n",outList,"\n");
      Print("Stopping log..\n");
      LogTo();
   fi;

  return outList;
end;


###########################################################
###########################################################
##
##  these are the tau such that
##  tau \leq omega_j for some j
##
##  WARNING! GAP does not use the Bourbaki ordering for the
##  fundamental weights in type F4. To get from GAP's
##  ordering to the Bourbaki ordering, apply the 3-cycle (142).
##
###########################################################
###########################################################

possibleTauRec := rec(
  e8:= [om(4,8),
          om(1,8)+om(6,8),
	  2*om(1,8)+om(8,8),om(2,8)+om(7,8),
          om(3,8)+om(8,8),2*om(7,8),
	  om(1,8)+om(2,8),om(6,8)+om(8,8),
	  om(5,8),om(1,8)+2*om(8,8),
	  om(1,8)+om(7,8),
	  2*om(1,8),om(2,8)+om(8,8),3*om(8,8),
	  om(3,8),om(7,8)+om(8,8),
	  om(6,8),
	  om(1,8)+om(8,8),
	  om(2,8),2*om(8,8),
	  om(7,8),
	  om(1,8),
	  om(8,8),
	  om(0,8)],
  e7:= [om(0,7),om(1,7),om(6,7),om(3,7),2*om(7,7),om(2,7)+om(7,7),
          2*om(1,7),om(1,7)+om(6,7),om(4,7),
	  om(7,7),om(2,7),om(1,7)+om(7,7),om(5,7)],
  e6:= [om(1,6),om(5,6),
          om(0,6),om(2,6),om(1,6)+om(6,6),om(4,6),
	  om(6,6),om(3,6)],
  f4:= [om(4,4),om(1,4)+om(2,4),2*om(1,4),om(3,4),om(2,4),om(1,4),om(0,4),],
  g2:= [om(0,2),om(1,2),om(2,2)],
  c4:= [om(0,4),om(1,4),om(2,4),om(3,4),om(4,4)],
  c6:= [om(0,6),om(1,6),om(2,6),om(3,6),om(4,6),om(5,6),om(6,6)]
);