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

read "C11.txt";

with(combinat):

# Return a random partition of n with largest part k.
RandPnk := proc(n, k)
    if k > n or k < 0 then
        return FAIL:
    fi:

    if n = 0 then
        return []:
    fi:

    if n = k then
        return [k]:
    fi:

    # What's the next largest part we should use?
    nextLargest := LD([seq(pnk(n - k, j), j=1..min(n-k, k))]):

    return [k, op(RandPnk(n - k, nextLargest))]:
end:

# Return a random partition of n.
RandPn := proc(n)
    # What size should we pick?
    k := LD([seq(pnk(n, k), k=1..n)]):
    RandPnk(n, k):
end:

EstStatAnalOfLargestPart := proc(n, K, N)
    sample := [seq(RandPn(n), k= 1..N)]:

    mean := add(nops(s), s in sample) / N:
    std := sqrt(add((nops(s) - mean)^2, s in sample) / N):

    moms := [mean, std^2]:

    for mom from 3 to K do
        moms := [op(moms),
                    add((nops(s) - mean)^mom, s in sample) / N / std^mom]:
    od:

    moms:
end:

# 3.

# The sum is

#   sum(b(n, k) * t^k * q^n, k, n),

# where b(n, k) = number of partitions of n with k parts.

# This is like a non-exponential version of the hand-enumerator formula that
# Wilf gives. In fact, Wilf also gives this formula! (See near the end of
# chapter 3.) So, yes, I am convinced.

# (In the language of Wilf, this is a "prefab" where each "unlabeled deck"
# consists of a single "unlabeled card" which represents a fixed integer. A
# partition of n is a hand of weight n drawn from this prefab.)

# 4.

pnSeqGFt := proc(N, t)
    # We can compute the coefficients up to q^N using only the first N terms of
    # the product.
    RHS := mul(1 / (1 - t * q^i), i=1..N):
    expr := taylor(RHS, q, N + 1):

    [seq(expand(coeff(expr, q, k)), k=0..N)]:
end:

# 6.

(*

    Using Maple's built-in profiler, I got the following results:

        PnClass
        PnClass := proc(n)
        local n1, L, L1, j;
             |Calls Seconds  Words|
        PROC |   31  15.871 317933838|
           1 |   31   0.001    183| if n = 0 then
           2 |    1   0.000      4|     RETURN({[]})
                                    elif n < 0 then
           3 |    0   0.000      0|     RETURN({})
                                    else
           4 |   30   0.000      0|     L := {};
           5 |   30   0.002      0|     for n1 to n do
           6 |  465   0.001   1902|         L1 := PnClass(n-n1);
           7 |  465   0.159      0|         for j to nops(L1) do
           8 |109350  15.708 317931689|             L := L union {sort([op(L1[j]), n1],`>`)}
                                            end do
                                        end do;
           9 |   30   0.000     60|     RETURN(L)
                                    end if
        end proc

        Pn
        Pn := proc(n)
        local k;
             |Calls Seconds  Words|
        PROC |    1   0.015 265176|
           1 |    1   0.015 265176| [seq(op(Pnk(n,n-k+1)),k = 1 .. n)]
        end proc

    What does this mean? The overwhelming amount of time in PnClass is spent
    sorting the partitions and creating the set.

    Maybe the problem is that we're sorting the partitions? Well, it's hard to
    tell, because we're also inserting them into a set. In order to not get
    repetitions, we *have* to sort all of them.

*)