#!/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.