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

# 1.

# Given a partition L of add(L) - a(j), where a(j) = (3j^2 + j) / 2, use the
# Bressoud-Zeilberger involution to generate a partition of add(L) - a(j'),
# where j' has the opposite parity of j.
BZ := proc(L, j)
    local t:
    t := nops(L):

    if t + 3 * j >= L[1] then
        [t + 3 * j - 1, seq(l - 1, l in L)]:
    else
        [seq(l + 1, l in L[2..]), seq(1, k=1..L[1] - 3 * j - t - 1)]:
    fi:
end:

# If t + 3j >= L[1], then BZ(L, j) takes L to a partition of
#       n + 3j - 1 - a(j) = n - a(j - 1).

# If t + 3j < L[1], then BZ(L, j) takes L to a partition of
#       n - 3j - 2 - a(j) = n - a(j + 1).

# 2.

Glaisher := proc(L)
    local counts, current, k, res:

    current := 1:
    counts := [[1, current]]:
    for k from 2 to nops(L) do
        if L[k] = L[current] then
            counts[current][1] := counts[current][1] + 1:
        else
            current := k:
            counts := [op(counts), [1, current]]:
        fi:
    od:

    res := []:

    for count in counts do
        x := L[count[2]]:
        freq := count[1]:
        binary := convert(freq, base, 2):
        usedDigits := select(p -> p[1] = 1, [seq([binary[k], k-1], k=1..nops(binary))]):
        res := [op(res), seq(2^d[2] * x, d in usedDigits)]:
    od:

    sort(res):
end:

# Experimenting gave me this conjecture:
#       nops(Glaisher([d $ n])) = A120(n) = number of 1's in binary expansion of n
# This is obvious in hindsight - but cool!

# What's something that would depend on d?

# Let a(d, n) be the largest part of Glaisher([d $ n]). This is clearly "most
# significant bit of n" * d. In other words,
#       a(d, n) = d * 2^(floor(lg(n))).

# In the OEIS?
# a(1, n) - yes, definition
# a(2, n) - yes, probably easy to see
# a(3, n) - yes, comment
# a(4, n) - no!
# a(n, n) - no!

# None of these entries mention the Glaisher map.

InvGlaisher := proc(L)
    local repeats, part, power, odd, count, res, p:
    # Each part is 2^m * odd.
    # Have "odd" repeat itself 2^(m_1) + ... 2^(m_n) times, where n is the
    # number of times odd appears multiplied by some power of 2.
    repeats := table():
    for part in L do
        if part mod 2 = 1 then
            power := 1:
        else
            power := 2^ifactors(part)[2][1][2]:
        fi:

        odd := part / power:
        if not assigned(repeats[odd]) then
            repeats[odd] := 0:
        fi:

        repeats[odd] := repeats[odd] + power:
    od:

    res := []:
    for p in entries(repeats, pairs) do
        odd := lhs(p):
        count := rhs(p):
        res := [op(res), odd $ count]:
    od:

    sort(res):
end:

# 3.

T := proc(L)
    if L = [] then
        return []
    fi:

    if max(L) = 1 then
        return [add(L)]:
    fi:

    m := 0:
    # Count the number of 1's.
    while L[-(m + 1)] = 1 do m := m + 1: od:

    prev := T([seq(x - 2, x in L[2..-(m + 1)])]):

    # r + m = nops(T)
    a := floor(L[1] / 2):
    [a + nops(L), a + nops(L) - m - 1, op(prev)]:
end:

S := proc(L)
    local prev, r, a1, m:
    if L = [] then
        return []:
    fi:

    if nops(L) = 1 then
        return [1 $ L[1]]:
    fi:

    prev := S(L[3..]):
    r := nops(prev) + 1:

    a1 := L[2] - r + 1:
    m := L[1] - L[2] - 1:

    [2 * a1 + 1, seq(x + 2, x in prev), 1 $ m]:
end: