#OK to post homework
#George Spahn, 3/7/2021, Assignment 12

#PnkC(n,k,S,m): The set of integer partitions of n with largest part k,
#where all parts are congruent to an element 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 k>n then
 RETURN({}):
fi:

if not(irem(k,m) in S) 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:=L union {seq([k, op(L1[j])],j=1..nops(L1))}:
od:
L:
end:

#Pn(n,S,m):  The set of integer 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 integer 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(irem(k,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 integer 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): random integer 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 Die,d,P:

if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN(FAIL):
fi:
if k>n then
 RETURN(FAIL):
fi:
if not(irem(k,m) in S) then
 RETURN(FAIL):
fi:
if k=n then
 RETURN([n]):
fi:

Die:=[seq(pnkC(n-k,k1,S,m),k1=1..min(n-k,k))]:
d:=LD(Die):
P:=RandPnkC(n-k,d,S,m):
[k, op(P)]:
end:


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


#####################################################################
#####################################################################



#PnkD(n,k,DI): The set of integer 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:=L union {seq([k, op(L1[j])],j=1..nops(L1))}:
 fi:
od:
L:
end:

#PnD(n,DI):  The set of integer 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 integer 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,sum:
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
 sum:=0:
 for k1 from min(n-k,k) by -1 to 1 do
  if not(k-k1 in DI) then
   sum:=sum+pnkD(n-k,k1,DI):
  fi:
 od:
fi:
sum:
end:

#pnD(n,DI) The NUMBER of integer 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): random integer partition of n with largest part k
#where the difference of two consecutive parts is NOT in DI
RandPnkD:=proc(n,k,DI) local Die,k1,d,P:

if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN(FAIL):
fi:
if k>n then
 RETURN(FAIL):
fi:
if k=n then
 RETURN([n]):
fi:

Die:=[seq(pnkD(n-k,k1,DI),k1=1..min(n-k,k))]:

for k1 from 1 to min(n-k,k) do
 if k-k1 in DI then
  Die[k1]:=0:
 fi:
od:
d:=LD(Die):
P:=RandPnkD(n-k,d,DI):
[k, op(P)]:
end:


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


#[seq(pnC(n,{1},2),n=1..30)]
#A000009

#[seq(pnC(n,{1,4},5),n=1..30)]
#A003114

#[seq(pnD(n,{0}),n=1..30)]
#A000009

#[seq(pnD(n,{0,1}),n=1..30)]
#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: