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

read "C6.txt":

# 1 + 5.

# I'll just write the multi-dimensional version first.

Nonnegative := l -> andmap(x -> x >= 0, l):

# I don't know what "subdiagonal" means in N > 2 dimensions, so I made up the
# following definition. Essentially, "subdiagonal" means you're "underneath"
# the hyperplane with normal vector (-1, 1, -1, ...).
subdiagonal := l -> add((-1)^k * l[k], k=1..nops(l)) <= 0:

GNuWalks := proc(A, pt)
    option remember:
    if nops({seq(nops(a), a in A)}) > 1 then
        error("atomic steps should be the same dimension");
        return FAIL
    fi:

    d := nops(A[1]):

    GnuWalksRec := proc(pt)
        if pt = [0 $ d] then
            return 1:
        end:

        total := 0:

        for step in A do
            if Nonnegative(pt - step) and subdiagonal(pt - step) then
                total := total + GNuWalks(A, pt - step):
            fi:
        od:

        total:
    end:

    GnuWalksRec(pt):
end:

NuWalks := proc(A, pt)
    option remember:
    if nops({seq(nops(a), a in A)}) > 1 then
        error("atomic steps should be the same dimension");
        return FAIL
    fi:

    d := nops(A[1]):

    NuWalksRec := proc(pt)
        if pt = [0 $ d] then
            return 1:
        end:

        total := 0:

        for step in A do
            if Nonnegative(pt - step) then
                total := total + NuWalks(A, pt - step):
            fi:
        od:

        total:
    end:

    NuWalksRec(pt):
end:

# 2.

# A984 (obvious!)
seq(NuWalks([[1,0],[0,1]],[i,i]),i=1..20);

# A108 (famous!)
seq(GNuWalks([[1,0],[0,1]],[i,i]),i=1..20);

# A36692 (an Apagodu-Zeilberger sequence!)
seq(NuWalks([[1,0],[0,1],[2,0],[0,2]],[i,i]),i=1..20);

# A122951 (literal description)
seq(GNuWalks([[1,0],[0,1],[2,0],[0,2]],[i,i]),i=1..20);

# A122680 (another Apagodu-Zeilberger sequence!)
seq(NuWalks([[1,0],[0,1],[2,0],[0,2],[3,0],[0,3]],[i,i]),i=1..20);

# A175883 (literal description)
seq(GNuWalks([[1,0],[0,1],[2,0],[0,2],[3,0],[0,3]],[i,i]),i=1..20);

# A1850 (well-known)
seq(NuWalks([[1,0],[0,1],[1,1]],[i,i]),i=1..20);

# A6318 (not obvious!)
seq(GNuWalks([[1,0],[0,1],[1,1]],[i,i]),i=1..20);

# Not in there!
seq(NuWalks(FootBall(),[i,i]),i=1..20);

# Not in there!
seq(GNuWalks(FootBall(),[i,i]),i=1..20);

# I was curious if a(n) = "number of A-walks from (0, 0) to (n, n)" is always
# C-finite. This is not true (see A1850), but here's the function I was playing
# with.

randomSteps := proc(d, N, K)
    r := rand(0..K):
    res := {}:

    while nops(res) < N do
        step := [seq(r(), j=1..d)]:
        if Nonnegative(step) and add(step) > 0
            then res := res union {step}:
        fi:
    od:

    res:
end:

# 3.

# (i)

S := {[0, 1], [1, 0]}:

seq(GNuWalks(S, [n, n]) / NuWalks(S, [n, n]), n=0..20);
seq(evalf(GNuWalks(S, [n, n]) / NuWalks(S, [n, n])), n=0..20);

# I conjecture that the probability is 1 / (n + 1) exactly.

# (ii)

# The subdiagonal (n, n) walks are enumerated by the Catalan numbers.  The (n,
# n) walks are enumerated by the central binomial coefficients.  Thus, the
# probability is

#       binomial(2n, n) / (n + 1) / binomial(2n, n) = 1 / (n + 1),

# as we predicted.

# 4.

footballPs := 0:

for n from 100 to 200 do
    d := NuWalks(FootBall(), [n, n]):
    if NuWalks(FootBall(), [n, n]) > 0 then
        p := GNuWalks(FootBall(), [n, n]) / d:
        footballPs := footballPs + p * x^n:
    fi:
od:

# If the probability is n^alpha for some alpha, then we can discover alpha
# experimentally by computing log(probability(n)) / log(n) for large n. Doing
# this here gives alpha ~ -0.891575.

evalf(log(lcoeff(footballPs, x)) / log(degree(footballPs, x)));

# 5. Procedure done above! (It may not be what everyone had in mind, though.)

# A6480 (literal description)
seq(NuWalks([[1,0,0],[0,1,0],[0,0,1]],[i,i,i]),i=1..10);

# None!
seq(GNuWalks([[1,0,0],[0,1,0],[0,0,1]],[i,i,i]),i=1..10);

# A8977 (obvious)
seq(NuWalks([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],[i,i,i,i]),i=1..10);

# None!
seq(GNuWalks([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],[i,i,i,i]),i=1..10);