with(combinat, partition):
# Try: [seq(AllNuYT3D(n,2),n=1..15)];
# Input: number of cell n and power k.
# Output: sum of (f_lambda)^k.
# Try: AllNuYT3D(5,1);
AllNuYT3D := proc(n,k) local p;
               add( NuYT3D(p)^k, p = Shape3D(n));
             end:

# Input: integer n
# Output: Set of shapes of Young diagram
# Try: Shape2D(3);
Shape2D := proc(n) local p;
             {seq(sort(p,`>`), p = partition(n))};
           end:

# Input: integer n
# Output: Set of shapes of Young diagram
# Try: Shape3D(3);
Shape3D := proc(n) option remember; local i,p;
             {seq(seq(op(Part3D([p],n-i)) ,p = Shape2D(i)),i=1..n)};
           end:

# Input: the shape of unfinished tableau in 3D
# with number of cell left.
# Output: the finished shape of Tableau 
# with more height. 
# Try: Part3D([[2,1]],2);
Part3D := proc(Tab,n) option remember; local t,m,i,p:
            if n=0 then return({Tab}); fi:
            t := add(m, m = Tab[nops(Tab)]);
            {seq(seq(`if`( IsGood(Tab[nops(Tab)],p),
              op(Part3D([op(Tab),p],n-i)),NULL) 
              ,p = Shape2D(i)),i=1..min(t,n))};         
          end:

# Input: List S and T
# Ouput: true if S[i]>=T[i] for all i
# Try: IsGood([5,3,2],[4,1,1]);
IsGood := proc(S,T) local i,n,m;
            n := nops(S);
            m := nops(T);
            if (m > n) or ({seq( `if`(S[i]<T[i],1,NULL),i=1..m)} <> {}) then
                       return(false); fi:
            true;
          end:

#####################################################
# Input: Shape of lambda
# Output: The number of standard Young Tableaux
# Try: NuYT3D([[2,1],[1]]);
NuYT3D := proc(S) option remember; local T;
            if S = [[1]] then return(1); fi:
            add(NuYT3D(T), T = SubShape3D(S));
          end:

# Input: Shape of lambda
# Output: the set of lambda after take out corner
# Try: SubShape3D([[5,2],[1]]);
SubShape3D := proc(S) option remember; local i,j,n,ret  ;
                ret := {};
                for i from 1 to nops(S) do
                      n := nops(S[i]):
                      for j from 1 to n do
                              if (j=n or S[i][j] > S[i][j+1]) and
                              (i=nops(S) or nops(S[i+1]) < j 
                               or S[i][j] > S[i+1][j]) then
                                      ret := ret union {Dog3D(S,i,j)};
                              fi:
                      od:
                od:
                ret;
               end:

# Output: S without the cell [i,j].
# Try: Dog3D([[5,2],[1]],2,1);
Dog3D := proc(S,k,l) option remember;
           if k = nops(S) and S[k]=[1] then [op(1..nops(S)-1,S)];
              else [op(1..k-1,S),Dog(S[k],l),op(k+1..nops(S),S)];
           fi:
         end:
# Try: Dog([5,2,1],3);
Dog := proc(S,k) option remember;
         if k = nops(S) and S[k]=1 then [op(1..nops(S)-1, S)]
       else [op(1..k-1, S), S[k]-1, op(k+1..nops(S), S)]
         fi
       end:
a:= n-> AllNuYT3D(n,2):
seq (a(n), n=1..15);