#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 5
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

# Stuff from class 5:

Help := proc(): print(`Di(L), GP(L, x), GQP1(L, x, m), GQP(L, x), GFtoSeq(f,q,m), pmn(m,n), F(m, x), T(n)`): end:

Di := proc(L):
    [seq(L[i] - L[i-1], i=2..nops(L))];
end:

GP := proc(L, x) local i, d, coeffs, diff:
    d := nops(L)-1:
    diff := L[1..d+1]:
    coeffs := []:
    if convert(diff, set) = {0} then:
        return 0:
    fi:
    for i from 0 to d do:
        coeffs := [op(coeffs), diff[1]]:
        diff := Di(diff):
        if convert(diff, set) = {0} then:
            return factor(expand(add(coeffs[j]*binomial(x, j-1), j=1..i+1))):
        fi:
    od:
    return FAIL:
end:

GQP1 := proc(L, x, m) local i, k, seqs:
    seqs := [seq([seq(L[i+k*m], k=0..floor((nops(L)-i)/m))], i=1..m)]:
    return [seq(GP(seqs[i], x), i=2..nops(seqs)), GP(seqs[1], x)]:
    # Our convention for quasipolynomials requires what probably
    # should be the first term to be at the end instead. Yet another
    # issue caused by the fact that lists do not start at 0.
end:

GQP := proc(L, x) local m, qp, rv:
    for m from 1 to nops(L/2) do
        qp := GQP1(L, x, m):
        if not FAIL in qp then return qp: fi:
    od:
    return FAIL:
end:

#############
# Problem 2 #
#############

# Remark 1: There is no sensible reason to start at 1. So I started at
# 0 instead. This will come in very handy in problems 6 and 7.

GFtoSeq := proc(f,q,m):
    [seq(coeff(taylor(f, q, m+1), q, i), i=0..m)]:
end:

#############
# Problem 3 #
#############

pmn := proc(m,n):
    GQP(GFtoSeq(product(1/(1-q^i), i=1..m), q, (m+4)*m!), n):
end:

# We have:
# pmn(2, n) = [n+1, n+1]
# pmn(3, n) = [(1+n)(3n+1), (1+n)(3n+2), 3(1+n)^2, (3n+4)(1+n),
# (3n+5)(1+n), 1+3n+3n^2].

# The next two are too long to include here.

#############
# Problem 4 #
#############

# Human mathematics is good enough to figure out this problem:

F := proc(m,x):
    [seq(x-i/m, i=1..m-1), x/m]:
end:

#############
# Problem 5 #
#############

T := proc(n):
    return round(n^2/12) - floor(n/4)*floor((n+2)/4):
end:

# Now we just run GQP([seq(T(i), i=0..48)], x) and get
# the quasipolynomial
# [x(2+3x)    , x(1+3x)    , 1+3x+3x^2,
#  x(2+3x)    , (1+x)(1+3x), 1+3x+3x^2,
#  (1+x)(2+3x), (1+x)(1+3x), 3(1+x)^2 ,
#  (1+x)(2+3x), (3x+4)(1+x), 3x^2      ].

#############
# Problem 6 #
#############

# Remark: The statement of the problem seems to be wrong. It should
# probably say a >= b >= c > 0, and b + c > a (and similarly in the
# next two statements of similar problems). Otherwise
# degenerate triangles/quadrilaterals are included, which is probably
# not what we want.

# In any event, suppose a >= b >= c >= d > 0, and b + c + d > a.

# We can write the generating function for these partitions as
# follows:

# \sum_{d > 0} \sum_{c >= d} \sum_{b >= c} \sum_{b+c+d >= a >=b} q^{a+b+c+d}.

# Maple can simplify this if we run

# simplify(sum(sum(sum(sum(q^(a+b+c+d), a=b..b+c+d-1), b=c..infinity),
# c=d..infinity), d=1..infinity));

# We end up with the generating function
# (q^3 - q^2 + 1)q^4 / ((q^4-1)(q^2-1)(q-1)(q^6-1))

# The form of the denominator indicates that the coefficients can be
# given by a quasipolynomial of order lcm(1,2,4,6) = 12. The degree of the
# quasipolynomial should be at most 4-1 = 3 (since there are 4 factors
# in the denominator). Hence we only need to check the first 49
# terms. Taking 200 or so just to be safe, we get the following
# quasipolynomial
# (by running GQP(GFtoSeq((q^3-q^2+1)q^4/((q^4-1)(q^2-1)(q-1)(q^6-1))))):

# AnswerToProblem6 := [1/2*x*(1+9*x+12*x^2), 3/2*x*(1+3*x+4*x^2),
#                      1/2*x*(5+15*x+12*x^2), 1+7/2*x+15/2*x^2+6*x^3,
#                      1/2*(1+x)*(12*x^2+9*x+2), 1+13/2*x+21/2*x^2+6*x^3,
#                      1/2*(1+x)*(12*x^2+15*x+4), 3/2*(1+x)*(4*x^2+5*x+2),
#                      1/2*(1+x)*(12*x^2+21*x+8), 1/2*(1+x)*(12*x^2+21*x+10),
#                      1/2*(1+x)*(12*x^2+27*x+14), 1/2*x*(1+3*x+12*x^2)]:

# That is a mess.
# If we really do want to allow degenerate quadrilaterals, then (by a
# very similar process) we instead get the quasipolynomial

# AnswerToProblem6D := [1/2*x*(5+15*x+12*x^2), 1+13/2*x+21/2*x^2+6*x^3,
#                       1/2*(1+x)*(12*x^2+9*x+2), 3/2*(1+x)*(4*x^2+5*x+2),
#                       1/2*(1+x)*(12*x^2+15*x+4), 1/2*(1+x)*(12*x^2+21*x+10),
#                       1/2*(1+x)*(12*x^2+21*x+8), 1/2*(1+x)*(12*x^2+27*x+16),
#                       1/2*(1+x)*(12*x^2+27*x+14), 3/2*(1+x)*(4*x^2+11*x+8),
#                       1/2*(1+x)*(12*x^2+33*x+22), 1+7/2*x+15/2*x^2+6*x^3]:

#############
# Problem 7 #
#############

# Now we would like to do the same thing we did in problem 6, but with
# one additional term. This is going to take more computer time but
# should otherwise be the same process.

# The generating function is produced by:

# simplify(sum(sum(sum(sum(sum(q^(a+b+c+d), a=b..b+c+d+e-1), b=c..infinity),
# c=d..infinity), d=e..infinity), e=1..infinity));

# which yields
# -(q^10+q^9+q^8+q^7+q^6+q^5+q^4+q^3+q^2+q+1)*q^5/(q^6-1)/(q^8-1)/(q-1)^3/(q+1)/(q^3+q^2+q+1)/(q^4+q^3+q^2+q+1)

# Massaging the denominator a little, we can put it in the form
#           (q^6-1)(q^8-1)(q^2-1)(q^4-1)(q^5-1)
# so we are looking for a quasipolynomial of order lcm(6,8,2,4,5) = 120
# and degree at most 4, and that will require at least 601 terms.

# I don't have the patience to nicely format an order-120
# quasipolynomial so I have put it in a separate file, hw5p7.txt.

# Note also that if we allow degenerate pentagons, then the whole
# process has to be done again. This time the generating function will
# be

# simplify(sum(sum(sum(sum(sum(q^(a+b+c+d), a=b..b+c+d+e), b=c..infinity),
# c=d..infinity), d=e..infinity), e=0..infinity))

# The resulting quasipolynomial is, of course, different, but the
# method is the same.