######################
# This is code implements my algorithm to compute the canonical basis for type A 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([1,1,-1,-1]);
#
# Then call the function 
#
# gap> A([1,2,2,1]);
#
# for example to compute the canonical basis b_{(1,2|2,1)}. Note this example is not multiplicity free!!!
#
# Obviously the code will only work for pretty small sig's... then it will start running out of memory or taking too long.
#
######################

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

Setup := function(sigma)
  sig := sigma; n := Length(sig);
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,a;
  for a in [Minimum(Minimum(term1), Minimum(term2))..Maximum(Maximum(term1), Maximum(term2))] do
    psi1 := 0; psi2 := 0;
    for i in [1..n] do
      if term1[i] > a then
        psi1 := psi1 + sig[i];
      fi;
      if term2[i] > a then
        psi2 := psi2 + sig[i];
      fi;
      if psi1 < psi2 then return(false);fi;
    od;
    if psi1 <> psi2 then return(false);fi;
  od;
  return(true);
end;


###############################
# Here are routines to act with the quantum group

# 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 sig[r]=1 then
      if this[r]=i then cf := cf * (q^power);fi;
      if this[r]=i+1 then cf := cf / (q^power);fi;
    else
      if this[r]=i then cf := cf / (q^power);fi;
      if 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 K_i^{-1} on the rth position

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 sig[r]=1 then
      if this[r]=i then cf := cf / (q^power);fi;
      if this[r]=i+1 then cf := cf * (q^power);fi;
    else
      if this[r]=i then cf := cf * (q^power);fi;
      if 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 on the rth position

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 sig[r]=1 and this[r] = i+1 then
      this[r] := i;
      ans := AddV(ans,V(this),vec[j][1]);
    elif sig[r]=-1 and this[r] = i then
      this[r] := i+1;
      ans := AddV(ans,V(this),vec[j][1]);
    fi;
  fi;od;
  return ans;
end;

# This acts by F_i on the rth position

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 sig[r]=1 and this[r] = i then
      this[r] := i+1;
      ans := AddV(ans,V(this),vec[j][1]);
    elif sig[r]=-1 and this[r] = i+1 then
      this[r] := i;
      ans := AddV(ans,V(this),vec[j][1]);
    fi;
  fi;od;
  return ans;
end;

# Given a vector in T^{n|m} 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;
  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;
  v := [];
  for position in [1..n] do
    v := AddV(v,ActMultipleF(vec, i,position),1);
  od;
  return(v);
end;

##############
# This applies Lemma 4.3 from paper with Nick

Lemma := function(b) local a,r,s,this;
  a := [b[1]];
  for s in [2..n] do
    this := b[s];
    if sig[s] = 1 then
      for r in [1..s-1] do
        if sig[r] = -1 and b[r]-1 < this then this := b[r]-1;fi;
        if sig[r] = 1 and a[r]-1 < this then this := a[r]-1;fi;
    od;
    else
      for r in [1..s-1] do
        if sig[r] = 1 and b[r]+1 > this then this := b[r]+1;fi;
        if sig[r] = -1 and a[r]+1 > this then this := a[r]+1;fi;
      od;
    fi;
    Add(a,this);
  od;
  return a;
end;


####

A := 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,A(nextleading), -p),leading);
end;

A := function(term) local i,newt,lasts,lastt,ans,r;
  n:=Length(sig);
  if Length(sig) = 1 then return(V(term));fi;
  newt := ShallowCopy(term);
  lasts := sig[n];
  lastt := term[n];
  Unbind(newt[n]);
  Unbind(sig[n]); n := n-1;
  ans := A(newt);
  for i in [1..Length(ans)] do
    if lasts=1 then
      for r in [1..n] do
        if sig[r] = 1 and ans[i][2][r] < lastt then lastt := ans[i][2][r];fi;
        if sig[r] = -1 and ans[i][2][r] <= lastt then lastt := ans[i][2][r]-1;fi;
      od;
    else
      for r in [1..n] do
        if sig[r] = -1 and ans[i][2][r] > lastt then lastt := ans[i][2][r];fi;
        if sig[r] = 1 and ans[i][2][r] >= lastt then lastt := ans[i][2][r]+1;fi;
      od;
    fi;
  od;
  Print("Tensoring with ", lastt, "\n");
  for i in [1..Length(ans)] do
    Add(ans[i][2], lastt);
  od;
  Add(sig,lasts); n := n+1;
  if lasts = 1 then
    for i in [lastt..term[n]-1] do
      ans := ActF(ans,i);
    od;
  else
    for i in Reversed([term[n]..lastt-1]) do
      ans := ActF(ans,i);
    od;
  fi;
  return Correct(ans,term);
end;