#C22.txt, Nov. 21, 2016, Math 640 (Rutgers, NB, Dr. Z.)
#Random generation of a standard Young tableau of a given shape
Help:=proc(): print(` Kids(L), F(L), LD(R), LDal(R) , RYT(L) `): end:

 Kids:=proc(L)  local i,S:
 option remember:
  if L=[] then
    RETURN([]):
  fi:

    S:=[]:

    for i from 1 to nops(L)-1 do
      if L[i]>L[i+1] then
        S:=[op(S), [[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)],i]]:
      fi:
    od:

      for i from nops(L) to nops(L) do
      if L[i]>1 then
        S:=[op(S), [[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)],i]]:
       else
         S:=[op(S),[[op(1..i-1,L)],i]]:
      fi:
    od:
   
   
 end:

 #F(L): the number of Young tableaux of shape L
  F:=proc(L) local K,k:
   option remember:

   if L=[] then
      RETURN(1):
   fi:

   K:=Kids(L):

     add(F(k[1]),k in K):
   
  end:
    
     #LDdz(R): inputs a list of positive integers R, and outputs a random integer from
     #1 to nops(R), such that the prob. op i is R[i]/convert(R,`+`): (depracated)
     LDdz:=proc(R) local ra,su,p,i,j:
      su:=convert(R,`+`):
       p:=rand(1..su)():
    
        for i from 1 to nops(R) while p>add(R[j],j=1..i) do od:
          i:
         end:


     #LDaz(R): Anthony Zaleski's inputs a list of positive integers R, and outputs a random integer from
     #1 to nops(R), such that the prob. op i is R[i]/convert(R,`+`):
      #Elegant  but inefficient
    LDaz:=proc(R) local i: [seq(i$R[i],i=1..nops(R))][rand(1..convert(R,`+`))()]:end:

   #LD(R): Andrew Lohr's version
   #inputs a list of positive integers R, and outputs a random integer from
     #1 to nops(R), such that the prob. op i is R[i]/convert(R,`+`):
     LD:=proc(R) local su,p,i:
      su:=convert(R,`+`):
       p:=rand(1..su)():
    
        for i from 1 to nops(R) while p>0 do p:=p-R[i] : od:
          i-1:
         end:
    #RYT(L): inputs a shape (alias partition) L, and outputs a (unfiormly at) random
    #Young tableau of shape L
    RYT:=proc(L) local n,K,R,p,k,L1,Y1,i1:
      if L=[] then
         RETURN([]):
      fi:
     n:=convert(L,`+`):
     
     K:=Kids(L):
     R:=[seq(F(K[i1][1]),i1=1..nops(K))]:
     
      p:=LD(R):
       k:=K[p][2]:
        L1:=K[p][1]:
        Y1:=RYT(L1):
         
          if k<=nops(Y1) then
        [op(1..k-1,Y1),[op(Y1[k]),n],op(k+1..nops(Y1),Y1)]:
          else
      
        [op(1..k-1,Y1),[n]]:
       fi:
     
     end:

        #PickCorner(L): One step of the Greene-Nijenhuis-Wilf algorithm
        #PickCorner:=proc(L) local c:
         #Under Construction
        #end: