with(combinat, partition):
# Input: number of cell n and power k.
# Output: sum of (f_lambda)^k.
AllNuYT2D := proc(n, k) local p;
               add(NuYT2D(p)^k, p = Shape2D(n));
             end:

# Input: Shape of lambda
# Output: The number of standard Young Tableaux
# Try: NuYT2D([2, 1]);
NuYT2D := proc(S) option remember; local T;
            if S <> sort(S,`>`) or S[nops(S)] < 0 then ERROR(S);
          elif S[nops(S)] = 0 then return(NuYT2D([op(1..nops(S)-1,S)]));
            fi:
            if S = [1] then return(1) fi:
            add (NuYT2D(T), T = SubShape(S))
          end:

# Input: Shape of lambda
# Output: the set of lambda after take out corner
# Try: SubShape([5,2,1]);
SubShape := proc(S) option remember; local i,n;
              n:= nops(S);
              {seq(`if`(i=n or S[i] > S[i+1], Dog(S, i), NULL), i = 1..n)}
            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:

# 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:

a:= n-> AllNuYT2D(n,4):
seq (a(n), n=1..20);