Help:=proc()
print(`Solid partition functions by Brian Nakamura`):
print(`The primary functions that this file contains:`):
print(`PPK(M,N,K,q), Solid22(N,M,q), Solid33(M,N,q), SolidKK(M,N,K,q)`):
end:



# PPK(M,N,K,q): generalization of PP2 and PP3 to K rows

PPK:=proc(M,N,K,q)
normal(PPKa([N$K],[],M+1,q)/q^(K*N)):
end:



# PPKa(L1,L2,N,q): similar to PP2a and PP3a, except...
# inputs lists L1=[a1, a2, ..., ai] and L2=[b1, b2, ..., bj] and pos
# integer N, and behaves as a dynamic for loop (with changing levels
# of recursion) for PPK.
# In particular, the initial call is PPK([a1, a2, .., ak], [], N, q).
# The function then does:
# PPK([a1, a2, .., ak], [], N, q) -> PPK([a2, .., ak], [b1], N, q)
#   -> PPK([a3, .., ak], [b1, b2], N, q) -> ..., where b1 is an index
#   between 0 and a1, b2 is an index between 0 and min(b1,a2), etc.
# Once it hits PPKa([], [b1, .., bk], N, q), it reacts as if it hit
# the K-th level of recursion and loops to PPKa([b1,..,bk], [], N-1, q).
# 
# For example, PP2a(a1,a2,N,q) = PPKa([a1,a2],[],N,q)
# and PP3a(a1,a2,a3,N,q) = PPKa([a1,a2,a3],[],N,q)

PPKa:=proc(L1, L2, N, q)
local k,min1,Lsum,fq:
option remember:


if nops(L1)=0 then
  Lsum:=0:
else
  Lsum:=add(L1[k],k=1..nops(L1)):
fi:


if N=1 then
  RETURN(q^Lsum):
fi:

if nops(L1)=0 then
  RETURN(PPKa(L2,[],N-1,q)):
fi:

if nops(L2)=0 then
  min1:=L1[1]
else
  min1:=min(L1[1],L2[nops(L2)]):
fi:


fq:=0:

for k from 0 to min1 do
  fq:=fq+PPKa(L1[2..nops(L1)], [op(L2), k], N, q):
od:

if nops(L2)=0 then
  fq:=fq*q^Lsum:
fi:

sort(expand(fq)):

end:



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


# Solid22(N,M,q): inputs pos integers M and N, and outputs the generating
# function (in q) for all solid partitions whose 4D Ferrers graph fits
# inside a 2x2xMxN (four-dimensional) box.

Solid22:=proc(N,M,q) 

normal(Solid22a(N,N,N,N,M+1,q)/q^(4*N)):
end:



# Solid22a(a11,a12,a21,a22,M,q): computes the generating function for top
# layer of a11,a12,a21,a22

Solid22a:=proc(a11,a12,a21,a22,M,q)
local fq, b11,b12,b21,b22:
option remember:

if M=1 then
  RETURN(q^(a11+a12+a21+a22)):
fi:

fq:=0:

for b11 from 0 to a11 do

  for b12 from 0 to min(a12,b11) do
    for b21 from 0 to min(a21,b11) do
      for b22 from 0 to min(a22,b12,b21) do
        fq:=fq+Solid22a(b11,b12,b21,b22,M-1,q):
      od:
    od:
  od:

od:

fq:=fq*q^(a11+a12+a21+a22):

sort(expand(fq)):

end:



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


# The equivalent of Solid22 except with a 3x3xMxN box...
# The coefficients of the generating function can also be
# interpreted as such: the coeff of q^t is the number of 
# ways t can be partitioned in a 3xMxN box such that each
# component is at most 3

Solid33:=proc(M,N,q)
normal(Solid33a([ [N,N,N], [N,N,N], [N,N,N] ], M+1, q)/q^(9*N)):
end:



# Computes Solid33 when the top (3x3) layer is already set

# Format of layer1=
#    [ [a11,a12,a13], [a21,a22,a23], [a31,a32,a33] ]

Solid33a:=proc(layer1, M, q)
local fq, Lsum, i, j, b11,b12,b13,b21,b22,b23,b31,b32,b33:
option remember:

Lsum:=0:

for i from 1 to 3 do
  for j from 1 to 3 do
    Lsum:=Lsum + layer1[i][j]:
  od:
od:

if M=1 then
  RETURN(q^Lsum):
fi:

fq:=0:

