# OK to post homework
# Robert Dougherty-Bliss
# 2021-02-14
# Assignment 10

read "C9.txt":

with(combinat):

permGF := proc(n, f, x)
    gf := 0:
    for p in permute(n) do
        gf := gf + x^(f(p)):
    od:

    gf:
end:

invs := proc(p)
    n := nops(p):
    total := 0:
    for pair in choose(n, 2) do
        if p[max(pair)] < p[min(pair)] then
            total := total + 1:
        fi:
    od:

    total:
end:

ExactGFinv := proc(n, x)
    permGF(n, invs, x):
end:

maj := proc(p)
    n := nops(p):
    total := 0:
    for i from 1 to n - 1 do
        if p[i] > p[i + 1] then
            total := total + i:
        fi:
    od:

    total:
end:

ExactGFmaj := proc(n, x)
    permGF(n, maj, x):
end:

desc := proc(p)
    n := nops(p):
    total := 0:
    for i from 1 to n - 1 do
        if p[i] > p[i + 1] then
            total := total + 1:
        fi:
    od:

    total:
end:

ExactGFdes := proc(n, x)
    permGF(n, desc, x):
end:

quots := proc(l)
    res := []:
    for k from 1 to nops(l) - 1 do
        res := [op(res), l[k + 1] / l[k]]:
    od:

    res:
end:

# The nth inversion generating function is the product of the first n partial
# sum of the geometric series sum(x^k, k >= 0).

# This is easy to see experimentally: Take adjacent quotients and simplify.

# This is also not hard to see combinatorially: Inserting n into position j of
# a permutation on [n - 1] gives n - j new inversions plus the old ones, which
# gives a nice recurrence.

# It feels like there should be some nice generatingfunctionological proof.
# The idea here would be to interpret each partial sum as the generating
# function for a smaller piece. For example, 1 + x + x^2 tells us how many
# inversions we get from placing 3 into the permutation [1, 2] based on the
# position. Somehow we can "stitch" this information together to get the
# product formula, but I would have to think about how exactly this works.

# The major index generating function seems exactly the same!
# This implies the conjecture "number of perms on [n] with k inversions" equals
# "number of perms on [n] with major index k."

# 2.

EstGFperm := proc(n, x, f, K)
    pgf := 0:
    for k from 1 to K do
        p := randperm(n):
        pgf := pgf + x^(f(p)):
    od:

    pgf / K:
end:

EstGFinv := proc(n, x, K)
    EstGFperm(n, x, invs, K):
end:

EstGFmaj := proc(n, x, K)
    EstGFperm(n, x, maj, K):
end:

EstGFdes := proc(n, x, K)
    EstGFperm(n, x, desc, K):
end:

# 3.

# Ah, now I see!

# The given formula is actually just what Wilf calls "the" exponential
# generating function theorem in his chapter on decks and hands. (See
# generatingfunctionology, chapter 3.) This is nothing more than a way to
# compactly write the two-variable hand-enumerator which is, of course, also a
# sum of "weighted polynomials" like we have here.

# So, yes, famously, the hand-enumerator for partition deck is exp(t (exp(x) -
# 1)), and setting t = 1 gives the standard egf for the bell numbers.

PnSeq := proc(n, t)
    n! * coeff(taylor(exp(t * (exp(x) - 1)), x, n + 1), x, n):
end:

PnSeqMean := proc(n)
    f := PnSeq(n, t):
    subs(t = 1, diff(f, t) / f):
end:

# It looks like the answer here is "yes."
# Can we figure this out symbolically?

# We know that P(n, t) = n! * [x^n] exp(t * (exp(x) - 1)).
# So, the derivative w.r.t. t is (probably)
#       n! [x^n] (exp(x) - 1) * exp(t * (exp(x) - 1)).
# Letting t = 1 and gives

#       n! [x^n] (exp(x) - 1) * exp(exp(x) - 1).

# Normalizing gives

#       n! [x^n] (exp(x) - 1) * exp(exp(x) - 1)
#   /   n! [x^n]                exp(exp(x) - 1),

# or

#       ([x^n] exp(x) exp(exp(x) - 1) - [x^n] exp(exp(x) - 1)) / [x^n] exp(exp(x) - 1)
#     = ((nth binomial transform of bell numbers) - bell(n)) / bell(n).