# OK to post homework
# Robert Dougherty-Bliss
# 2021-04-12
# Assignment 22

read "C22.txt":

# Tally:

#   Attempted challenge 1 ($5)
#   Attempted challenge 2 ($5)
#   Attempted challenge 3 ($10)
#   Attempted challenge 5 ($15)

# In honor of how difficult it was to solve these problems without it, I hereby
# donate any potential rewards to the OEIS.

# 4.

#####
# FIRST PROBLEM: Guess syt([a, b]).
# Answer: binomial(a + b, b) * (a - b + 1) / (a + 1)

# Equivalent: Guess a(n, k) = syt([n + k, n]).

A := (n, N) -> seq(syt([n + k, n]), k=0..N-1):

# These are definitely polynomials. Let's ask Maple to guess them.

with(gfun):
for n from 1 to 5 do
    factor(rsolve(listtorec([A(n, 20)], a(k))[1], a(k)));
od;

# Looks like the kth column is binomial(2n + k, n) (k + 1) / (n + k + 1).

# Thus, our guess is

#       syt([n + k, n]) = binomial(2n + k, n) (k + 1) / (n + k + 1),

# or

#       syt([a, b]) = binomial(a + b, b) (a - b + 1) / (a + 1)

# If we take a = b, then we get the Catalan numbers. This must be right!

(*

    Proof?

    Let a(n, k) = syt([n + k, n]).

    Recurrence:

        syt([n + k, n]) = [k > 0] syt([n + k - 1, n]) + [n > 1] syt(n + k, n - 1)

        a(n, k) = [k > 0] a(n, k - 1) + [n > 1] a(n - 1, k + 1)

    Base cases:

        a(n, 0) = C(n) (Catalan numbers. Easy bijection with Dyck paths!)
        a(1, k) = k + 1

    This recurrence + base cases uniquely determine the sequence.

*)

# Check the recurrence:

guess := (n, k) -> binomial(2 * n + k, n) * (k + 1) / (n + k + 1):

R := (guess(n, k - 1) + guess(n - 1, k + 1)) / guess(n, k):

simplify(convert(R, factorial));

# Base cases:

simplify(guess(n, 0) / binomial(2 * n, n) * (n + 1));
simplify(guess(1, k) / (k + 1));

# QED!

#### SECOND PROBLEM
# Guess syt([a, b, c]).

# Answer: (b - c + 1)(a - b + 1)(a - c + 2) choose(a + b + c, a, b, c) / (a + 2) / (a + 1) / (b + 1)

# Or, in terms of syt([n + i + j, n + i, n]),

#   (i + 1) (j + 1) (j + i + 2) choose(3n + 2i + j, n, n + i, n + i + j) / (n + i + j + 2) / (n + i + j + 1) / (n + i + 1)

#### BEGIN EXPLORATION

B := (i, j, N) -> seq(syt([n + i + j, n + i, n]), n=1..N-1):

# Wild guess: For fixed i, the columns of B'(j, n) are polynomials. Let's try
# to guess.

denoms := []:

maxI := 2:

for i from 0 to maxI do
    denoms := [op(denoms), []]:

    BjCols := (n, N) -> seq(syt([n + i + j, n + i, n]), j=0..N-1):
    print(`i is`, i);
    for n from 1 to 10 do
        print(`n is`, n);
        print(BjCols(n, 5));
        soln := factor(rsolve(listtorec([BjCols(n, 30)], a(j))[1], a(j))):
        soln := simplify(convert(soln, factorial)):
        denoms[-1] := [op(denoms[-1]), denom(soln)]:
        print(soln):
    od:
od:

# Guess denominators:
for i from 0 to maxI do
    rsolve(listtorec(denoms[i + 1], a(u))[1], a(u));
od;

(*

    Let b(i, j, n) = sym([n + i + j, n + i, n]).

    Guesses:

        b(0, j, n) = (j + 1) (j + 2) (j + 3n)! / (j + n + 2))! / (n! * (n + 1)!)
        b(1, j, n) = (j + 1) (j + 3) (j + 2 + 3n)! / (j + 3 + n)! / (n! * (n + 2)! / 2)
        b(2, j, n) = (j + 1) (j + 4) (j + 4 + 3n)! / (j + 4 + n)! / (n! * (n + 3)! / 3)

    Pattern seems pretty clear!

        b(i, j, n) = (j + 1) * (j + i + 2) * (j + 2 * i + 3 * n)! / ((j + i + 2 + n)! * (n! * (n + i + 1)! / (i + 1))):

    We can definitely clean this up.

        b(i, j, n) = (j + 1) * (j + i + 2) * (j + 2 * i + 3 * n)! / ((j + i + 2 + n)! * (n! * (n + i + 1)! / (i + 1))):
                   = (i + 1) (j + 1) (j + i + 2) choose(3n + 2i + j, n, n + i, n + i + j) / (n + i + j + 2) / (n + i + j + 1) / (n + i + 1)

    Translating:
        n = c
        i = b - c
        j = a - b

            (b - c + 1)(a - b + 1) choose(a + b + c, a, b, c) / (a + 2) / (a + 1) / (b + 1)

*)

with(combinat):
guess := (i, j, n) -> (j + 1) * (j + i + 2) * (j + 2 * i + 3 * n)! / ((j + i + 2 + n)! * (n! * (n + i + 1)! / (i + 1))):
guess2 := (i, j, n) -> (i + 1) * (j + 1) * (j + i + 2) * multinomial(3 * n + 2 * i + j, n, n + i, n + i + j) / (j + i + n + 2) / (n + i + j + 1) / (n + i + 1):

# Try some random combinations.

ff := (x, n) -> expand(pochhammer(x - n + 1, n)):

for i from 0 to 4 do
    for j from 0 to 5 do
        for n from 1 to 10 do
            actual := syt([n + i + j, n + i, n]):
            if actual <> guess(i, j, n) then
                print(i, j, n);
            fi:
            if guess(i, j, n) <> guess2(i, j, n) then
                print(i, j, n);
            fi:
        od:
    od:
od:

# END EXPLORATION

#### THIRD PROBLEM
# Guess syt([a(1), a(2), ..., a(n)]).

# Previous two answers:

#   (a - b + 1) binomial(a + b, b) / (a + 1)
#   (b - c + 1)(a - b + 1)(a - c + 2) choose(a + b + c, a, b, c) / (a + 2) / (a + 1) / (b + 1)

Guess2 := (a, b) -> (a - b + 1) * binomial(a + b, b) / (a + 1):
Guess3 := (a, b, c) -> (a - b + 1) * (a - c + 2) * (b - c + 1) * multinomial(a + b + c, a, b, c) / (a + 2) / (a + 1) / (b + 1):

# Let's guess just ONE more first: syt([a, b, c, d]).

Guess4 := (a, b, c, d) -> (a - b + 1) * (a - c + 2) * (a - d + 3) * (b - c + 1) * (b - d + 2) * (c - d + 1) * multinomial(a + b + c + d, a, b, c, d) / (a + 3) / (a + 2) / (a + 1) / (b + 2) / (b + 1) / (c + 1):

# This looks correct!

# I believe that the formula is as follows:

#     (a(1) - a(2) + 1)(a(1) - a(3) + 2) ... (a(1) - a(n) + n - 1)   *
#                      (a(2) - a(3) + 1) ... (a(2) - a(n) + n - 2)   *
#                                   ...                              *
#                                            (a(n - 1) - a(n) + 1)   *
#       multinomial(a(1) + a(2) + ... + a(n), a(1), a(2), ..., a(n)) /
#       ((a(1) + n - 1) * (a(1) + n - 2) * ... * (a(1) + 1))         /
#       ((a(2) + n - 2) * (a(2) + n - 3) * ... * (a(2) + 1))         /
#                                   ...
#                                             (a(n - 1) + 1)