#!/usr/local/bin/maple

# Nathaniel Shar
# HW 2
# Experimental Mathematics

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

Help := proc(): print("GC(x,L,n), GLP(x,L,n), GLP3(x, L, n), CW(r), CWgen(r, n), Parent(r), kCousins(r, k)"): end:

# One iterate of the generalized Collatz problem

GC := proc(x,L,n) local m, r:
    m := nops(L):
    r := n mod m:
    if r = 0 then
        r := m:
    fi:
    subs(x=n, L[r]):
end:

# Returns a 2-element list. The second element is 1 if a cycle resulted
# within maxdepth iterations, and 0 if it didn't. Also, if a cycle was
# found, the the number of iterations required to reach that cycle is
# put in the first element. Otherwise, maxdepth is put in the first
# element.

GLP := proc(x,L,n, maxdepth) local i,y,currentval, seenvals:

    y:=n:
    seenvals := []:
    for i from 0 while (not(y in seenvals) and i < maxdepth) do
        seenvals := [op(seenvals), y]:
        y:=GC(x,L,y):
    od:
    if i = maxdepth then
        return [i, 0];
    else
        return [i, 1];
    fi:
end:

# GLP3 returns a 3-element list, where the first two elements are as
# with GLP, while the last element is the value obtained on the last
# iteration. (This is so I can check whether two cycles are the same.)

GLP3 := proc(x,L,n, maxdepth) local i,y,currentval, seenvals:

    y:=n:
    seenvals := []:
    for i from 0 while (not(y in seenvals) and i < maxdepth) do
        seenvals := [op(seenvals), y]:
        y:=GC(x,L,y):
    od:
    if i = maxdepth then
        return [i, 0, y];
    else
        return [i, 1, y];
    fi:
end:

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

# Simple experiment 1: Collatz functions of the form [a*x+1, x/2], a odd.
# Observations: If a = 1 or a = 3, a loop is always reached. In the
# case that a = 1, we can even prove this fact!

# If a = 5, a loop does not appear to be reached when a = 7. I can't
# prove this. However, no loop was found after 40000 iterations. The
# 40000th number in the sequence had 1239 digits, so it seems likely
# that the iterates are growing without bound.

# If a = 7, only a few initial values fall into cycles quickly (most obviously
# the powers of 2). For a >= 9 and not 2^k-1 for some k, it appears
# that there are likely not to be any cycles at all.

# Simple experiment 2: Try functions of the form [a*x+b, x/2], a and b
# odd, b <> 1.

# Once again, we get a similar result. For small values of a (that is,
# less than or equal to 9), some numbers fall into loops. For large
# values of a, it appears that no numbers do, unless an + b = n2^k for
# some positive integer n. In that case n falls into a loop. It is
# easy to see that this trivial form of loop (with only one ascent)
# does not occur for some combinations of values. For example, if b is
# prime, then an+b = n2^k means n = b/(2^k-a). Therefore 2^k-a is in
# {+1, -1, +b, -b}. For example, if a = 19 and b = 11, 2^k is in the
# set {20, 18, 30, -8}, which is a contradiction, so no such loop
# exists.

# However, it's hard to say a lot more than that.
# For example, if a = 7 and b = 5, then there is a loop of length 43
# that includes 7. This breaks a lot of conjectures.


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

# The sequence of iterates appears to stabilize at 0 for all n <= 10^6.
# It is very plausible that it always returns to 0, since
# it only increases at the multiples of 3, and in those cases it only
# increases by a factor of 4/3 while the rest of the time it decreases
# by 3.

# However, I can't prove this. For any k, there are clearly values at which the
# iterates increase k times in a row (for example, 3^k). This wrecks a
# lot of naive attempts at a proof.

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

# CW: from class2a

CW := proc(r) local a,b:
    a := numer(r):
    b := denom(r):
    return [a/(a+b),(a+b)/b]:
end:

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

# CWgen: Returns the nth generation of the CW tree rooted at r.
# I don't actually need this function but I wrote it as practice.

CWgen := proc(r, n) local i:
    if n <= 0 then return [r]: fi:
    for i from 1 to n do:
        return [seq(op(CW(s)), s = CWgen(r, n-1))]:
    od:
end:

# Return the parent of the rational number r in the Calkin-Wilf
# construction. 1 has no parent.

Parent := proc(r) local a, b:
    if r = 1 then return FAIL: fi:
    a := numer(r):
    b := denom(r):
    if (b >= a) then
        return a/(b-a):
    else:
        return (a-b)/b:
    fi:
end:

# kCousins(r, k): returns the k-cousins of r in the Calkin-Wilf
# construction. 0-cousins are siblings.

kCousins := proc(r, k):
    if r = FAIL then
        return []:
    elif r = 1 then:
        return []:
    elif k = 0 then
        return convert(CW(Parent(r)), set) minus {r}:
    else:
        return [seq(op(CW(s)), s = kCousins(Parent(r), k-1))]:
    fi:
end:

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

#### NOTICE: I have seen this sequence before. Although I didn't know the
#### generating function, I did know enough to come up with it significantly
#### more easily than I would have otherwise. As a result I do not believe
#### I am eligible for a share of the $5.

##################################


# Observe that the sequence is produced by the recurrence w(2n) =
# w(n), w(2n+1) = w(n) + w(n+1), and the initial value w(1) = 1
# This is immediate from the Calkin-Wilf construction.

# However, it will be easier to handle the generating function of the
# sequence v(n) = w(n+1). In this case the recurrence is v(2n+2) = v(n)
# + v(n+1), v(2n+1) = v(n), and the initial condition is v(0) = 1.

# Translating this into the language of generating functions, if
# V(x) is the generating function of v(n), then
#           V(x) = (1 + x^2)V(x^2) + xV(x^2).
# In other words V(x) = (1+x+x^2)V(x^2).

# Now, we can simply put
#        V(x) = (1 + x + x^2)(1 + x^2 + x^4)(1 + x^4 + x^8)...,
# It is evident that each factor is equal to the previous factor
# evaluated at x^2, so that V(x^2) is equal to the product of all but
# the first factor. We can also see that [x^0]V(x) = 1, so V(x) is the
# desired generating function.

# Of course W(x) = xV(x), so
#        W(x) = x(1+x+x^2)(1+x^2+x^4)(1+x^4+x^8)(1+x^8+x^16)...

# This should be a sufficiently nice form for the generating function.

# As a check, we may note that the first few terms are
# x + x^2 + 2x^3 + x^4 + 3x^5 + 2x^6 + 3x^7 + x^8 + ...
# which is correct.