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

with(combinat):

Help:=proc():
print(`BZ(L,j),BZa(L,j),Glaisher(L),InvGlaisher(L),Sylv(L),InvSylv(L)`):
end:

# 1.

#BZ(L,j): inputs a partition L and an integer j, and outputs a pair [L', j'] where j' has
#opposite parity than j and |L|+(3*j-1)*j/2 =|L'|+(3*j'-1)*j'/2
BZ:=proc(L,j) local L1,j1,n,t,i:
n:=add(op(L))+j*(3*j-1)/2:
t:=nops(L):
if t+3*j<L[1] then
  j1:=j+1:
  L1:=[seq(L[i]+1,i=2..t),1$(L[1]-3*j-t-1)]:
else
  j1:=j-1:
  L1:=[t+3*j-1,seq(L[i]-1,i=1..t)]:
fi:
[L1,j1]:
end:

#BZa(L,j): inputs a partition L and an integer j, and outputs a pair [L', j'] where j' has
#opposite parity than j and |L|+(3*j+1)*j/2 =|L'|+(3*j'+1)*j'/2
#NOTE: this implentation matches the Bressoud-Zeilberger paper from 1985. The only difference
#is in the calculation of n (i.e. whether we use a(j)=j*(3*j-1)/2 or a(j)=j*(3*j+1)/2)
#The results seem to be identical in either case.
BZa:=proc(L,j) local L1,j1,n,t,i:
n:=add(op(L))+j*(3*j+1)/2:
t:=nops(L):
if t+3*j<L[1] then
  j1:=j+1:
  L1:=[seq(L[i]+1,i=2..t),1$(L[1]-3*j-t-1)]:
else
  j1:=j-1:
  L1:=[t+3*j-1,seq(L[i]-1,i=1..t)]:
fi:
[L1,j1]:
end:


# 2.

#Glaisher(L): input a partition into odd parts, and output a partition into distinct parts
#using Glaisher's bijection
Glaisher:=proc(L) local L1,pos,num,cur,pwr:
if nops(L)=0 then RETURN([]): fi:
L1:={}:
pos:=1:
for cur from L[1] to 1 by -2 do
  for num from 0 while (pos<=nops(L) and L[pos]=cur) do pos:=pos+1: od:
  while num>0 do:
    pwr:=2^trunc(log[2](num)):
    num:=num-pwr:
    L1:=L1 union {pwr*cur}:
  od:
od:
sort(convert(L1,list),`>`):
end:

#InvGlaisher(L): input a partition into distinct parts, and output a partition into odd parts
#using Glaisher's bijection
InvGlaisher:=proc(L) local L1,pos,cur,num:
if nops(L)=0 then RETURN([]): fi:
L1:=[]:
for pos from 1 to nops(L) do
  num:=1:
  cur:=L[pos]:
  while cur mod 2=0 do:
    cur:=cur/2:
    num:=num*2:
  od:
  L1:=[op(L1),cur$num]:
od:
sort(L1,`>`):
end:

# 3.

#Sylv(L): input a partition into odd parts, and output a partition into distinct parts
#using Sylvester's bijection
Sylv:=proc(L) local n,i,r,m,L1,L2:
n:=nops(L):
if n=0 then RETURN([]): fi:
if max(L)=1 then RETURN ([n]): fi:
m:=0:
r:=1: #Because L[1] is known to be >1
L1:=[]:
for i from 2 to n do
  if L[i]=1 then
    m:=m+1:
  else
    r:=r+1:
    L1:=[op(L1),L[i]-2]:
  fi:
od:
L2:=Sylv(L1):
[(L[1]-1)/2+r+m,(L[1]-1)/2+r-1,op(L2)]:
end:

#InvSylv(L): input a partition into distinct parts, and output a partition into odd parts
#using Sylvester's bijection
InvSylv:=proc(L) local a1,n,i,r,m,L1:
n:=nops(L):
if n=0 then RETURN([]): fi:
if n=1 then RETURN ([1$(L[1])]): fi:
if n=2 then
  L1:=[]:
else
  L1:=InvSylv([op(3..n,L)]):
fi:
r:=nops(L1)+1:
m:=L[1]-L[2]-1:
a1:=L[2]-r+1:
[2*a1+1,seq(L1[i]+2,i=1..nops(L1)),1$m]:
end: