#!/usr/local/bin/maple # -*- maplev -*- # Nathaniel Shar # HW 4 # Experimental Mathematics # It is okay to link to this assignment on the course webpage. Help := proc(): print(`Di(L), NGP(L, x), SumPol(P, x), F(a,b), G(a,b), GGuess(a,b), Gk(k,a,b), Gk1k2(k1, k2, a, b), F3(a,b,c), G3(a,b,c), G3Guess(a,b,c)`): end: Di := proc(L): [seq(L[i] - L[i-1], i=2..nops(L))]; end: # Find the unique polynomial of degree at most nops(L)-1 that fits L. NGP := proc(L, x) local i, d, coeffs, diff: d := nops(L)-1: if d = FAIL then return FAIL: fi: diff := L[1..d+1]: coeffs := []: for i from 0 to d do: coeffs := [op(coeffs), diff[1]]: diff := Di(diff): od: return factor(expand(add(coeffs[i]*binomial(x, i-1), i=1..d+1))): end: # Return the polynomial expression for the indefinite sum add(P(i), i=0..n). SumPol := proc(P, x) local d, L, S, i, j: d := degree(P, x): L := [seq(subs(x=i, P), i=0..d+1)]: S := [seq(add(L[i], i=1..j), j=1..nops(L))]: factor(expand(NGP(S, x))): end: # Number of walks in manhattan F := proc(a, b) option remember: if a <= 0 or b <= 0 then: return 1: else: return F(a-1, b) + F(a, b-1): fi: end: # Number of walks south of Broadway. Note: "east" = positive a # direction, "north" = positive b direction. G := proc(a, b) option remember: if b = 0 then: return 1: elif a < b then: return 0: else: return G(a-1, b) + G(a, b-1): fi: end: ############# # Problem 2 # ############# # We conjecture that G(a,b) = GGuess(a,b): GGuess := proc(a,b) option remember: if (type(a, numeric) and type(b, numeric)) then: if b < 0 then: return 0: elif a < b then: return 0: fi: fi: return (1/b!)*(a-b+1)*product(a+i, i=2..b): end: # Now all we have to do is run # simplify(GGuess(a-1,b) + GGuess(a, b-1) - GGuess(a, b)); (1) # This returns 0, which completes the proof. ############# # Problem 3 # ############# # Number of random walks south of a-b = -k. Gk := proc(k, a, b) option remember: if a = 0 and b = 0 then: return 1: elif b < 0 then: return 0: elif a < 0 then: return 0: elif a < b-k then: return 0: else: return Gk(k,a-1,b) + Gk(k,a,b-1): fi: end: # It is easy to see that Gk(1,a,b) = G(a+1, b). Why? Because the # number of ways to reach (a,b) from (0,0) staying south of a-b=-1 is # the same as the number of ways to reach (a,b) from (-1, 0) staying # south of a-b=-1. But then, shifting Manhattan one unit to the east , # we see that this is the same as the number of ways to reach (a+1, b) # from (0,0) staying south of a-b=0. # Similarly, consider Gk(2,a,b). Let N(a,b) be the number of ways to # reach (a,b) from (0,0) staying south of a-b=-2. Let M(a,b) be the # number of ways to reach (a,b) from (-1, 0) under the same # constraint. Now, from (-1,0) we could either start with a step east # or a step north. If we start by going east then there are N(a,b) # ways to reach (a,b). If we start by going north then we are forced # to go east next. At this point, shifting the entire grid south and # east one unit, we see that there are Gk(1,a,b-1) ways to reach # (a,b). Also, by shifting the entire grid east one unit, we see that # M(a,b) = Gk(1, a+1, b). # Combining these facts, we have Gk(1,a+1,b) = Gk(1,a,b-1) + N(a,b), # so that N(a,b) = Gk(1,a+1,b) - Gk(1,a,b-1) # = G(a+2,b) - G(a+1,b-1) # This can be written: # (1/b!)*(a-b+3)*product(a+2+i, i=2..b) - # (1/(b-1)!)*(a-b+2)*product(a+1+i, i=2..b-1) # Simplifying (using Maple) gives: # Gk(2, a, b) = (1/b!)(a+3-b)(a^2 + 3a + ab + b^2 + 2)(a+b)!/(a+3)! ############# # Problem 4 # ############# Gk1k2 := proc(k1, k2, a, b) option remember: if a = 0 and b = 0 then: return 1: elif b < 0 then: return 0: elif a < 0 then: return 0: elif a < b-k1 then: return 0: elif a > b+k2 then: return 0: else: return Gk1k2(k1,k2,a-1,b) + Gk1k2(k1,k2,a,b-1): fi: end: # When k = 1, Gk1k2(k,k,i,i) is sequence A000079 in OEIS. # The formula for the ith term is 2^i. This is obvious # since there are exactly 2 ways to get from (i,i) to (i+1, i+1). # When k = 2, Gk1k2(k,k,i,i) is sequence A025192 in OEIS. # The formula for the ith term is 2*3^(n-1), except that the 0th term # is 1. This is also fairly easy to prove. # When k = 3, Gk1k2(k,k,i,i) is sequence A006012 in OEIS. # A formula for the ith term is not immediately obvious. However, we # may guess one. # Since we expect the sequence to be produced by a linear recurrence # of moderate order (3 or less), we may try to find the largest # eigenvalue of the recurrence by looking at the ratio of consecutive # terms as i increases. This ratio appears to approach 3.4142135... = # 2 + sqrt(2). # We next divide the ith term by (2+sqrt(2))^i. This gives a sequence # that rapidly approaches 1/2. Therefore, the coefficient of # (2+sqrt(2))^i is 1/2. # Now, we could subtract this from the sequence and again take ratios to # find the next largest eigenvalue. Unfortunately, the floating point # numbers start to lose precision at this stage. Instead, it is easier # to just note that 2 - sqrt(2) must also be an eigenvalue of the # recurrence (since it should have rational coefficients), and since # the 0th term is 1 the coefficient of (2-sqrt(2))^i must also be 1/2. # We now guess that the formula for the ith term is # (1/2)((2+sqrt(2))^i + (2-sqrt(2))^i). # It is easy to check that this does in fact produce the sequence. # We can also write down the conjectured recurrence relation. Since # (x-2+sqrt(2))(x-2-sqrt(2)) = x^2 - 4x + 2, the recurrence relation # should be # a_i = 4a_{i-1} - 2a_{i-2}. ############# # Problem 5 # ############# # Number of walks in 3D Manhattan F3 := proc(a,b,c) option remember: if a = 0 and b = 0 and c = 0 then: return 1: elif a < 0 or b < 0 or c < 0 then: return 0: else: return F3(a-1,b,c) + F3(a,b-1,c) + F3(a,b,c-1): fi: end: ############# # Problem 6 # ############# # Number of walks in 3D Manhattan remaining in the region a >= b >= c G3 := proc(a,b,c) option remember: if a = 0 and b = 0 and c = 0 then: return 1: elif a < 0 or b < 0 or c < 0 then: return 0: elif a < b or a < c or b < c then: return 0: else: return G3(a-1, b, c) + G3(a, b-1, c) + G3(a, b, c-1): fi: end: ############# # Problem 7 # ############# # Doing a few examples, it quickly became clear that the result could # be expressed in the form 1/C * P(x), where P(x) is a monic # polynomial with b+c distinct roots and C is an integer. Furthermore, # it also became clear that the b+c factors of P(x) were (a-b+1), # (a-c+2), and a bunch of other factors of the form (a+3), (a+4), ..., # (a+b+c). # Getting the constant C was harder. The values, for small b and c, # are: # # b -> # c 1 2 3 4 5 6 # # |1 2 3 8 30 144 840 # V # 2 12 24 80 360 2016 # # 3 144 360 1440 7560 # # 4 2880 8640 40320 # # 5 86400 302400 # OEIS says the main diagonal is of the form c!(c+1)!, and the next # diagonal is of the form c!*(c+2)!/2. It isn't so hard now to guess # that in general, the diagonals are c!(b+1)!/(b-c+1). # So we guess the following: G3Guess := proc(a,b,c): return (a-b+1)*(a-c+2)*product(a+3+i, i=0..b+c-3)/(c!*(b+1)!/(b-c+1)); end: # And notice that running # # simplify( # G3Guess(a,b,c) # - G3Guess(a,b,c-1) - G3Guess(a,b-1,c) - G3Guess(a-1,b,c) # ); # # yields 0, as desired. (We still have to check the boundary # conditions, but those are easy. For example, we check that # G3Guess(a, a+1, c) = G3Guess(a, b, b+1) = 0.)