#OK to post homework
#Blair Seidler, 2021-03-07 Assignment 12

with(combinat):

Help:=proc():
print(`PnkC(n,k,S,m), PnC(n,S,m), RandPnkC(n,k,S,m), RandPnC(n,S,m), pnkC(n,k,S,m), pnC(n,S,m)`):
print(`PnkD(n,k,DI), PnD(n,DI), RandPnkD(n,k,DI), RandPnD(n,S,DI), pnkC(n,k,DI), pnD(n,DI)`):
end:

# 1.

#PnkC(n,k,S,m): 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 k>n or not (k mod 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
 if k1 mod m in S then
   L1:=PnkC(n-k,k1,S,m):
   L:=[op(L), seq([k, op(L1[j])],j=1..nops(L1))]:
  fi:
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):
pnkC:=proc(n,k,S,m) local k1,p:
option remember:
if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN(0):
fi:
if k>n or not (k mod m in S) then
 RETURN(0):
fi:
if n=k then
 RETURN(1):
fi:
p:=0:
for k1 from min(n-k,k) by -1 to 1 do
 if k1 mod m in S then
   p:=p+pnkC(n-k,k1,S,m):
  fi:
od:
p:
end:

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 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
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 all parts are congruent to a member of S mod m
PnD:=proc(n,DI) local k:option remember:[seq(op(PnkD(n,n-k+1,DI)),k=1..n)]:end:

#pnkD(n,k,DI):
pnkD:=proc(n,k,DI) local k1,p:
option remember:
if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN(0):
fi:
if k>n then
 RETURN(0):
fi:
if n=k then
 RETURN(1):
fi:
p:=0:
for k1 from min(n-k,k) by -1 to 1 do
 if not(k-k1 in DI) then
   p:=p+pnkD(n-k,k1,DI):
  fi:
od:
p:
end:

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 all parts are congruent to a member of S mod m.
RandPnkD:=proc(n,k,DI) local die1, k1,k2,i:
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))]:
for i from 1 to min(n-k,k) do
  if k-i in DI then
    die1[i]:=0:
  fi:
od:
k2:=LD(die1):
[k, op(RandPnkD(n-k,k2,DI))]:

end:

#RandPnD(n,DI): A random partition of n
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:

# 3.

#[seq(pnC(n,{1},2),n=1..30)] yields
#[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: Expansion of Product_{m >= 1} (1 + x^m);
# number of partitions of n into distinct parts; number of partitions of n into odd parts.

#[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: Number of partitions of n into parts 5k+1 or 5k+4.

#[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: Expansion of Product_{m >= 1} (1 + x^m);
# number of partitions of n into distinct parts; number of partitions of n into odd parts.
#Once again, this is comforting because the number of partitions into odd parts and the
#number of partitions into distinct parts are supposed to be the same.

#[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: Number of partitions of n into parts 5k+1 or 5k+4.

#Yes, I see that this is the same as [seq(pnC(n,{1,4},5),n=1..30)]. I haven't figured out why yet.

#### From C4.txt ####

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