#OK to post
#Yuxuan Yang, March 7th, Assignment 13

with(combinat):


pnE:=proc(n) local j,i,nu:
option remember:
if n=0 or n=1 then 
RETURN(1):
fi:
if n<0 then 
RETURN(0):
fi:
nu:=0:
for j from 1 while n-(3*j-1)*j/2>=0 do
nu:=nu-(-1)^j*add(pnE(n-(3*j-1)*j/2)+pnE(n-(3*j+1)*j/2)):
od:
nu:
end:

pnSeqE:=proc(N) local i:
[seq(pnE(i),i=0..N)]:
end:

#restart;read `hw13YuxuanYang.txt`;time(pnSeq(1000));
#39.859
#restart;read `hw13YuxuanYang.txt`;time(pnSeqGF(1000));
#3.281
#restart;read `hw13YuxuanYang.txt`;time(pnSeqE(1000));
#0.031
#restart;read `hw13YuxuanYang.txt`;time(pnSeqGF(3000));
#79.000
#restart;read `hw13YuxuanYang.txt`;time(pnSeqE(3000));
#0.203

Franklin:=proc(P) local P1,k,m,s,i,j:
k:=nops(P):
m:=P[k]:
s:=1:
for i from 1 to k-1 do 
 if P[i]=P[i+1]+1 then 
  s:=s+1:
 else
  break:
 fi:
od:

if s=m and s=k then 
RETURN(FAIL):
fi:

if s=k and s=m-1 then 
RETURN(FAIL):
fi:
P1:=P:
if s<m then
 for j from 1 to s do
  P1[j]:=P[j]-1:
 od:
RETURN([op(P1),s]):
fi:

if s>=m then
RETURN([op(1..m,P)+(1$m),op(m+1..k-1,P)]):
fi:

end:


FranklinOldMaids:=proc(n) local i,allM:
allM:=PnD(n,{0}):
for i from 1 to nops(allM) do
 if Franklin(allM[i])=FAIL then 
  RETURN(allM[i]):
 fi:
od:
end:

#Those two situations in wiki of Pentagonal Number Theorem
#s=k=m or s=k=m-1




















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