###################### # This is code implements my algorithm to compute the canonical basis for gl(n|m) # # Load it up using GAP4. There shouldn't be any errors. # Once successfully loaded you set global variables n and m. It defaults to # # gap> n:=2;m:=2; # # Then call the function # # gap> B([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 n and m... then it will start running out of memory or taking too long. # ###################### if not IsBound(n) then startside := 1; n := 2; m := 2; Indeterminate(Rationals,"q"); q := LaurentPolynomialByCoefficients(FamilyObj(1), [1], 1); fi; # 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; if n > 0 then Print("|");fi; for i in [n+1..n+m] do Print(term[2][i]); if i < n+m 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 + 1; fi; if term2[i] >= a then psi2 := psi2 + 1; fi; if psi1 < psi2 then return(false);fi; od; for i in [1..m] do if term1[n+i] >= a then psi1 := psi1 - 1; fi; if term2[n+i] >= a then psi2 := psi2 - 1; fi; if psi1 < psi2 then return(false);fi; od; if psi1 <> psi2 then return(false);fi; od; return(true); end; ########################### # Here is the correction routine # You pass it the leading guy which must indeed be leading with coefficient 1 Bsub := function() end; # Forward reference Correct := function(vec,leading,truncate) 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("Correcting by ", p, "*", nextleading, "\n"); return Correct(AddV(vec,Bsub(nextleading,truncate), -p),leading,truncate); end; ########################### # You can act on vecs... ActC := function(vec,i) local ans,this,j,new; ans := []; for j in [1..Length(vec)] do if IsBound(vec[j]) then this := vec[j][2]; if i >= 1 and i < n then if this[i] = this[i+1] then ans := AddV(ans, V(this), vec[j][1]*(q+q^(-1))); elif this[i] < this[i+1] then new := ShallowCopy(this); new[i] := this[i+1]; new[i+1] := this[i]; ans := AddV(ans, V(new), vec[j][1]); ans := AddV(ans, V(this), vec[j][1]*q^(-1)); else new := ShallowCopy(this); new[i] := this[i+1]; new[i+1] := this[i]; ans := AddV(ans, V(new), vec[j][1]); ans := AddV(ans, V(this), vec[j][1]*q); fi; elif i > n and i <= n+m-1 then if this[i] = this[i+1] then ans := AddV(ans, V(this), vec[j][1]*(q+q^(-1))); elif this[i] < this[i+1] then new := ShallowCopy(this); new[i] := this[i+1]; new[i+1] := this[i]; ans := AddV(ans, V(new), vec[j][1]); ans := AddV(ans, V(this), vec[j][1]*q); else new := ShallowCopy(this); new[i] := this[i+1]; new[i+1] := this[i]; ans := AddV(ans, V(new), vec[j][1]); ans := AddV(ans, V(this), vec[j][1]*q^(-1)); fi; else Error("Bad index!");fi; fi;od; return ans; end; # 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 r <= n 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 r <= n 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 ActE := 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 r <= n and this[r] = i+1 then this[r] := i; ans := AddV(ans,V(this),vec[j][1]); elif r >= n+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 ActF := 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 r <= n and this[r] = i then this[r] := i+1; ans := AddV(ans,V(this),vec[j][1]); elif r >= n+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,positions) local mult,r,j; mult := 0; for r in Reversed([1..n+m]) do if r in positions then vec := ActE(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,positions) local mult,r,j; mult := 0; for r in [1..n+m] do if r in positions then vec := ActF(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; ActDividedE := function(vec, i, s) local v,positions; v := []; for positions in Combinations([1..n+m],s) do v := AddV(v,ActMultipleE(vec,i,positions),1); od; return(v); end; ActDividedF := function(vec, i, s) local v,positions; v := []; for positions in Combinations([1..n+m],s) do v := AddV(v,ActMultipleF(vec, i,positions),1); od; return(v); end; ############################ # This is a step in the bumping algorithm # Given a word thisside and i, it finds i...j all consecutive and adds one to them all, # tacking on [i,side*mult of i],...,[j,side*mult of j] to monomial and returning j+1 Bump := function(thisside,side,i,monomial) local oldside,j,mult,this; j := i; oldside := ShallowCopy(thisside); while j in thisside do j := j+1;od; while i < j do mult := 0; for this in [1..Length(oldside)] do if oldside[this] = i then thisside[this]:=thisside[this]+1;mult:=mult+1;fi; od; Add(monomial,[i,side*mult]); i := i+1; od; return j; end; SetTruncate := function(vec) local v,this,truncate; truncate := fail; for v in vec do this := v[2]; if truncate = fail or Maximum(this) > truncate then truncate := Maximum(this);fi; od; return truncate; end; Truncate := function(vec,truncate) local v,this,ans; if truncate = fail then return vec;fi; ans := []; for v in vec do this := v[2]; if Maximum(this) <= truncate then ans := AddV(ans, V(this), v[1]); else Print("Truncating!\n");fi; od; return ans; end; ############################ # This is the main routine # To optimize things I work in a finite interval # after the first bump. # This guarantees that the algorithm terminates theoretically. # Actually I suspect it always terminates anyway!!! # In other words the code never prints "Truncating!" which it does only if the truncation actually occurs Bsub := function(term,truncate) local i,new,intersect,top,side,thisside,left,right,vec,monomial,ans; # Find biggest match left := []; right := []; for i in [1..n] do Add(left,term[i]); od; for i in [n+1..n+m] do Add(right,term[i]); od; intersect := Intersection(left,right); if intersect <> [] then top := Maximum(intersect); monomial := []; # This is stored as a list of [i,m] where m > 0 is multiplicitiy of e and m < 0 is minus multiplicicity of f side := startside; # 1 for left hand side, -1 for right hand side if side = 1 then thisside := left; else thisside := right;fi; while top in thisside do top := Bump(thisside,side,top,monomial); side := -side; if side = 1 then thisside := left; else thisside := right;fi; od; Append(left,right); vec := Truncate(Bsub(left,fail),truncate); # Here to recurse with lower atypicality for i in monomial do if i[2] > 0 then vec := ActDividedE(vec,i[1],i[2]); else vec := ActDividedF(vec,i[1],-i[2]); fi; od; # And here when we correct we can be truncated ans := Correct(vec,term,SetTruncate(vec)); return ans; fi; # Here if typical for i in [1..n-1] do if term[i] < term[i+1] then new := ShallowCopy(term); new[i] := term[i+1]; new[i+1] := term[i]; ans := Correct(ActC(Bsub(new,fail),i), term,fail); return ans; fi; od; for i in [n+1..n+m-1] do if term[i] > term[i+1] then new := ShallowCopy(term); new[i] := term[i+1]; new[i+1] := term[i]; ans := Correct(ActC(Bsub(new,fail),i), term,fail); return ans; fi; od; # Here if typical and weakly dominant return V(term); end; # This just prints the output in a reasonable format... B := function(term) PrintV(Bsub(term,fail)); end;