######################
# This is code implements my algorithm to compute the canonical basis for type C blocks
#
# Load it up using GAP4. There shouldn't be any errors.
# Once successfully loaded  you set global variables giving the sign sequence sigma. It defaults to 
#
# gap> Setup(4);
#
# Then call the function 
#
# gap> C([1,2,-2,-1]);
#
# for example to compute the canonical basis b_{(1,2,-2,-1)}.
#
# Obviously the code will only work for pretty small n's... then it will start running out of memory or taking too long.
#
######################

if not IsBound(n) then
  n := 4;
  Indeterminate(Rationals,"q");
  q := LaurentPolynomialByCoefficients(FamilyObj(1), [1], 1);
fi;

Setup := function(nn)
  n := nn;
end;

# The basic data type is a term, which is an n|m-tuple, and then a vec, which
# is a list of pairs [cf,term]

V := function(term) local string,i;
  return [[q^0,ShallowCopy(term)]];
end;

# This returns v1 + c * v2
# It corrupts v2 but not v1

AddV := function(v1,v2,c) local a,i,cf,newv;
  newv := [];
  for a in v1 do 
    cf := a[1];
    for i in [1..Length(v2)] do if IsBound(v2[i]) then
      if v2[i][2] = a[2] then cf := cf + c * v2[i][1];Unbind(v2[i]);fi;
    fi;od;
    if cf <> 0*q then Add(newv,[cf,a[2]]);fi;
  od;
  for a in v2 do
    Add(newv,[a[1]*c,a[2]]);
  od;
  return newv;
end;

# Print a vector in a reasonable format

PrintV := function(v) local term,i;
  for term in v do
    Print("V(");
    for i in [1..n] do
      Print(term[2][i]);
      if i < n then Print(",");fi;
    od;
    Print(")");
    if term[1] <> q^0 then Print(" * [",term[1],"]\n");else Print("\n");fi;
  od;
end;

###########################
# This is the Bruhat ordering
# Returns true if term1 >= term2

Dominates := function(term1,term2) local i,psi1,psi2,r,biggie;
  biggie := Maximum(Maximum(Maximum(term1), Maximum(term2)), 1-Minimum(term1), 1-Minimum(term2));
  for i in [0..biggie] do
    psi1 := 0; psi2 := 0;
    for r in [1..n] do
      if term1[r] > i then
        psi1 := psi1 + 1;
      elif term1[r] <= -i then
        psi1 := psi1 - 1;
      fi;
      if term2[r] > i then
        psi2 := psi2 + 1;
      elif term2[r] <= -i then
        psi2 := psi2 - 1;
      fi;
      if psi1 < psi2 then return(false);fi;
      if i = 0 and EuclideanRemainder(psi1-psi2,2) <> 0 then return false;fi;
    od;
    if psi1 <> psi2 then return(false);fi;
  od;
  return(true);
end;

#############
# This acts by K_i in the rth position of the tensor
# (r = 1,...,n).

# This acts by K_i in the rth position

ActKp := function(vec,i,r,power,rescale) local ans,j,cf,this;
  ans := [];
  for j in [1..Length(vec)] do if IsBound(vec[j]) then
    this := vec[j][2];
    cf := vec[j][1]*rescale;
    if i = 0 then
      if this[r] = 0 then cf := cf * q^(2*power);
      elif this[r] = 1 then cf := cf * q^(-2*power);
      fi;
    else
      if this[r] = i or this[r] = -i then  cf := cf * q^(power);
      elif this[r] = 1-i or this[r] = i+1 then cf := cf * q^(-power);
      fi;
    fi;
    ans := AddV(ans,V(this), cf);
  fi;od;
  return ans;
end;

ActKm := function(vec,i,r,power,rescale) local ans,j,cf,this;
  ans := [];
  for j in [1..Length(vec)] do if IsBound(vec[j]) then
    this := vec[j][2];
    cf := vec[j][1]*rescale;
    if i = 0 then
      if this[r] = 0 then cf := cf * q^(-2*power);
      elif this[r] = 1 then cf := cf * q^(2*power);
      fi;
    else
      if this[r] = i or this[r] = -i then  cf := cf * q^(-power);
      elif this[r] = 1-i or this[r] = i+1 then cf := cf * q^(power);
      fi;
    fi;
    ans := AddV(ans,V(this), cf);
  fi;od;
  return ans;
end;

# This acts by E_i in the rth position of the tensor

ActEsub := function(vec,i,r) local ans,j,this;
  ans := [];
  for j in [1..Length(vec)] do if IsBound(vec[j]) then
    this := ShallowCopy(vec[j][2]);
    if i = 0 and this[r] = 1 then
      this[r] := 0;
      ans := AddV(ans,V(this),vec[j][1]);
    elif this[r] = i+1 or this[r] = 1-i then
      this[r] := this[r]-1;
      ans := AddV(ans,V(this),vec[j][1]);
    fi;
  fi;od;
  return ans;
end;

# This acts by F_i in the rth position of the tensor

ActFsub := function(vec,i,r) local ans,j,this;
  ans := [];
  for j in [1..Length(vec)] do if IsBound(vec[j]) then
    this := ShallowCopy(vec[j][2]);
    if i = 0 and this[r] = 0 then
      this[r] := 1;
      ans := AddV(ans,V(this),vec[j][1]);
    elif this[r] = i or this[r] = -i then
      this[r] := this[r]+1;
      ans := AddV(ans,V(this),vec[j][1]);
    fi;
  fi;od;
  return ans;
end;

# Given a vector in T^{n} and a list of r positions, this acts
# K_i^{-r}\otimes...\otimes K_i^{-r}\otimes q^r K_i^{-r} E_i\otimes K_i^{1-r}\otimes ...\otimes K_i^{1-r}\otimes q^{r-1} K_i^{1-r} E_i
# \otimes...\otimes E_i \otimes 1 \otimes... \otimes 1 

ActMultipleE := function(vec,i,position) local mult,r,j;
  mult := 0;
  for r in Reversed([1..n]) do
    if r = position then
      vec := ActEsub(vec,i,r);
      vec := ActKm(vec,i,r,mult,q^mult);
      mult := mult + 1;
    else
      vec := ActKm(vec,i,r,mult,1);
    fi;
  od;
  return(vec);
end;

ActMultipleF := function(vec,i,position) local mult,r,j;
  mult := 0;
  for r in [1..n] do
    if r = position then
      vec := ActFsub(vec,i,r);
      vec := ActKp(vec,i,r,mult,q^mult);
      mult := mult + 1;
    else
      vec := ActKp(vec,i,r,mult,1);
    fi;
  od;
  return(vec);
end;

ActE := function(vec, i) local v,position,t;
  v := [];
  for position in [1..n] do
    v := AddV(v,ActMultipleE(vec,i,position),1);
  od;
  return(v);
end;

ActF := function(vec, i) local v,position,t;
  v := [];
  for position in [1..n] do
    v := AddV(v,ActMultipleF(vec, i,position),1);
  od;
  return(v);
end;


#############
# Here's the implementation of the main algorithm explained in our paper

C := function() end;	      # Forward reference

Correct := function(vec,leading) local nextleading,p,leadingok,nextcf,t,poly,i,this;

# Scan through all the terms.
# Check that the leading coefficient is 1, and find the 
# next lowest term in dominance that involves 1+(stuff) as its coeff.

  nextleading := false;
  leadingok := false;
  for t in [1..Length(vec)] do if IsBound(vec[t]) then
    this := vec[t][2];
    if not Dominates(this, leading) then Error("Not lower!");fi;
    if this = leading then
      if vec[t][1] <> q^0 then Error("Bad leading coefficient!");fi;
      leadingok := true;
    elif nextleading = false or (this <> nextleading and Dominates(nextleading,this)) then
      poly := CoefficientsOfLaurentPolynomial(vec[t][1]);	# See if it has coefficients in constant or negative powers of q
      if poly[2] <= 0 then
        nextleading := this;
        nextcf := poly;
      fi;
    fi;
  fi;od;

  if leadingok = false then Error("Wrong leading term!");fi;
  if nextleading = false then return(vec);fi;

  p := 0*q;
  for i in [nextcf[2]..0] do
    if i < 0 and IsBound(nextcf[1][i+1-nextcf[2]]) then
      p := p + (q^i + q^(-i))*(nextcf[1][i+1-nextcf[2]]);
    fi;
    if i = 0 and IsBound(nextcf[1][i+1-nextcf[2]]) then
      p := p + q^0*(nextcf[1][i+1-nextcf[2]]);
    fi;
  od;
#  Print("\nCorrecting by ", p, "*", nextleading, "\n");
  return Correct(AddV(vec,C(nextleading), -p),leading);
end;

###
# This is Lemma A.5 from the paper
# The routine CheckLemma displays the answer in a user-friendly way

Lemma := function(b) local a,r,s,this,done;
  a := [b[1]];
  for s in [2..n] do
    this := b[s];
    for r in [1..s-1] do
      if this >= a[r] then this := a[r]-1;fi;
      if this > -b[r] then this := -b[r];fi;
    od;
    done := false;
    while not done do
      done := true;
      for r in [1..s-1] do
        if a[r]+this = 1 then
	  this := this-1;
	  done := false;
	fi;
      od;
    od;
    Add(a,this);
  od;
  return a;
end;

Word := function(a,b) local w,i;
  w := [];
  if a >= 0 then
    for i in [a..b-1] do Add(w,i);od;
  elif b <= 0 then
    for i in Reversed([1-b..-a]) do Add(w,i);od;
  else
    for i in Reversed([1..-a]) do Add(w,i);od;
    for i in [0..b-1] do Add(w,i);od;
  fi;
  return w;
end;

CheckLemma := function(b) local a,r,s,i;
  a := Lemma(b);
  for s in Reversed([1..n]) do
    Print("(");
    for i in Reversed(Word(a[s],b[s])) do Print("f_",i);od;
    Print(")");
    for r in [1..s-1] do
      if a[s] in Word(a[r],b[r]) or -a[s] in Word(a[r],b[r]) then
        Error("Here!");
      fi;
    od;
  od;
  for i in a do Print("v_", i);od;
  Print("\n");
end;


###
# Here is the main routine

C := function(b) local a, smaller,i,ans,lasts,r;
  if n = 1 then return V(b);fi;
  a := Lemma(b);
  smaller:=ShallowCopy(b);
  Unbind(smaller[n]);
  n := n-1;
  ans:= C(smaller);
  n:=n+1;
  lasts := b[n];
  for i in [1..Length(ans)] do
    for r in [1..n-1] do
      if lasts > -AbsInt(ans[i][2][r]) then lasts := -AbsInt(ans[i][2][r]);fi;
    od;
  od;
#  Print("Tensoring with ", lasts, "\n");
  for i in [1..Length(ans)] do
    Add(ans[i][2], lasts);
  od;
  for i in Word(lasts,b[n]) do
    ans := ActF(ans,i);
  od;
  return Correct(ans,b);
end;