for b11 from 0 to layer1[1][1] do
for b12 from 0 to min(layer1[1][2], b11) do
for b13 from 0 to min(layer1[1][3], b12) do
for b21 from 0 to min(layer1[2][1], b11) do
for b22 from 0 to min(layer1[2][2], b12, b21) do
for b23 from 0 to min(layer1[2][3], b13, b22) do
for b31 from 0 to min(layer1[3][1], b21) do
for b32 from 0 to min(layer1[3][2], b22, b31) do
for b33 from 0 to min(layer1[3][3], b23, b32) do
  fq:=fq + Solid33a([ [b11,b12,b13], [b21,b22,b23], [b31,b32,b33] ], M-1, q):
od:od:od:od:od:od:od:od:od:

fq:=fq*q^Lsum:

sort(expand(fq)):

end:




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


# SolidKK(M,N,K,q): generalizes Solid22 and Solid33 and inputs pos integers
# M, N, and K and outputs the generating function for all solid partitions
# whose 4D Ferrers diagram fits inside a K x K x M x N box.
# Likewise, the coeff of q^t (in the generating function) will be the number
# of solid partitions in a KxMxN box with components not more than K.
SolidKK:=proc(M,N,K,q)
normal(SolidKKa([[N$K]$K], N*K^2, [], 1, 1, M+1, q)/q^(N*K^2)):
end:




# SolidKKa(layer1, Lsum, layer2, i, j, M, q): this is the analogue of the
# above PPk function. It produces the generating function when the top 
# KxK layer (layer1) is already set. The value Lsum is merely the sum over
# all values in layer1 and is used for (possible) efficiency purposes.
# The layer2 is the "next step" in the recursive call that gets built up,
# and indices i and j represent the next value to add to layer2.

# Format examples:
# layer1 = [ [a11, .. ,a1k], [a21, .. ,a2k], ... [ak1, .. ,akk] ]
# Lsum = add(add(layer1[s][t], s=1..k), t=1..k)

SolidKKa:=proc(layer1, Lsum, layer2, i, j, M, q)
local k,s,t,fq,inew,jnew,L2len,L2sum:
option remember:


if layer2=[] and M=1 then
  RETURN(q^Lsum):
fi:

k:=nops(layer1):

fq:=0:


# these cases are only used when not both i and j are k
if j=k then
  inew:=i+1:
  jnew:=1:
else
  inew:=i:
  jnew:=j+1:
fi:


if i=k and j=k then

  L2sum:=add(add(layer2[s][t], s=1..k-1), t=1..k) + add(layer2[k][t], t=1..k-1):

  for t from 0 to min(layer1[i][j], layer2[i-1][j], layer2[i][j-1]) do
    fq:=fq + SolidKKa([ op(1..k-1,layer2), [op(layer2[k]),t] ], L2sum+t, [], 1, 1, M-1, q):
  od:

elif i=1 and j=1 then
  for t from 0 to layer1[1][1] do
    fq:=fq + SolidKKa(layer1, Lsum, [[t]], 1, 2, M, q):
  od:
elif i=1 then
  for t from 0 to min(layer1[1][j], layer2[1][j-1]) do
    fq:=fq + SolidKKa(layer1, Lsum, [[op(op(layer2)),t]], inew, jnew, M, q):
  od:  
elif j=1 then
  for t from 0 to min(layer1[i][1], layer2[i-1][j]) do
    fq:=fq + SolidKKa(layer1, Lsum, [op(layer2), [t]], i, j+1, M, q):
  od:
else
  L2len:=nops(layer2):
  for t from 0 to min(layer1[i][j], layer2[i-1][j], layer2[i][j-1]) do
    fq:=fq + SolidKKa(layer1, Lsum, [ op(1..L2len-1,layer2), [op(layer2[L2len]),t] ], inew, jnew, M, q):
  od:
fi:


if layer2=[] then
  fq:=fq*q^Lsum:
fi:

sort(expand(fq)):


end:



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

# NOTES:

# Unfortunately, while SolidKK appears to function correctly, it is too slow
# (and the approach probably too naive) to produce substantial data. It has been
# verified to coincide with Solid22 and Solid33 for small values. In addition,
# the values generated from these were checked against Sloane to verify their correctness.

# p2:=Solid22(2,2,q): subs(q=1,p2);
#    yields 168
# p3:=Solid33(3,3,q): subs(q=1,p3);
#    yields 17792748

# and checked versus Sloane, we get the following sequence:
#    A007760 -  Solid partitions of height n in cube of side n.