############################################################################
##
##  GramDets.g                                               (Frank Luebeck)
##
##  Implementation of the recursive formula in Thm 4.11 of
##          Hebin Rui and Mei Si, Discriminants of Brauer algebras
##
##  To try this out start GAP, read this file and use the
##  function ComputeAllGramDets:
##
##  gap> Read("GramDets.g");
##  gap> for n in [1..20] do
##  gap>   ComputeAllGramDets(n, 2);
##  gap> od;
##
##  and so on.



##  The results of these computations would be too big, even for a computer,
##  if they are not computed in factorized form. But this is easy to do because
##  the recursive results are all multiplied.
##
##  The leading coefficient of the resulting polynomials is always a product of
##  small primes p < 2n. We store a product
##     s := (-1)^a * p1^a1 * p2^a2 * ... * pk^ak,   p1 < p2 < ... < pk primes
##  as a list
##     [a,...,ak]
##  which has length pk, first entry a, entry pi is ai for i=1..k, and the other
##  entries are zero. Note that then a power s^r with r some integer is
##  represented by
##    r * [a,...,ak].
##
##  Similarly, we handle the polynomial part. We are all the time only
##  multiplying factors of the form (delta + c) where c is an integer in
##  the range -2n < c < 2n. We store a polynomial
##    t := (delta+c1)^a1 * (delta+c2)^a2 * ... * (delta+ck)^ak
##  as a list of length 4m with 2m bigger than the biggest absolute value of any
##  ci, i=1..k∴
##    [0,...,ai, ...]
##  Here an entry ai is stored in position ci+2m, the other entries are 0.
##  Note again, that multiplying two such expressions, written with respect to
##  the same m, are easily multiplied by adding the lists component-wise.
##  Powering is again multiplying the list with the exponent.
##
##
##  The main function is detGfla(f, lambda) which computes the determinant of
##  the Gram matrix for (f, lambda). If the result has the form
##    c * pol, c a number and pol a monic polynomial in delta
##  then detGfla returns the result in the form
##    [clist, pollist]
##  where clist is a factored number and pollist a factored polynomial as
##  described above.
##
##  During the recursion is happens that some results from smaller n are
##  needed again and again. Therefore, the results which are computed
##  should be stored and used again when they are needed next time.
##  This is done with a list CACHEdetGfla of the form
##       [flambdas, results]
##  where 'flambdas' is a list of pairs [f,lambda], 'results' contains the
##  result of detGfla(f, lambda) in the same position as [f,lambda] is stored.
##

############################################################################
##  Here is first a function which multiplies two lists representing
##  factored numbers as described above. And second a function for
##  multiplying such a list with an integer (which is then factored).
# multiply factored numbers, store exponent of prime p in position p
# (only small primes < 2n occur), and sign by an exponent in position 1
multFactoredInts := function(a, b)
  local i;
  # if a and b have different lengths we first append some zeros to the
  # shorter list - then we can easily add component-wise
  if Length(a) < Length(b) then
    a := ShallowCopy(a);
    for i in [Length(a)+1..Length(b)] do
      Add(a,0);
    od;
  fi;
  if Length(b) < Length(a) then
    b := ShallowCopy(b);
    for i in [Length(b)+1..Length(a)] do
      Add(b,0);
    od;
  fi;
  return a+b;
end;
# multiply with not yet factored number r
multFactoredWithInt := function(a, r)
  local coll, i, x;
  if r = 0 then
    Error("factor cannot be zero");
  fi;
  if Length(a) > 0 then
    a := ShallowCopy(a);
  else
    a := [0];
  fi;
  if r < 0 then
    a[1] := a[1]+1;
    r := -r;
  fi;
  if r > 1 then
    # this computesѕ factorization as list of pairs [pi, ai]
    coll := Collected(Factors(r));
    # initialize with zeros (first entry of last entry in coll is largest
    # prime factor)
    for i in [Length(a)+1..coll[Length(coll)][1]] do
      a[i] := 0;
    od;
    for x in coll do
      a[x[1]] := a[x[1]] + x[2];
    od;
  fi;
  return a;
end;

############################################################################
##  Here is a function which prints a result of detGfla in ore human readable
##  form. (It first builds it in a string which is then be printed.)
##
stringFactoredDet := function(res)
  local str, nn, a, i;
  str := "";
  if res[1] <> 0*res[1] then
    if res[1][1] mod 2 <> 0 then
      Add(str, '-');
    fi;
    for i in [2..Length(res[1])] do
      if res[1][i] <> 0 then
        Append(str, String(i));
        if res[1][i] <> 1 then
          Add(str, '^');
          Append(str, String(res[1][i]));
        fi;
        Append(str, " * ");
      fi;
    od;
  fi;
  nn := Length(res[2])/4;
  for i in [1..4*nn] do
    if res[2][i] <> 0 then
      a := i-2*nn;
      if a = 0 then
        Add(str, 'd');
      elif a > 0 then
        Append(str, "(d+");
        Append(str, String(a));
        Add(str, ')');
      else
        Append(str, "(d");
        Append(str, String(a));
        Add(str, ')');
      fi;
      if res[2][i] <> 1 then
        Add(str, '^');
        Append(str, String(res[2][i]));
      fi;
      Append(str," * ");
    fi;
  od;
  if Length(str) = 0 then
    str := "1";
  elif str[Length(str)] = '-' then
    Add(str, '1');
  elif Length(str) > 2 and str[Length(str)-1] = '*' then
    Unbind(str[Length(str)]);
    Unbind(str[Length(str)]);
    Unbind(str[Length(str)]);
  fi;
  return str;
end;

printFactoredDet := function(res)
  Print(stringFactoredDet(res),"\n");
end;

############################################################################
##  This computes rank(Delta(f, lambda)). GAP can compute lambda' with
##  the command AssociatedPartition. In the formula before 2.8 we use
##  (2f)!  instead of  2f!
rankDelta := function(f, la)
  local lap, prodhook, i, j, n;
  # lambda'
  if Length(la)>0 then
    lap := AssociatedPartition(la);
  fi;
  n := Sum(la)+2*f;
  prodhook := 1;
  for i in [1..Length(la)] do
    for j in [1..la[i]] do
      prodhook := prodhook * (la[i]+lap[j]-i-j+1);
    od;
  od;
  return Factorial(n) * Product([1,3..2*f-1])/prodhook/Factorial(2*f);
end;

############################################################################
##  R(lambda). This returns the list of removable nodes of a partition lambda,
##  written as pairs [i,j] of coordinates in the Young diagram.
RemovableNodes := function(la)
  local res, i;
  res := [];
  for i in [1..Length(la)-1] do
    if la[i+1] < la[i] then
      Add(res, [i,la[i]]);
    fi;
  od;
  if Length(la) > 0 then
    Add(res, [Length(la),la[Length(la)]]);
  fi;
  return res;
end;
# A(lambda). And here the addable nodes.
AddableNodes := function(la)
  local res, i;
  if Length(la) = 0 then
    return [[1,1]];
  fi;
  res := [[1,la[1]+1]];
  for i in [2..Length(la)] do
    if la[i-1] > la[i] then
      Add(res, [i,la[i]+1]);
    fi;
  od;
  Add(res, [Length(la)+1,1]);
  return res;
end;


############################################################################
##  Here are functions for computing the expressions for \gamma_{\lambda/\mu}
##  in the cases of Propositions 4.5, 4.8 and 4.9.
# prop 4.5
gamma1 := function(la, p)
  local res, q;
  # this is -1 in factored form
  res := [1];
  for q in AddableNodes(la) do
    # the if statement checks if an addable node is in A(lambda)^{<p}
    if q[1] > p[1] then
      res := multFactoredWithInt(res, -p[2]+p[1]+q[2]-q[1]);
    fi;
  od;
  for q in RemovableNodes(la) do
    if q[1] > p[1] then
      # must divide here (inverse of factored integer is described by
      # negative list
      res := - multFactoredWithInt(-res, -p[2]+p[1]+q[2]-q[1]);
    fi;
  od;
  return res;
end;

# partial product from 4.8 and 4.9 over AR(la)
gammaAR := function(la, p, n)
  local resc, resp, ind, q;
  # scalar and polynomial part of result
  resc := [0];
  # we store a factor (delta + a)^e by writing e in resp[a+2n]
  resp := 0*[1..4*n];
  # first big product
  for q in AddableNodes(la) do
    if q <> p and q[1] <= p[1] then
      resc := - multFactoredWithInt(-resc, -p[1]+p[2]+q[1]-q[2]);
      ind := -1-p[1]+p[2]-q[1]+q[2]+2*n;
      resp[ind] := resp[ind]+1;
    fi;
  od;
  for q in RemovableNodes(la) do
    if q[1] <= p[1] then
      resc := multFactoredWithInt(resc, -p[1]+p[2]+q[1]-q[2]);
      ind := -1+p[2]-p[1]-q[1]+q[2]+2*n;
      resp[ind] := resp[ind]-1;
    fi;
  od;
  return [resc, resp];
end;

# prop 4.8
gamma2 := function(la, p, n)
  local res, resp, ind;
  res := gammaAR(la, p, n);
  resp := res[2];
  # mu_k = 1 or not
  if p[2] = 1 then
    ind := -2*p[1]+2+2*n;
    resp[ind] := resp[ind]+1;
  else
    ind := p[2]-1-2*p[1]+2*n;
    resp[ind] := resp[ind]+1;
    ind := 2*p[2]-2*p[1]+2*n;
    resp[ind] := resp[ind]+1;
  fi;
  return res;
end;

