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

# 1.

UDbetter := proc(n)
    if n = 0 then return
        [[]]:
    fi:

    if n = 1 then return
        [[1]]:
    fi:

    res := []:

    # Insert n into position 2 * k.
    for k from 1 to floor(n / 2) do
        all_left := UDbetter(2 * k - 1):
        all_right := UDbetter(n - 2 * k):
        for L in combinat[choose](n - 1, 2 * k - 1) do
            R := select(k -> not k in L, [seq(1..n)]):
            for left_ud in all_left do
            for right_ud in all_right do
                new_left := Expa(left_ud, L):
                new_right := Expa(right_ud, R):
                new := [op(new_left), n, op(new_right)]:
                res := [op(res), new]:
            od:
            od:
        od:
    od:

    res:
end:

#Expa(pi,S): inputs a permutation pi of 1,..., nops(pi) and
#a sorted list of distinct numbers S, replaces 1 by
#the smallest member of S, etc. 
#For example
#Expa([2,1,3],[11,13,15])= [13,11,15]
Expa:=proc(pi,S) local i:
subs({seq(i=S[i],i=1..nops(S))},pi):
end:

# 2.

# Let U(z) be the egf of the up-down numbers. We claim that

#       U(z) = sec(z) + tan(z),

# or equivalently that

#       cos(z) * U(-z) = 1 - sin(z).

# cos(z) is the (-1)^n-weighted egf of increasing sequences of even length, and
# U(z) is the (-1)^n-weighted egf of up-down sequences. Their product is the
# weighted egf of pairs of sequences of the form

#       a1 < a2 < ... < ak; b1 < b2 > ... (<>) b(n - k),

# where [n] = {a1, a2, ..., ak, b1, b2, ..., b(n - k)}, k is even, and the
# weight is (-1)^n.

# The idea behind the killer involution here is to move an element from both
# sides of the a's to the b's to extend the up-down sequence. (But I can't
# actually remember how to do this, and I waited too long to think about it :))

# 3.

piAppx := proc(n)
    total := 0:
    front := 1:
    power := 1:
    for k from 1 to n do
        front := -1 * front * (2 * k - 1) / (2 * k):
        power := power * 2^10:
        top := 20 * k * (50 + 41 * (k - 1)) + 13:
        total := total + (front)^5 * top / power:
    od:

    total := 13 + total:
    sqrt(128 / total):
end:

Digits := 1000:

# Happy Pi day!
piAppx(1000);
evalf(piAppx(100));