## Homework 2, done by Patrick Devlin. 25-January-2012 ##
### I don't have any preferences about keeping my homework private. ###
Help:=proc(): print(` GC(x,L,n), GCFindCycle(x,L,n,M), FindConjectures(modWhat,maxCoeffs,howFar), FindConjecturesMod2(maxCoeffs,howFar), testACase(x,L,n), checkProblem4(n), Parent(x), kCousins(x,k), kthGeneration(x,k) `); end proc:
### Problem 1 ###
# Not done since it is for novices
### Problem 2 ###
GC:=proc(x,L,n) local i:
i:= (n mod nops(L)):
if(i = 0) then i:=nops(L): fi:
return eval(L[i], x=n):
end proc:
### Problem 3 ###
## As is fitting for this course, the process of formulating analogs to the 3n+1 conjecture has been automated below
## Type FindConjectures(modWhat, maxCoeffs, howFar) to systematically run through stuff, and the program spits out conjectures
## For instance, it would say stuff like "L = [...] seems to work for all x < ..., and all its cycles are [...].". and
## "L=[...] doesn't work. It looks like x = ... is the smallest counterexample."
## In running through all possible cases, it uses maybeInteresting(x,L) to skip the functions where it is "obvious" what will happen
## (e.g., it skips functions that are [at least eventually] strictly increasing/decreasing, as well as
## functions such as [x+2, x/2], where it is `obvious` that it will increase without bound for every odd)
## You can also use the function testACase(x,L,n) to have it spit out a conjecture for a particular case.
## This program has so far generated hundreds of thousands of (seemingly interesting!) Collatz-type
## conjectures (both for "it looks like THIS ONE works" and "it looks like THIS ONE doesn't").
## For one of my own more human conjectures, I would say [ax + 1, x/2] works if and only if a=1 or a=3.
## This particular conjecture can be tested with FindConjecturesMod2(maxCoeffs, howFar)
## Another conjecture (that seems false now) is as follows:
# There exists a constant C(n) [depending only on n] such that
# any Collatz function F:=[f1, f2, f3, ... , fn] works if and only if
# it works for all inputs up to C(n)
##
## If it were true, then I now ``know" [believe experimentally] that C(2) > 197, and C(3) > 7.
##
## The conjecture can be weakened to something like...
# For all n, there exists a function p(u1, v1, u2, v2, ... , un, vn)
# such that any Collatz function F(x) := [a1 * (x-1)/n + b1, a2 * (x-2)/n + b2, ... , an * x/n + bn]
# works if and only if F(x) works for all x up to p(a1, b1, a2, b2, ... , an, bn).
##
## This conjecture has the benefit of being "trivially true" (since the
## function p could be anything under the sun), but the conjecture at least implies that
## the function p is somehow nice.
#
## The conjecture has the downside of being hopelessly vague and thus untestable.
#
## Hopefully the family of functions corresponding to various values of n
## all sort of look the same in some nice sense.
#*
## Future steps would be to have the program determine all possible fixed points and do some analysis there.
## For example, if we knew that any cycle of [7x + 3, x/2] needed to contain 1, then we could prove
## that the function has no cycles whatsoever (i.e., fails at every input)
## since 7x+3 = 2^n has no integer solutions [because 2^n = 3 (mod 7) has no solutions]
#*
GCFindCycle:=proc(x::uneval,L,n,M,howOftenToCheck) local i, i2, C, n2, L2: # This iterates GC(x,L,n) until either the function cycles (in which case it returns the smallest element in the cycle and the number of iterations needed to reach that element) or it's called more than M times (in which case it gives up)
n2:=n:
C:=[]:
i2:=0:
for i from 1 to M do
C:=[op(C), n2]:
n2:=GC(x,L,n2):
i2:=i2+1:
if(i2 > howOftenToCheck) then
L2:=GCIsCycle(C,n2):
i2:=0:
if (nops(L2) >1) then return L2: fi:
fi:
od:
return GCIsCycle(C,n2):
end proc:
GCIsCycle:=proc(C,n) local i, i2, minEl, minInd: # tests if there is a cycle in the list C or not
for i from 1 to nops(C) do # check to see if there's a cycle
if(n=C[i]) then # then there's a cycle!
minEl:=C[i]:
minInd:=i:
for i2 from i to nops(C) do # find smallest element in the cycle
if(C[i2] < minEl) then minEl:=C[i2]: minInd:=i2: fi:
od:
return [minEl, minInd]:
fi:
od:
return []:
end proc:
FindConjectures:=proc(modWhat, maxCoeffs, howFar) local x, L, coeffs, consts, i, i2, i3, tempCoef, tempCon: # This systematically suggests and tests various generalized Collatz functions
for i from 0 to ((maxCoeffs-1)^(modWhat)-1) do
for i2 from 0 to ((maxCoeffs)^(modWhat)-1) do
tempCoef:=i:
tempCon:=i2:
coeffs:=[]: consts:=[]:
for i3 from 1 to modWhat do
coeffs:=[op(coeffs), ( (tempCoef mod (maxCoeffs-1)) + 1)]:
consts:=[op(consts), (tempCon mod maxCoeffs)]:
tempCoef:= floor(tempCoef / (maxCoeffs-1)):
tempCon:= floor(tempCon / maxCoeffs):
od:
L:=[]:
for i3 from 1 to (modWhat-1) do:
L:=[op(L), (coeffs[i3]*(x-i3)/modWhat+consts[i3])]:
od:
L:=[op(L), coeffs[modWhat]*x/modWhat+consts[modWhat]]:
if(maybeInteresting(x,L)) then
testACase(x,L,howFar):
fi:
od:
od:
end proc:
maybeInteresting:=proc(x,L) local allIncreasing, allDecreasing, i, t, k: # does a few quick sanity checks to see if the function determined by L is ``obviously" stupid
allIncreasing:=true:
allDecreasing:=true:
for i from 1 to nops(L) do
t:=eval(diff(L[i],x),x=0): # let t be slope at zero
if t>1 then allDecreasing:=false: fi: # one coordinate is increasing
if t<1 then allIncreasing:=false: fi: # one coordinate is decreasing
if t=1 then # some careful thought needed
t:=eval(L[i],x=0): # now t is y-intercept
k:=(t+i mod nops(L)):
if(k=0) then k:=nops(L): fi:
if(k=i) then
if(t<>0) then return false: else allIncreasing:=false: fi: # it is not interesting because this coordinate will increase (or decrease!) indefinitely
else
t:=eval(diff(L[k],x),x=0): # checking to see if it's [eventually] increasing/decreasing after one iteration
if t<1 then allIncreasing:=false: fi: # it's eventually decreasing
if t>1 then allDecreasing:=false: fi: # it's eventually increasing
if t=1 then # if t=1, then it's x+c, and then the next iteration is x+d, calling for some analysis
t:=eval(L[k],x=0)+eval(L[i],x=0): # then starting at x, F(x) = x+c, and F(F(x)) = x+c+d = x+t
k:=(t+i mod nops(L)): if(k=0) then k:=nops(L): fi: # k is what gets called NEXT
if((k = i) and (t<>0)) then return false: fi: # it gets stuck in an infinite loop between these two states, and it's unbounded
if((t>0) and (eval(L[i],x=0)>0)) then allDecreasing:=false: fi: # then F(F(x)) > x and F(x) > x
if((t<=0) and (eval(L[i],x=0)<=0)) then allIncreasing:=false: fi: # then F(F(x)) <= x and F(x) <= x. (so it's not strictly increasing)
fi:
fi:
fi:
od:
if ( (allIncreasing) or (allDecreasing)) then return false: else return true: fi:
end proc:
FindConjecturesMod2:=proc(maxCoeffs, howFar) local x, L, coeffs, i, i2, i3, tempCoef, tempCon: # This systematically suggests and tests various generalized Collatz functions all mod 2
for i from 0 to maxCoeffs do
L:=[(2*i+1)*x+1, x/2];
testACase(x,L,howFar):
od:
end proc:
testACase:=proc(x,L,n) local depth, cycles, i, i2, newCycle, tempRes: # this tests one case and sees what cycles seem to happen. If it looks like the iterated function settles into a few cases, it outputs that as a conjecture. Else, it conjectures that the function determined by L can go to infinity.
cycles:=[]:
depth:=1000:
tempRes:=[]:
for i from 1 to n do
tempRes:=GCFindCycle(x,L,i,depth, 10):
if (nops(tempRes) < 2) then
tempRes:=GCFindCycle(x,L,i,10000, 10000):
fi:
if (nops(tempRes) < 2) then
print(`Fail at `, i, `for L=`, L):
return:
fi:
newCycle:=true:
for i2 from 1 to nops(cycles) do
if(cycles[i2]=tempRes[1]) then newCycle:=false: fi:
od:
if(newCycle) then cycles:=[op(cycles), tempRes[1]]: fi:
od:
if(nops(cycles) < 10) then
print(`Pass! It seems like L=`, L, ` works for n<= `, n, ` and all go into one of the cycles `, cycles):
else
print(`Pass! It seems like L=`, L, ` works for n<=`, n, ` but there are `, nops(cycles), `cycles so far`):
fi:
return:
end proc:
### Problem 4 ###
# For L:=[(x-1)/3, (x-2)/3, 4 x/3], it appears that every n is sent to 0.
# In fact, suppose n is the smallest positive number that is not sent to 0.
# Then after some brute-force checking, we have that n >= 1000.
# Moreover, the smallest counter example would need to be divisible by 3^5, so
# we need only check multiples of 3^5. For this purpose, the following
# procedure was made:
## (Running this for a little bit shows the smallest counter example is
## greater than 243 million)
## This could probably be pushed further by restructuring the below program, but perhaps not by a whole lot
checkProblem4:=proc(n) local depth, tempRes, i, i2,L,x: #this is a very slightly modified version of testACase
L:=[(x-1)/3, (x-2)/2, 4*x/3]:
depth:=1000:
tempRes:=[]:
for i from 1 to n do
tempRes:=GCFindCycle(x,L,i*(3^5),depth, 10):
if (nops(tempRes) < 2) then
tempRes:=GCFindCycle(x,L,i*(3^5),100000, 100000):
fi:
if (nops(tempRes) < 2) then
print(`Fail at `, i*(3^5), `for L=`, L):
return:
fi:
od:
print(`Pass! It seems like L=`, L, ` works for n<= `, n*(3^5)):
return:
end proc:
### Problem 5 ###
Parent:=proc(x):
if((x=1) or x<=0) then return -1; fi; # returns -1 as a flag for "x has no parent"
return max(x/(1-x), x-1);
end proc: # since x > -1 and x neq 1, both of these expressions are rational numbers, and exactly one of them is positive. Moreover, by simplifying the terms a/(a+b) sort of representations, 1/(1-x) and x-1 are the inverse functions for LeftChild(x) and RightChild(x) respectively
### Problem 6 ###
kCousins:=proc(x,k) local i, y, errorMessages: # inputs rational number x and finds its k^th cousins in order from left to right
errorMessages:=false:
if(x <= 0) then if(errorMessages) then print(`Cannot work for nonpositive inputs of x`): fi: return []: fi:
if(k<0) then if(errorMessages) then print(`Cannot work for negative inputs of k`): fi: return []: fi:
y:=x;
for i from 1 to k do
y:=Parent(y):
od:
if((y<=0) or (y=1)) then if(errorMessages) then print(`Input does not have enough ancestors`): fi: return []: fi:
return kthGeneration(max((1-1/y), (1/(1-y))),k): # since y>0 and y neq 1, exactly one of the elements in max is positive, which coincides with y's brother
end proc:
kthGeneration:=proc(x,k): #returns kth generation of the number x (note this generalizes to non-rational inputs for x)
if(x<=0) then return []: fi:
if(k<=0) then return [x]: fi:
return [op(kthGeneration(x/(x+1), k-1)), op(kthGeneration(x+1,k-1))]:
end proc:
### Problem 7 ###
## Generating function for w[n] as a (convergent) infinite product ##
# Let a[n] = w[n-1] for n=1, 2, ...
# Then from the tree structure of the construction, for all n>=1 we have
# a[n] = a[2n], and a[n] + a[n+1] = a[2n+1].
# Thus, it follows that for all n>=0, we have w[n] = w[2n+1], and w[n] + w[n+1] = w[2n+2].
# Therefore, if U(x) = w[0] + w[1] x + w[2] x^2 + ... , then we have
# x U(x^2) = w[0] x + w[1] x^3 + w[2] x^5 + w[3] x^7 + ... = w[1] x + w[3] x^3 + w[5] x^5 + w[7] x^7 + ....
#
# Moreover, we also have
# (1+ x^2) U(x^2) = w[0] + (w[0]+w[1]) x^2 + (w[1]+w[2]) x^4 + (w[2]+w[3]) x^6 + ... = w[0] + w[2] x^2 + w[4] x^4 + w[6] x^6 + ...
# Therefore, it follows that:
# U(x) = x U(x^2) + (1+x^2) U(x^2) = U(x^2) (1 + x + x^2).
# Then by substituting y=x^2, and using U(y) = (1 + y + y^2) U(y^2) , we obtain:
# U(x) = (1 + x + x^2) (1 + x^2 + x^4) U(x^4), and in general, we have
# U(x) = (1 + x + x^2) (1 + x^2 + x^4) (1 + x^4 + x^8) (1 + x^8 + x^16) (1 + x^16 + x^32) ... (1 + x^(2^k) + x^(2^(k+1))) U(x^(2^(k+1))).
# Then since U(x^m) = 1 + w[1] x^m + ... , we can "take limits as k->infinity" in the above expression
# to obtain that [as formal algebraic power series]:
#
# U(x) = product( (1 + x^(2^k) + x^(2^(k+1))), k=1..infinity).
# = (1 + x + x^2) (1 + x^2 + x^4) (1 + x^4 + x^8) (1 + x^8 + x^16) ...
# If you care (which you probably don't) we know that the right hand side of the above
# equation converges to a nonzero holomorphic function in the unit disc |x| < 1,
# which holds since the partial products converge uniformly on any compact subset of the open
# unit disc. Moreover, the limit function has no zeros in the open unit disc
# by Hurwitz's theorem. Finally U(x) has a singularity at 1, so it does not
# converge in any larger open disc centered at 0.
#
# Thus, the above equation for U(x) also holds as a statement about holomorphic functions on the unit disc, and the generating function
# U(x) also exists as a holomorphic function as well, which has no zeroes for
# any complex numbers |x| < 1.
## Bounds and asymptotics for w[n] ##
# In other news, the sequence a[n] = w[n-1] must satisfy
# a[n] <= Fib[1+log[2](n)] < C*n^0.7, where Fib[..] is the Fibonacci sequence,
# since if A[n] = max(a[k], k=1..n), then we have
# A[n] <= Fib[k+1] A[n/2^k] for all k, and A[1] = 1.
# This inequality follows from estimates such as:
# If n is odd, then a[n] = a[(n-1)/2] + a[(n+1)/2], which equals either 2 a[(n-1)/4] + a[(n+3)/4] or a[(n-1)/4] + 2 a[(n+3)/4]
# depending on whether (n-1)/2 or (n+1)/2 is even. Continuing this argument
# gives the above bound involving Fibonacci numbers.
# If n is even, then we repeat with a[n/2] etc. until one of these is odd.
#
# But on the other hand, since the k^th generation of numbers contains the Fibonacci number Fib[k+1], (proven by an easy induction argument)
# it follows that the above upper bound is best possible in the sense that there
# are infinitely many n for which the inequality is tight.
#
# As corollaries, the largest numerator in the k^th generation is exactly Fib[k+1].
# Moreover, A[n] is asymptotic to Fib[log[2](n)+1], which is something like C*n^(0.69...) + o(1).