#OK to post
#Yuxuan Yang, March 6th, Assignment 12

with(combinat):

#1
#The set of partitions of n with largest part k where all parts are congruent to a member of S mod m.
PnkC:=proc(n,k,S,m) local k1,L,L1:
option remember:
if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN([]):
fi:
if not (k mod m) in S then
 RETURN([]):
fi:
if k>n then
 RETURN([]):
fi:
if k=n then
 RETURN([[n] ]):
fi:
L:=[]:

for k1 from min(n-k,k) by -1 to 1 do
 L1:=PnkC(n-k,k1,S,m):
 L:=[op(L), seq([k, op(L1[j])],j=1..nops(L1))]:
od:
L:
end:


#PnC(n,S,m): The set of partitions of n where all parts are congruent to a member of S mod m
PnC:=proc(n,S,m) local k:option remember:[seq(op(PnkC(n,n-k+1,S,m)),k=1..n)]:end:



#pnkC(n,k,S,m):The NUMBER of partitions of n with largest part k where all parts are congruent to a member of S mod m.
pnkC:=proc(n,k,S,m) local k1:
option remember:
if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN(0):
fi:
if not (k mod m) in S then
 RETURN(0):
fi:
if n=k then
 RETURN(1):
else
 RETURN(add(pnkC(n-k,k1,S,m),k1=1..min(k,n-k))):
fi:
end:

#pnC(n,S,m):The NUMBER of partitions of n where all parts are congruent to a member of S mod m.
pnC:=proc(n,S,m) local k: option remember: 
if n=0 then
 RETURN(1):
fi:
add(pnkC(n,k,S,m),k=1..n):end:




#RandPnkC(n,k,S,m): A random partition of n with largest part k where all parts are congruent to a member of S mod m.
RandPnkC:=proc(n,k,S,m) local die1, k1,k2:
if not (k mod m) in S then
 RETURN(FAIL):
fi:
if n=1 then 
RETURN([1]):
fi:
if n=k then 
RETURN([k]):
fi:
die1:=[seq(pnkC(n-k,k1,S,m),k1=1..min(n-k,k))]:
k2:=LD(die1):
[k, op(RandPnkC(n-k,k2,S,m))]:
end:

#RandPnC(n,S,m): A random partition of n where all parts are congruent to a member of S mod m
RandPnC:=proc(n,S,m) local k,i,die1,k1:
die1:=[seq(pnkC(n,k,S,m),k=1..n)]:
k1:=LD(die1):
RandPnkC(n,k1,S,m):
end:

#2
#PnkD(n,k,DI): The set of partitions of n with largest part k where the difference of two consecutive parts is NOT in DI
PnkD:=proc(n,k,DI) local k1,L,L1:
option remember:
if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN([]):
fi:
if k>n then
 RETURN([]):
fi:
if k=n then
 RETURN([[n]]):
fi:
L:=[]:

for k1 from min(n-k,k) by -1 to 1 do
 if not k-k1 in DI then 
  L1:=PnkD(n-k,k1,DI):
  L:=[op(L), seq([k, op(L1[j])],j=1..nops(L1))]:
 fi:
od:
L:
end:

#PnD(n,DI): The set of partitions of n where the difference of two consecutive parts is NOT in DI
PnD:=proc(n,DI) local k:option remember:[seq(op(PnkD(n,n-k+1,DI)),k=1..n)]:end:



#pnkD(n,k,DI):The NUMBER of partitions of n where the difference of two consecutive parts is NOT in DI
pnkD:=proc(n,k,DI) local k1,nu:
option remember:
if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN(0):
fi:
if n=k then
 RETURN(1):
else
 nu:=0:
 for k1 from 1 to min(k,n-k) do
  if not k-k1 in DI then
   nu:=nu+pnkD(n-k,k1,DI):
  fi:
 od:
 RETURN(nu):
fi:
end:

#pnD(n,DI):The NUMBER of partitions of n where the difference of two consecutive parts is NOT in DI
pnD:=proc(n,DI) local k: option remember: 
if n=0 then
 RETURN(1):
fi:
add(pnkD(n,k,DI),k=1..n):end:




#RandPnkD(n,k,DI): A random partition of n with largest part k where the difference of two consecutive parts is NOT in DI
RandPnkD:=proc(n,k,DI) local die1, k1,k2,klist:
if pnkD(n,k,DI)=0 then 
 RETURN(FAIL):
fi:
if n=1 then 
RETURN([1]):
fi:
if n=k then 
RETURN([k]):
fi:
die1:=[seq(pnkD(n-k,k1,DI),k1=1..min(n-k,k))]:
k2:=LD(die1):
if not k-k2 in DI then
RETURN([k, op(RandPnkD(n-k,k2,DI))]):
else

RETURN(RandPnkD(n,k,DI)):
fi:
end:

#RandPnD(n,DI): A random partition of n where the difference of two consecutive parts is NOT in DI
RandPnD:=proc(n,DI) local k,i,die1,k1:
die1:=[seq(pnkD(n,k,DI),k=1..n)]:
k1:=LD(die1):
RandPnkD(n,k1,DI):
end:



#[seq(pnC(n,{1},2),n=1..30)]
#[1, 1, 2, 2, 3, 4, 5, 6, 8, 10, 12, 15, 18, 22, 27, 32, 38, 46, 54, 64, 76, 89, 104, 122, 142, 165, 192, 222, 256, 296]
#A000009

#[seq(pnC(n,{1,4},5),n=1..30)]
#[1, 1, 1, 2, 2, 3, 3, 4, 5, 6, 7, 9, 10, 12, 14, 17, 19, 23, 26, 31, 35, 41, 46, 54, 61, 70, 79, 91, 102, 117]
#A003114

#[seq(pnD(n,{0}),n=1..30)]
#[1, 1, 2, 2, 3, 4, 5, 6, 8, 10, 12, 15, 18, 22, 27, 32, 38, 46, 54, 64, 76, 89, 104, 122, 142, 165, 192, 222, 256, 296]
#A000009       amazing!

#[seq(pnD(n,{0,1}),n=1..30)]
#[1, 1, 1, 2, 2, 3, 3, 4, 5, 6, 7, 9, 10, 12, 14, 17, 19, 23, 26, 31, 35, 41, 46, 54, 61, 70, 79, 91, 102, 117]
#A003114








#LD(L): Inputs a list of positive integers L (of n:=nops(L) members)
#outputs an integer i from 1 to n with the prob. of i being
#proportional to L[i]
#For example LD([1,2,3]) should output 1 with prob. 1/6
#output 2 with prob. 1/3
#output 3 with prob. 3/6=1/2
LD:=proc(L) local n,i,su,r:
n:=nops(L):
r:=rand(1..convert(L,`+`))():
su:=0:

for i from 1 to n do
 su:=su+L[i]:
  if r<=su then
    RETURN(i):
  fi:
od:
end: