with(combinat):
read RSK;

#-------------------------- Intro ----

print(`Welcome to NPS.  For help type ``Help();```);
print(``);

#-------------------------- Help -----------------

Help:=proc()
   if nargs = 0 then
     print(`Welcome to NPS, a Maple package implementing the`);
     print(`Novelli-Pak-Stoyanovskii bijective proof of the`);
     print(`Hook-Length Formula.`);
     print(``);
     print(`Version of May 4, 2003 `);
     print(``);
     print(`Procedures available: NPS, SPN, YT, RandNPS`);
   fi;
   if nargs = 1 then
     if args[1] = NPS then
       print(`NPS(T): Inputs a (not necessarily standard) Young 
tableau,`);
       print(`and outputs the pair [P, J] given by the 
Novelli-Pak-Stroyanovskii`);
       print(`bijection.  P is a stardard Young tableau and J is a "hook 
tableau."`);
       print(`NPS may also be applied to sets of tableau.`);
       print(`The inverse of NPS is SPN.`);
     fi;
     if args[1] = SPN then
       print(`SPN([P, J]): Inputs a standard Young tableau P and a "hook 
tableau" J`);
       print(`and outputs the tableau given by the inverse of the`);
       print(`Novell-Pak-Stroyanovskii bijection.`);
       print(`SPN may also be applied to sets of pairs of tableau.`);
       print(`The inverse of SPN in NPS.`);
     fi;
     if args[1] = YT then
       print(`YT(lambda): Inputs a partition and outputs all (not 
necessarily`);
       print(`standard) Young tableau of shape lambda.`);
     fi;
     if args[1] = RandNPS then
       print(`RandNPS(lambda): Inputs a partition, creates a random 
tableau T`);
       print(`of that shape, and checks whether T = SPN(NPS(T)).`);
       print(`Prints an error message if that is not the case, 
otherwise`);
       print(`outputs nothing.`);
     fi;
   fi;
   print(``);
end:

#-------------------------- Procedures for testing --------------------

Perm2YT := proc(pi, lambda)
   # Given the permutation pi and partition (of nops(pi)) lambda, returns 
#the Young
   # Tableau determined by pi and lambda.
   if pi = [] then
     return [];
   else
     return [[op(1..lambda[1], pi)],op(Perm2YT([op(lambda[1]+1..nops(pi), 
pi)], [op(2..nops(lambda), lambda)]))]
   fi;
end:

YT := proc(lambda)
   # Returns all the Young Tableau of shape lambda
   local pi, Perms, i, n, YT;
   n := convert([seq(lambda[i], i=1..nops(lambda))], `+`);
   Perms := permute(n);
   YT := {};
   for pi in Perms do
     YT := YT union {Perm2YT(pi, lambda)};
   od;
   return YT;
end:

RandPerm := proc(n)
   local pi, N, v;
   pi := [];
   N := {$1..n};
   while nops(N) > 0 do
     v := N[rand(1..nops(N))()];
     pi := [op(pi), v];
     N := N minus {v};
   od;
   return pi;
end:

RandNPS := proc(lambda)
   local n, pi, T1, T2;
   n := convert(lambda, `+`);
   pi := RandPerm(n);
   T1 := Perm2YT(pi, lambda);
   T2 := SPN(NPS(T1));
   if T1 <> T1 then
     ERROR(`Bijection failed for tableau `, T1);
   fi;
end:

#-------------------------- Helper procedures --------------------

LowerOrderTableau := proc(T, c)
   # Returns the tableau T^{\le c}, as described on p. 126.
   # T should be a list of lists, as usual, but the returned
   # tableau will be a matrix.
   local M, i, j, x, y;
   M := linalg[matrix](nops(T)+1, nops(T[1])+1, 0);
   # The "+1" above is just for convenience later...
   x := c[1];
   y := c[2];
   for i from 1 to nops(T) do
     for j from 1 to nops(T[i]) do
       if (j > y) or (j = y and i >= x) then
         M[i,j] := T[i][j];
       fi;
     od;
   od;
   return op(M);
end:

