### Homework 8. Patrick Devlin. February 14, 2012 ### ### I have no preferrences about keeping my work private ### ### Programs from homework 7 and class 6 at the bottom ### Help:=proc(): print(`From homework 8: findPolySpecialProblem(), GuessMPR(L, X, ORDER), RecToPR(C, DEGREE, ORDER, f, n)`): print(`From class 8: GuessPR(L, X, ORDER, DEGREE)`): print(`From homework 7: GenHomPol(X,d,a,co:=1), GenNonHomPol(X,d,a,co:=1)`): print(`From class 6: GuessRec1(L,r), GuessRec(L) `): print(`From homework 6: RecToSeq(C,n)`): end proc: ############################### # Problem 0 (Special Problem) # ############################### findPolySpecialProblem:=proc() local L, g, tempAns, x, y, A, var, solvedStuff: tempAns:=GenNonHomPol([x, y], 4, A): g:=tempAns[1]: var:=tempAns[2]: L:=[[8, 29, 104, 293, 680], [61, 112, 237, 496, 973], [216, 309, 504, 861, 1464], [557, 704, 989, 1472, 2237], [1192, 1405, 1800, 2437, 3400]]: g:=eval(g, solve( {seq(seq( eval(g,[x=i, y=j])=L[i][j], j=1..nops(L[i])) ,i=1..nops(L))}, var)): with(plots, textplot, display): display(plots[implicitplot](g,x=-6..6,y=-5..5,axes=none, style = POINT, symbol=point, numpoints=10000), textplot([-3.2, 1.5, `proc():`]), textplot([-2.02, 1.15, `print(``I love Experimental Math!``)`]), textplot([-3.2, .8, `g:=`]), textplot([-1.5,.8,g]), textplot([-1.59, .475, `display(plots[implicitplot](g, x=-6..6, y=-5..5, ...`]), textplot([-3.1,.17,`end proc:`])): end proc: ############# # Problem 1 # ############# GuessMPR:=proc(L, X, ORDER) local DegreeCounter, P: for DegreeCounter from 0 to (nops(L)-6) do P:=GuessPR(L, X, ORDER, DegreeCounter): if(P <> FAIL) then return P: fi: od: return FAIL: end proc: ############# # Problem 2 # ############# RecToPR:=proc(C, maxDeg, maxOrd, f, n) local degreeOrderSum, tempDeg, tempOrd, L, maxNumberOfCoeffs, varList, answer: maxNumberOfCoeffs:=binomial(maxDeg+maxOrd+1, maxOrd+1): # see note below for this count, we double it and add 10 for good measure L:=RecToSeq(C,(maxNumberOfCoeffs * 2 + 10)): # we won't need a list bigger than this if(L=FAIL) then return FAIL: fi: for degreeOrderSum from 0 to (maxDeg+maxOrd) do for tempDeg from max(0, degreeOrderSum-maxOrd) to min(degreeOrderSum, maxDeg) do tempOrd:=degreeOrderSum-tempDeg: varList:=[seq(f(n+i), i=0..tempOrd)]: answer:=GuessPR(L, varList, tempOrd, tempDeg, 1): if(answer<> FAIL) then return answer: fi: od: od: return FAIL: end proc: ## Note! # # Let C(D,R) be the number of coefficients in a general polynomial of total degree at most D in at most R variables. # Then C(D,R) is the same as the number of terms x1 ^d1 * x2 ^d2 * ... * # xr^dR, where each di is a nonnegative integer, and d1 + d2 + d3 + ... + dR <= D. # This is the same as how many ways there are to put AT MOST D # indistinguishable objects into R distinguishable containers, # which (if you add a new container called "I'm not using this object right now") is the same as # the number of ways to be EXACTLY D indistinguishable objects into R+1 distinguishable containers. # # By the usual "stars and bars" argument, this gives us C(D,R) = binomial(D+R, R). ### That is to say, let D stars: ****** be the objects, and use R bars | |... | to divide up the objects into the R+1 containers, ### then every ordering looks something like ***|*|**||****| , which would be for instance 3 in box 1, 1 in box 2, 2 in box 3, 0 in box 4, 4 in box 5, and 0 in box 6. ### There are binomial(D+R, R) ways of arranging these stars and bars # # Therefore! Since a recurrence of order ORD needs ORD+1 variables, we know # that the largest possible polynomial will have binomial(maxDeg+maxOrd+1, maxOrd+1) coefficients ############# # Problem 3 # ############# # Running RecToPR on C:=[[-1,-1,-1], [0,0,1]] with various values of degree and order reveals: # # The recurrence with smallest order [then smallest degree] has order 2 and degree 3. # The recurrence with smallest degree [then smallest order] has degree 1 and order 3. # # The recurrence with (degree + order) minimized has degree 1 and order 3 # # The "actual" relation of the tribonacci numbers is most naturally: # [[-1,-1,-1], [1,1,1]], which is to say, f(1) = f(2) = f(3) = 1, and f(n) + f(n+1) + f(n+2) = f(n+3) ############# # Problem 4 # ############# # Honestly, I don't see a pattern other than the (perhaps "trivially obvious"?) pattern that # [[-1, -1, -1, ... , -1], [0, 0, 0, ... , 0, 1]] defines a shift of the sequence # [[-1, -1, -1, ... , -1], [1, 1, 1, ... , 1, 1]] # # This seems to be the most natural recurrence for the k-bonacci numbers ## (It's possible there is some simple-looking quasipolynomial that fits well, but I didn't check that) ###################### # Stuff from Class 8 # ###################### # GuessPR(L,X, ORDER, DEGREE) ## Inputs a sequence of numbers L, a formal symbol X (for indexed variables X[0], X[1], ..., X[ORDER]), and nonnegative integers, ORDER and DEGREE ## Outputs a polynomial P of total degree at most DEGREE in the variables x[0], x[1], ..., x[ORDER] # such that P(L[n], L[n+1], L[n+2], ... , L[n+order]) = 0 for all relevant n # # If it fails to do so, then it returns FAIL GuessPR:=proc(L,X, ORDER, DEGREE, offset:=0) local ans, P, var, x, eqns, iter, k, a: x:=[seq(X[iter+offset], iter=0..ORDER)]: ans:=[GenNonHomPol(x,DEGREE,a,0)]: P:=ans[1]: var:=ans[2]: # The following line tries to solve a linear system that is as big as possible eqns:=solve([seq(0 = eval(P, [seq(X[k+offset]=L[iter+k], k=0..ORDER)] ),iter=1..(nops(L)-ORDER))], var): if(eqns <> NULL) then try: P:=eval(P,eqns): catch: return FAIL: end try: if(P=0) then return FAIL: fi: return factor(P): fi: return FAIL: end: ######################### # Stuff from Homework 7 # ######################### # Takes a list of variables X, the total degree d, a list of coefficients a, # and an optional parameter co, which keeps track of depth in recursive calls. # # A modified version of GenPol. It returns a polynomial in the variables # X[1], ... , X[nops(X)], where each term is a monomial of total degree d. GenHomPol:=proc(X,d,a,co:=1) local P, var, co1,m,i,x,X1,P1: option remember: m:=nops(X): x:=X[m]: if m=1 then RETURN(x^d * a[co], {a[co]},co+1): fi: co1:=co: X1:=[op(1..m-1,X)]: P:=0: var:={}: for i from 0 to d do P1:=GenHomPol(X1,d-i,a,co1): P:=expand(P+x^i * P1[1]): var:=var union P1[2]: co1:=P1[3]: od: return (P,var,co1): end: # Takes a list of variables X, a maximum total degree d, a variable name for # the list of coefficients, a, and an option argument co, indicating which # index of the list of variables should be the first used for the coefficients. # # Outputs a general polynomial in X, whose maximum total degree is exactly d GenNonHomPol:= proc(X,d,a, co:=1) local P, var, i, tempAns, co1: P:=0: var:={}: co1:=co: for i from 0 to d do tempAns:=GenHomPol(X,i,a,co1): P:=expand(P+tempAns[1]): var:=var union tempAns[2]: co1:=tempAns[3]: od: return (P, var, co): end proc: ###################### # Stuff from Class 6 # ###################### # RecToSeq(C,n) takes a C-finite sequence C and returns a list of length n # generated according to the rules dictated by C # # Example: RecToSeq([[-1,-1], [1,1]], 10) --outputs--> # [1,1,2,3,5,8,13,21,34,55] RecToSeq:=proc(C,n) local i, i2, L: if (nops(C) <> 2) then return FAIL: fi: if (nops(C[2]) < nops(C[1])) then return FAIL: fi: # two checks that the input is good L:=[seq(C[2][i], i=1..min(n, nops(C[2])))]: for i from (1+nops(C[2])) to n do L:=[op(L), add(-1*(C[1][j])*(L[i-j]) , j=1..nops(C[2]))]: od: return L: end proc: # GuessRec1(L,r): given a list of numbers and a positive integer, r, # it tries to guess a linear recurrence relation of order r # that fits the list # Looking for a recurrence of the form L[n] + c_1 L[n-1] + ... + c_r L[n-r] = 0 for all relevant n # We need r data points to determine this system of (c_1, c_2, ... , c_{r}) unknowns # So we'll output [ [c_1, c_2, c_3, ... , c_r], [initial_conditions]]. # (i.e., Fibonacci would be F[n] - F[n-1] - F[n-2] = 0 -> [[-1, -1], [1,1]] GuessRec1:=proc(L,r) local var, c, i, n, eqns: if (nops(L) <= 2*r + 6) then return FAIL; fi: # The list isn't long enough, so return fail var:={seq(c[i], i=1..r)}: eqns:=[ seq( (L[n] + sum(c[i] * L[n-i],i=1..r))=0, n=(r+1)..nops(L))]: var:=solve(eqns,var); if(var = NULL) then return FAIL: fi: return [[seq(eval(c[i],var),i=1..r)],L[1..r]]; end proc: # GuessRec(L) returns a linear recurrence relation of minimal order that fits L GuessRec:=proc(L) local r, ans: for r from 1 to (nops(L)-6)/2 do ans:=GuessRec1(L,r): if (ans <> FAIL) then return ans: fi: od: return FAIL: end proc: