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

# Due to being harassed by the government and having exams to grade, my
# homework this week is not so good.

# 1.

# Given a hypergeometric term F and polynomial P in n such that P(n) F(n) -> 0,
# compute the series of (-1)^n P(n) F(n) within d digits of accuracy.
ComputeCons := proc(F, P, n, d)
    local R, total, term, k:
    # P * F is a hypergeometric term, and therefore
    #   P(n + 1) * F(n + 1) = R(n) * P(n) F(n)
    # for some fixed rational function R.
    # (Note the conversion to factorials, so that maple will actually simplify
    # the quotient correctly.)
    R := simplify(convert(subs(n = n + 1, F * P) / F / P, factorial)):

    total := 0:
    term := eval(F * P, n = 0):

    for k from 0 while term >= 10^(-d) do
        total := total + (-1)^k * term:
        term := eval(R, n = k) * term:
    od:

    total:
end:

naiveAdd := (k, N) -> add((-1)^n * (n^3 + 1) * k^n / binomial(2 * n, n), n=0..N):
myAdd := (k, d) -> ComputeCons((n^3 + 1), k^n / binomial(2 * n, n), n, d):

evalf(myAdd(2, 500) - naiveAdd(2, 1800));
time(naiveAdd(2, 1800));
time(myAdd(2, 500));

# myAdd is better when the terms don't decay *too* quickly. Otherwise it isn't
# so clear.

# (Also, I think that series might have a closed form...)

# 2.

# We want to construct the reduced n by k latin trapezoids.
# Obviously, there's only one n by 1 trapezoid: [1, 2, ..., n].
# For k > 1, we need to add another row onto an n by (k - 1) latin trapezoid.
# So just recursively construct those and check all the possible rows we could
# add. That's what LT does.

# It's going to have not very good runtime.

read "C16.txt":

# Should be all 1's.
seq(nops(LT(k, 1)), k=1..10);

# 4, 236 never appears in the OEIS, so it's unlikely that this sequence is in
# there.
seq(nops(LT(k, k)), k=1..6);

# It seems like this sequence *is* in the OEIS!
# Probably A127548. The explicit description is "ogf = sum(n! * (x / (1 +
# x)^2)^n, n >= 0)."
seq(nops(LT(k, 2)), k=1..7);

# I'm told that George has a proof of this fact - huzzah! He doesn't think
# that it generalizes, but I'm more interested in how the *generating function*
# generalizes. Surely we could guess the form of the ogf for LT(n, 3)?

# Also, I discovered an "error" in the OEIS. A013999 is described as "the
# binomial transform" of A000271, but this is not true. If this were true, then
# our sequence would equal A000271. Instead, A013999 is the binomial transform
# of A000271 after dropping the initial 1.