isStandard := proc(M)
   # Inputs a tableau in matrix form (such as LowerOrderTableau outputs), 
#and returns TRUE if it is
   # standard, FALSE otherwise.
   local i, j;
   for i from 1 to linalg[rowdim](M) do
     for j from 1 to linalg[coldim](M) do
       if M[i,j] <> 0 then
         if i > 1 then
           if M[i-1, j] > M[i, j] then
             return false;
           fi;
         fi;
         if j > 1 then
           if M[i, j-1] > M[i, j] then
             return false;
           fi;
         fi;
       fi;
     od;
   od;
   return true;
end:

ModifiedForwardSlide := proc(T, c)
   # T is a tableau, c is a cell.
   # Returns the output tableau, along with the end point of the path.
   local T1, M, i, j, c1, cprime, cvalue;
   c1 := c;
   M := LowerOrderTableau(T, c1);
   T1 := T;
   while not isStandard(M) do
     i := c1[1];
     j := c1[2];
     if M[i+1, j] = 0 then
       cprime := [i, j+1];
     elif M[i, j+1] = 0 then
       cprime := [i+1, j];
     elif M[i+1, j] > M[i, j+1] then
       cprime := [i, j+1];
     elif M[i, j+1] > M[i+1, j] then
       cprime := [i+1, j];
     fi;
     # Now we exchange cells c1 and cprime.
     cvalue := M[op(c1)];
     M[op(c1)] := M[op(cprime)];
     M[op(cprime)] := cvalue;
     # Now we must also swap c1 and cprime in T1.
     cvalue := T1[c1[1]][c1[2]];
     T1[c1[1]][c1[2]] := T1[cprime[1]][cprime[2]];
     T1[cprime[1]][cprime[2]] := cvalue;
     c1 := cprime;
   od;
   return T1, c1;
end:

CellOrder := proc(T)
   # Given a tableau T, returns a list of the cells in the order defined 
#on p. 126.
   local T1, i;
   if T = [] then
     return [];
   fi;
   i := 1;
   while i < nops(T) and nops(T[i+1]) = nops(T[1]) do
     i := i + 1;
   od;
   if nops(T[i]) > 1 then
     T1 := T;
     T1[i] := [op(1..nops(T[i])-1, T[i])];
   else
     T1 := [op(1..nops(T)-1, T)];
   fi;
   return [[i,nops(T[i])], op(CellOrder(T1))];
end:

CanidateCells := proc(P, J, c)
   local C, i0, j0, ip, jp;
   i0 := c[1];
   j0 := c[2];
   C := {};
   for ip from i0 to nops(P) do
     if j0 <= nops(J[ip]) and J[ip][j0] >= 0 then
       C := C union {[ip, j0 + J[ip][j0]]};
     fi;
   od;
   return C;
end:

ModifiedBackwardSlide := proc(Tp, Jp, ck, c)
   local i0, j0, c1, i, j, cp, cvalue, Tp1, Path;
   Tp1 := Tp;
   i0 := ck[1];
   j0 := ck[2];
   c1 := c;
   Path := [c1];
   while c1 <> ck do
     i := c1[1];
     j := c1[2];
     if i - 1 <= 0 then
       cp := [i, j-1];
     elif j - 1 < j0 then
       cp := [i-1, j];
     elif Tp1[i-1][j] > Tp1[i][j-1] then
       cp := [i-1, j];
     elif Tp1[i][j-1] > Tp1[i-1][j] then
       cp := [i, j-1];
     fi;
     cvalue := Tp1[c1[1]][c1[2]];
     # Now we exchange Tp1[c1] and Tp1[cp] and swap c1 and cp.
     Tp1[c1[1]][c1[2]] := Tp1[cp[1]][cp[2]];
     Tp1[cp[1]][cp[2]] := cvalue;
     c1 := cp;
     Path := [op(Path), c1];
   od;
   return Tp1, Code(Path);
end:

Code := proc(Path)
   local i, C;
   C := [];
   for i from 1 to nops(Path) - 1 do
     if Path[i+1][1] = Path[i][1] - 1 then
       C := [op(C), `N`];
     else
       C := [op(C), `W`];
     fi;
   od;
   return [seq(C[nops(C) - i + 1], i=1..nops(C))];
end:

CodeEntryL  := proc(E1, E2)
   if E1 = `N` then
     if E2 = `N` then
       return false;
     else
       return true;
     fi;
   fi;
   if E1 = 0 then
     if E2 = `W` then
       return true;
     else
       return false;
     fi;
   fi;
   if E1 = `W` then
     return false;
   fi;
end:

CodeLE := proc(C1, C2)
   # Returns TRUE if the code C1 <= C2, false otherwise.
   local i, n, C11, C21;
   n := max(nops(C1), nops(C2));
   C11 := [op(C1), 0$(n - nops(C1))];
   C21 := [op(C2), 0$(n - nops(C2))];
   for i from 1 to n do
     if CodeEntryL(C11[i], C21[i]) then
       return true;
     fi;
     if CodeEntryL(C21[i], C11[i]) then
       return false;
     fi;
   od;
   return true;
end:

isCodeMax := proc(Codes, C)
   # Returns true if C is the maximum code in Codes, false otherwise.
   local D;
   for D in Codes do
     if not CodeLE(D, C) then
       return false;
     fi;
   od;
   return true;
end:

#-------------------------- The maps ---------------------------

NPS := proc(T)
   # Performs the NPS bijection.
   # Returns [P, J], where P is a standard tableau and J is a "hook 
#tableau."
   local n, P, J, J1, CO, c, k, h, Tem, ep, i, j, ip, jp;

   # NPS can also be applied to sets:
   if type(T, set) then
     return {seq(NPS(T[i]), i=1..nops(T))};
   fi;

   n := convert([seq(nops(T[i]), i=1..nops(T))], `+`);
   P := T;
   J := T;
   # J has the right dimensions, but we want it to start out at 0.
   for i from 1 to nops(T) do
     for j from 1 to nops(T[i]) do
       J[i][j] := 0;
     od;
   od;
   CO := CellOrder(T);
   for k from 1 to n do
     # We now build what Bruce calls (T_k, J_k).
     i := CO[k][1];
     j := CO[k][2];
     Tem := ModifiedForwardSlide(P, CO[k]);
     P := Tem[1];
     ep := Tem[2];
     ip := ep[1];
     jp := ep[2];
     J1 := J;
     for h from i to ip - 1 do
       J1[h][j] := J[h+1][j] - 1;
     od;
     J1[ip][j] := jp - j;
     J := J1;
   od;
   return [P, J];
end:

SPN := proc(X)
   local i, j, ip, jp, cp, c, h, k, n, CO, Codes, CodesWithCells, 
Canidates, Tem, P1, J1, J2;

   # SPN can also be applied to sets:
   if type(X, set) then
     return {seq(SPN(X[i]), i=1..nops(X))};
   fi;

   P1 := X[1];
   J1 := X[2];
   CO := CellOrder(P1);
   n := convert([seq(nops(P1[i]), i=1..nops(P1))], `+`);
   for k from n to 1 by -1 do
     Canidates := CanidateCells(P1, J1, CO[k]);
     Codes := {};
     CodesWithCells := {};
     for c in Canidates do
       Tem := ModifiedBackwardSlide(P1, J1, CO[k], c);
       Codes := Codes union {Tem[2]};
       CodesWithCells := CodesWithCells union {[Tem[2], c]};
     od;
     # Now we find the element with maximum code.
     for Tem in CodesWithCells do
       if isCodeMax(Codes, Tem[1]) then
         cp := Tem[2];
       fi;
     od;
     # Now cp is the cell with maximum code.
     Tem := ModifiedBackwardSlide(P1, J1, CO[k], cp);
     P1 := Tem[1];
     J2 := J1;
     i := CO[k][1];
     j := CO[k][2];
     ip := cp[1];
     jp := cp[2];
     for h from i + 1 to ip do
       J2[h][j] := J1[h-1][j] + 1;
     od;
     J2[i][j] := 0;
     # Now we swap J1 and J2.
     J1 := J2;
   od;
   return P1;
end: