#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 5
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

Help := proc(): print(`BallWithPaths(G,i,k), SphereWithPaths(G,i,k), Knk(n,k), Mnk(n,k,r), MMChildren(n, pos), MMpos2int(n,pos), MMint2pos(n,i), MMDigraph(n), MMSolve(n), PegSoliChildren(M), PegSolipos2int(n,k,pos), PegSoliint2pos(n,k,r), isPowerOfTwo(n), PegSoliDigraph(n, k), PegSoliSolve(M)`): end:

#############
# Problem 2 #
#############

BallWithPaths := proc(G,i,k) local path, v, seen, B, S: option remember:
    if k = 0 then:
        return {[i]}:
    else:
        B := BallWithPaths(G,i,k-1):
        S := B:
        seen := {}:
        for path in B do:
            for v in G[path[-1]] do:
                if  not v in seen then:
                    S := S union {[op(path), v]}:
                    seen := seen union {v}:
                fi:
            od:
        od:
        return S:
    fi:
end:

SphereWithPaths := proc(G,i,k) local B, L:
    if k = 0 then:
        return BallWithPaths(G,i,k):
    else:
        B := BallWithPaths(G, i, k-1):
        L := map(x->x[-1], B):
        return select(x->(not x[-1] in L), BallWithPaths(G,i,k)):
    fi:
end:

#############
# Problem 3 #
#############

Knk := proc(n,k) local S, ill, Sr, i,j,l,M: option remember:
    M := {seq(i,i=1..n)}:
    if k = 0 then:
        return {[]}:
    fi:
    S := Knk(n,k-1):
    Sr := {}:
    for i in S do:
        ill := {seq(i[j], j=1..k-1)}:
        if k >= 2 then:
            ill := ill union {i[k-1]-1, i[k-1]+1}:
        fi:
        Sr := Sr union {seq([op(i), l], l in M minus ill)}:
    od:
    return Sr:
end:

# The sequence is A2464.

#############
# Problem 4 #
#############

Mnk := proc(n,k,r) local S, ill, Sr, i,j,l,M: option remember:
    M := {seq(i,i=1..n)}:
    if k = 0 then:
        return {[]}:
    fi:
    S := Mnk(n,k-1,r):
    Sr := {}:
    for i in S do:
        ill := {seq(i[j], j=1..k-1)}:
        for j from 1 to min(k-1, r) do:
            ill := ill union {i[k-j]-j, i[k-j]+j}:
        od:
        Sr := Sr union {seq([op(i), l], l in M minus ill)}:
    od:
    return Sr:
end:

# The beginning of the sequence Mnk(i,i,r) is in OEIS for all integers
# r, because all the sequences share r initial terms with
# A170 :)

# However, the sequences only appear for r = 0, 1. When r=0, the
# sequence is A142.

#############
# Problem 5 #
#############

# pos = [number of missionaries on left bank, number of cannibals on
# left bank, position of boat]

MMChildren := proc(n, pos) local ml, cl, mr, cr, boat, S:
    S := {}:
    (ml, cl, boat) := op(pos):
    mr := n-ml:
    cr := n-cl:
    if boat = 0 then:
        if (ml >= cl + 2 or ml=2) and mr + 2 >= cr then:
            S := S union {[ml-2, cl, 1]}:
        fi:
        if ml >= cl and mr >= cr and cl >= 1 then:
            S := S union {[ml-1, cl-1, 1]}:
        fi:
        if cl >= 2 and (mr >= cr+2 or mr = 0) then:
            S := S union {[ml, cl-2, 1]}:
        fi:
        if (ml >= cl+1 or ml=1) and mr+1 >= cr then:
            S := S union {[ml-1, cl, 1]}:
        fi:
        if cl >= 1 and (mr >= cr+1 or mr = 0) then:
            S := S union {[ml, cl-1, 1]}:
        fi:
    fi:
    if boat = 1 then:
        if (mr >= cr + 2 or mr=2) and ml + 2 >= cl then:
            S := S union {[ml+2, cl, 0]}:
        fi:
        if mr >= cr and ml >= cl and cr >= 1 then:
            S := S union {[ml+1, cl+1,0]}:
        fi:
        if cr >= 2 and (ml >= cl+2 or ml = 0) then:
            S := S union {[ml, cl+2, 0]}:
        fi:
        if (mr >= cr+1 or mr=1) and ml+1 >= cl then:
            S := S union {[ml+1, cl, 0]}:
        fi:
        if cr >= 1 and (ml >= cl+1 or ml = 0) then:
            S := S union {[ml, cl+1, 0]}:
        fi:
    fi:
    return S:
end:

MMpos2int := proc(n,pos):
    return 2*(pos[1]*(n+1)+pos[2])+pos[3]+1:
end:

MMint2pos := proc(n,i) local j, k, l:
    l := (i-1) mod 2:
    k := ((i-l-1)/2) mod (n+1):
    j := (((i-l-1)/2-k))/(n+1):
    return [j,k,l]:
end:

MMDigraph := proc(n):
    return [seq(map(x->MMpos2int(n, x), MMChildren(n,MMint2pos(n,j))), j=1..2*(n+1)^2)]:
end:

MMSolve := proc(n) local first, last, i, p, G, S:
    first := MMpos2int(n, [n,n,0]):
    last := MMpos2int(m, [0,0,1]):
    G := MMDigraph(n):
    for i from 0 do:
        S := SphereWithPaths(G, first, i):
        if S = {} then:
            return FAIL:
        else:
            if last in {seq(p[-1], p in S)} then:
                return map(x->MMint2pos(n, x), op(select(x->x[-1]=last, S))):
            fi:
        fi:
    od:
end:

# Unfortunately, after all this effort it appears the problem is
# unsolvable for n >= 3.

#############
# Problem 6 #
#############

PegSoliChildren := proc(M) local n, k, i, j, S, Mnew:
    k := nops(M):
    if k = 0 then:
        return FAIL:
    fi:
    n := nops(M[1]):
    S := {}:
    for i from 1 to k do:
        for j from 1 to n do:
            if M[i][j] = 1 then:
                if j >= 3 and M[i][j-1] = 1 and M[i][j-2] = 0 then:
                    Mnew := M:
                    Mnew[i][j] := 0:
                    Mnew[i][j-1] := 0:
                    Mnew[i][j-2] := 1:
                    S := S union {Mnew}:
                fi:
                if j <= k-2 and M[i][j+1] = 1 and M[i][j+2] = 0 then:
                    Mnew := M:
                    Mnew[i][j] := 0:
                    Mnew[i][j+1] := 0:
                    Mnew[i][j+2] := 1:
                    S := S union {Mnew}:
                fi:
                if i >= 3 and M[i-1][j] = 1 and M[i-2][j] = 0 then:
                    Mnew := M:
                    Mnew[i][j] := 0:
                    Mnew[i-1][j] := 0:
                    Mnew[i-2][j] := 1:
                    S := S union {Mnew}:
                fi:
                if i <= n-2 and M[i+1][j] = 1 and M[i+2][j] = 0 then:
                    Mnew := M:
                    Mnew[i][j] := 0:
                    Mnew[i+1][j] := 0:
                    Mnew[i+2][j] := 1:
                    S := S union {Mnew}:
                fi:
            fi:
        od:
    od:
    return S:
end:

#############
# Problem 7 #
#############

PegSolipos2int := proc(n,k,pos) local exp, val, i, j:
    val := 0:
    exp := 0:
    for i from 1 to k do:
        for j from 1 to n do:
            if pos[i][j] = 1 then:
                val := val + 2^exp:
            fi:
            exp := exp+ 1:
        od:
    od:
    return val+1:
end:

PegSoliint2pos := proc(n,k,r) local v, M, a,b, i, j:
    v := r-1:
    M := [seq([seq(0, a=1..n)], b=1..k)]:
    for i from 1 to k do:
        for j from 1 to n do:
            if v mod 2 = 1 then:
                M[i][j] := 1:
            fi:
            v := (v-(v mod 2))/2:
        od:
    od:
    return M
end:

isPowerOfTwo := proc(n):
    if type(log[2](n), integer) then:
        return true:
    else:
        return false:
    fi:
end:

PegSoliDigraph := proc(n, k):
    return [seq(map(x->PegSolipos2int(n,k, x), PegSoliChildren(PegSoliint2pos(n,k,j))), j=1..2^(n*k))]:
end:

PegSoliSolve := proc(M) local n,k, first, i, p, G, S:
    k := nops(M):
    if k = 0 then:
        return FAIL:
    fi:
    n := nops(M[1]):
    first := PegSolipos2int(n, k, M):
    G := PegSoliDigraph(n,k):
    for i from 0 do:
        S := SphereWithPaths(G, first, i):
        if S = {} then:
            return FAIL:
        else:
            if select(x->isPowerOfTwo(x-1), {seq(p[-1], p in S)}) <> {} then:
                return map(x->PegSoliint2pos(n, k,x), op(select(x->isPowerOfTwo(x[-1]-1), S))):
            fi:
        fi:
    od:
end:

# The grid with 1 peg missing from the corner is not solvable.
# This methodology cannot be reasonably extended to the English
# board. If we want to work on the English board, we had better not
# pregenerate the digraph.