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;

gamma1 := function(la, p, n)
  local res, resc, resq, q, i,resp1,resp2;
  res :=[0*[1..n],0*[1..3], 0*[1..4*n], 0*[1..4*n] ];
  resc :=res[1];
  resq :=res[2];
  resp1 := res[3];
  resp2 := res[4];
  resq := resq+[p[1]*2, 0, 0];
  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
       resq[1] := resq[1]+q[2]*2-q[1]*2;
      if p[2]-p[1]-q[2]+q[1] > 1 then
         resc[p[2]-p[1]-q[2]+q[1]] := resc[p[2]-p[1]-q[2]+q[1]]+1;
      fi;
    fi;
  od;
  for q in RemovableNodes(la) do
    if q[1] > p[1] then
      resq[1]:= resq[1]+q[1]*2-q[2]*2;
        if p[2]-p[1]-q[2]+q[1] > 1 then
         resc[p[2]-p[1]-q[2]+q[1]] :=resc[p[2]-p[1]-q[2]+q[1]]-1;
      fi;
    fi;
  od;
  return res;
end;

gammaAR := function(la, p, n)
  local resc, resp1, resp2, resq, ind, q;
  resc := 0*[1..n];
  resp1 := 0*[1..4*n];
  resp2 := 0*[1..4*n];
  resq :=0*[1..3];
  for q in AddableNodes(la) do
    if q <> p and q[1] < p[1] then
      resq[1]:= resq[1]+2*q[1]-2*q[2];
      ind := -p[1]+p[2]-q[1]+q[2]+2*n;
      resp1[ind] := resp1[ind]+1;
      resp2[ind] := resp2[ind]+1;
         if p[1]-p[2]-q[1]+q[2] > 1 then
            resc[p[1]-p[2]-q[1]+q[2]] := resc[p[1]-p[2]-q[1]+q[2]]-1;
         fi;
     fi;
  od;
  for q in RemovableNodes(la) do
    if q[1] <= p[1] then
      resq[1]:= resq[1]+2*q[2]-2*q[1];
        if p[1]-p[2]+q[2]-q[1] > 1  then
           resc[p[1]-p[2]+q[2]-q[1]] := resc[p[1]-p[2]+q[2]-q[1]]+1 ;
         fi;
      ind := p[2]-p[1]-q[1]+q[2]+2*n;
      resp1[ind] := resp1[ind]-1;
      resp2[ind] := resp2[ind]-1;
    fi;
  od;
  return [resc, resq, resp1, resp2];
end;

gamma2 := function(la, p, n)
  local res, resp1, resp2, ind1, ind2, ind3, resc, resq;
  res := gammaAR(la, p, n);
  resc :=res[1];
  resq :=res[2];
  resp1 := res[3];
  resp2 := res[4];
  # mu_k = 1 or not
  if p[2] = 1 then
    resq[1] := resq[1]+4*p[1]-4*p[2]+1;
    resq[2] := resq[2]-1;
    resq[3] := resq[3]-1;
    ind1 := 2*p[2]-2*p[1]-1+2*n;
    resp1[ind1] := resp1[ind1]+1;
    ind2 := 2*p[2]-2*p[1]+1+2*n;
    resp2[ind2] := resp2[ind2]+1;
  else
    resq[1] := resq[1]+6*p[1]-2*p[2]-1;
    resq[2] := resq[2]-1;
    resq[3] := resq[3]-1;
    ind1 := 2*p[2]-2*p[1]-1+2*n;
    resp1[ind1] := resp1[ind1]+1;
    ind2 := 2*p[2]-2*p[1]+1+2*n;
    resp2[ind2] := resp2[ind2]+1;
    ind3 :=p[2]-2*p[1]+2*n;
    resp1[ind3] :=resp1[ind3]+1;
    resp2[ind3] :=resp2[ind3]+1;
  fi;
  return res;
end;

gamma3 := function(la, p, n)
  local res, resc, resq, resp1, resp2,mu, ind, q,ind1, ind2;
  res := gammaAR(la, p, n);
  resc := res[1];
  resq := res[2];
  resp1 := res[3];
  resp2 := res[4];
  mu := ShallowCopy(la);
  mu[p[1]] := p[2];
  for q in AddableNodes(mu) do
    if q[1] > p[1] then
      ind := -(p[1]-p[2]+q[1]-q[2])+2*n;
      resp1[ind] := resp1[ind]+1;
      resp2[ind] := resp2[ind]+1;
    fi;
  od;
  for q in RemovableNodes(la) do
    if q[1] > p[1] then
      ind := -(p[1]-p[2]+q[1]-q[2])+2*n;
      resp1[ind] := resp1[ind]-1;
      resp2[ind] := resp2[ind]-1;
    fi;
  od;
    resq[1]:= resq[1]+6*p[1]-2*p[2]-1;
    resq[2]:= resq[2]-1;
    resq[3]:= resq[3]-1;
  ind1 := 2*p[2]-2*p[1]+1+2*n;
  resp2[ind1] := resp2[ind1]+1;
  ind2 := 2*p[2]-2*p[1]-1+2*n;
  resp1[ind2] := resp1[ind2]+1;
  if not [p[1],p[2]-1] in RemovableNodes(la) then
    resq[1]:= resq[1]+2*p[2]-2*p[1]-2;
    ind := -1+2*p[2]-2*p[1]+2*n;
    resp1[ind] := resp1[ind]-1;
    resp2[ind] := resp2[ind]-1;
  fi;
  return res;
end;


CACHEdetGfla := [[],[]];
adjustedCACHEdetGfla := function(res, nn)
  local m1,m2, v1,v2, a1,a2, i;
  m1 := Length(res[3])/4;
  if m1 = nn then
    return res;
  fi;
  v1 := 0*[1..4*nn];
  for i in [1..4*m1] do
    if res[3][i] <> 0 then
      a1 := i - 2 * m1 + 2 * nn;
      if a1 < 1 or a1 > 4*nn then
        Error("Cannot adjust length of cache result.");
      fi;
      v1[a1] := res[3][i];
    fi;
  od;
   m2 := Length(res[4])/4;
  if m2 = nn then
    return res;
  fi;
  v2 := 0*[1..4*nn];
  for i in [1..4*m2] do
    if res[4][i] <> 0 then
      a2 := i - 2 * m2 + 2 * nn;
      if a2 < 1 or a2 > 4*nn then
        Error("Cannot adjust length of cache result.");
      fi;
      v2[a2] := res[4][i];
    fi;
  od;
   return [res[1],res[2],v1, v2];
end;

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;


detGfla := function(arg)
  local f, la, n, nn, res, mu, fac1, fac2, exp, p, pos,i,j,k;
  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*[1..n],0*[1..3], 0*[1..4*nn], 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,nn);
    exp := rankDelta(f,mu);
    res[1] := res[1]+exp*fac2[1]+fac1[1];
    res[2] := res[2]+exp*fac2[2]+fac1[2];
    res[3] := res[3]+fac1[3];
    res[4] := res[4]+fac1[4];
  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] := res[1]+fac1[1]+exp*fac2[1];
      res[2] := res[2]+fac1[2]+exp*fac2[2];
      res[3] := res[3]+fac1[3]+exp*fac2[3];
      res[4] := res[4]+fac1[4]+exp*fac2[4];
    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;


############################################################################
stringFactoredDet := function(res)
  local str, nn1, nn2,a1,a2, i,j;
  str := "";
  if res[1] <> 0*res[1] then
    if res[1][1] mod 2 <> 0 then
      Add(str, '-');
    fi;
   fi;
   for i in [2..Length(res[1])] do
      if res[1][i] <> 0 then
        Add(str, '[');
        Append(str, String(i));
        Add(str, ']');
        if res[1][i] <> 1 then
          Add(str, '^');
          Append(str, String(res[1][i]));
        fi;
        Append(str, " * ");
      fi;
    od;


  if res[2] <> 0*res[2] then
     if res[2][1] <> 0 then
       Append(str, "q^");
       Append(str, String(res[2][1]));
      fi;
      if res[2][2] <> 0 then
       Append(str, "r^");
       Append(str, String(res[2][2]));
      fi;
     if res[2][3] <> 0 then
       Append(str, "(q^2 - 1)^");
       Append(str, String(res[2][3]));
     fi;
   fi;

  nn1 := Length(res[3])/4;
  for i in [1..4*nn1] do
    if res[3][i] <> 0 then
      a1 := i-2*nn1;
        if a1 = 0 then
            Append(str, "(r+1)");
         elif a1 = 1 then
            Append(str, "(rq+1)");
         elif a1 > 1 then
            Append(str, "(rq^");
            Append(str, String(a1));
            Append(str, "+1)");
         else
         Append(str, "(rq^(");
         Append(str, String(a1));
         Append(str, ")+1)");
        fi;

       if res[3][i] <> 1 then
         Add(str, '^');
        Append(str, String(res[3][i]));
        fi;
        Append(str," * ");
    fi;
  od;

 nn2 := Length(res[4])/4;
  for j in [1..4*nn2] do
    if res[4][j] <> 0 then
      a2 := j-2*nn2;
        if a2 = 0 then
            Append(str, "(r-1)");
        elif a2 > 0 then
            Append(str, "(rq^");
            Append(str, String(a2));
            Append(str, "-1)");
        else
         Append(str, "(rq^(");
         Append(str, String(a2));
         Append(str, ")-1)");
        fi;

       if res[4][j] <> 1 then
         Add(str, '^');
        Append(str, String(res[4][j]));
        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;



#############################################




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;