# prop 4.9
gamma3 := function(la, p, n)
  local res, resc, resp, mu, ind, q;
  res := gammaAR(la, p, n);
  resc := res[1];
  resp := res[2];
  mu := ShallowCopy(la);
  mu[p[1]] := p[2];
  for q in AddableNodes(mu) do
    if q[1] > p[1] then
      resc[1] := resc[1]+1;
      ind := -(1+p[1]-p[2]+q[1]-q[2])+2*n;
      resp[ind] := resp[ind]+1;
    fi;
  od;
  for q in RemovableNodes(la) do
    if q[1] > p[1] then
      resc[1] := resc[1]+1;
      ind := -(1+p[1]-p[2]+q[1]-q[2])+2*n;
      resp[ind] := resp[ind]-1;
    fi;
  od;
  resc[1] := resc[1]+1;
  ind := 2*p[2]-2*p[1]+2*n;
  resp[ind] := resp[ind]+1;
  if not [p[1],p[2]-1] in RemovableNodes(la) then
    ind := -2+2*p[2]-2*p[1]+2*n;
    resp[ind] := resp[ind]-1;
  fi;
  return [resc, resp];
end;

# a cache for the  det G_{f,la}
CACHEdetGfla := [[],[]];
# function for adjusting the length of second (polynomial part) entry in result
adjustedCACHEdetGfla := function(res, nn)
  local m, v, a, i;
  m := Length(res[2])/4;
  if m = nn then
    return res;
  fi;
  v := 0*[1..4*nn];
  for i in [1..4*m] do
    if res[2][i] <> 0 then
      a := i - 2 * m + 2 * nn;
      if a < 1 or a > 4*nn then
        Error("Cannot adjust length of cache result.");
      fi;
      v[a] := res[2][i];
    fi;
  od;
  return [res[1],v];
end;


############################################################################
##  The main function implements the recursion formula for det G_{f,lambda}
##  from 4.11 using gamma1, gamma2 and gamma3 above. The optional third argument
##  says for which 'm' the factorized polynomial in the result schoud be
##  written.
##
# arguments f, la[, nn]
detGfla := function(arg)
  local f, la, n, nn, res, mu, fac1, fac2, exp, p, pos;
  f := arg[1];
  la := arg[2];
  n := 2*f + Sum(la);
  if Length(arg)=2 then
    nn := n;
  else
    # in recursion we may have nn > n
    nn := arg[3];
  fi;
  # first check cache, if result known just return it
  pos := PositionSorted(CACHEdetGfla[1], [f,la]);
  if IsBound(CACHEdetGfla[1][pos]) and CACHEdetGfla[1][pos] = [f,la] then
    return adjustedCACHEdetGfla(CACHEdetGfla[2][pos], nn);
  fi;

  # ok, not yet in cache, we must compute it
  res := [[0], 0*[1..4*nn]];
  if n = 1 then
    # this is the start of the recursion
    return res;
  fi;
  for p in RemovableNodes(la) do
    mu := ShallowCopy(la);
    if p[2] = 1 then
      Unbind(mu[Length(mu)]);
    else
      mu[p[1]] := p[2]-1;
    fi;
    fac1 := detGfla(f,mu,nn);
    fac2 := gamma1(la,p);
    exp := rankDelta(f,mu);
    res[1] := multFactoredInts(multFactoredInts(res[1], exp*fac2), fac1[1]);
    res[2] := res[2]+fac1[2];
  od;
  if f > 0 then
    for p in AddableNodes(la) do
      mu := ShallowCopy(la);
      mu[p[1]] := p[2];
      fac1 := detGfla(f-1,mu,nn);
      if p[1] >= Length(la) then
        fac2 := gamma2(la,p,nn);
      else
        fac2 := gamma3(la,p,nn);
      fi;
      exp := rankDelta(f-1,mu);
      res[1] := multFactoredInts(multFactoredInts(res[1],fac1[1]), exp*fac2[1]);
      res[2] := res[2]+fac1[2]+exp*fac2[2];
    od;
  fi;
  # before returning the result we store it in cache
  # (here the 'Immutable' and 'MakeImmutable' are not necessary, but
  # they make the program faster for some technical reasons)
  pos := PositionSorted(CACHEdetGfla[1], [f,la]);
  MakeImmutable(res);
  Add(CACHEdetGfla[1], Immutable([f,la]), pos);
  Add(CACHEdetGfla[2],res, pos);
  return res;
end;


############################################################################
##  Finally a function which computes all cases for fixed n, the second
##  argument 'printlevel' can be 0,1, or 2 for printing nothing, only [f.la],
##  or fl and the result, respectively.
ComputeAllGramDets := function(n, printlevel)
  local f, pp, res, la;
  f := 0;
  while 2*f <= n do
    pp := Partitions(n-2*f);
    for la in pp do
      if printlevel > 0 then
        Print([f, la], ": \c");
      fi;
      res := detGfla(f,la);
      if printlevel > 1 then
        Print(stringFactoredDet(res), "\n");
      fi;
    od;
    f := f + 1;
  od;
end;