# Matthew Russell # Experimental Math # Homework 8 # I give permission to post this. ##################################################################################### # Using GuessPR(L,x,ORDER,DEGREE), write a procedure # GuessMPR(L,x,ORDER) # That finds the (non-zero!) polynomial relation of order ORDER of # smallest possible degree, compatible with the amount of data in L, or returns FAIL. GuessMPR:=proc(L,x,ord) local deg,G: for deg from 1 to nops(L) do G:=GuessPR(L,x,ord,deg): if G<>0 then return G: fi: od: return `FAIL`: end: ##################################################################################### # Using procedure # RecToSeq(C,n) # that you hopefully wrote in hw6.txt, write a procedure # RecToPR(C,DEG,ORD,f,n) # that tries to find a (minimal) polynomial relation of degree DEG and # order ORD, satisfied by the members of the C-finite sequence given by C. For example # RecToPR([[-1,-1],[1,1]],4,2,f,n); # should give # (f(n)^2+f(n)f(n+1)-f(n+1)^2-1) (f(n)^2+f(n)f(n+1)-f(n+1)^2+1)=0 # [Note: I'll be very proud of you if you would figure the number of terms, # n, needed in the function-call to RecToSeq(C,n) to enable guessing for the # given DEG and ORD, using human reasonings, i.e. what is the number of terms # in a generic polynomial with ORD variables and total degree DEG, plus a few terms to spare]. RecToPR:=proc(C,deg,ord,f,n) local m,L,G,i: m:=nops(GP([seq(x[i],i=0..ord)],deg,a)[2])+10+ord: L:=RecToSeq(C,m): G:=GuessPR(L,x,ord,deg): for i from 0 to ord do G:=subs(x[i]=f(n+i),G): od: end: ##################################################################################### # What is the minimal order and degree, and what is the actual relation, satisfied by the Tribonacci numbers? Trib:=proc(n) option remember: if n<2 then return 0: elif n=2 then return 1: else return Trib(n-1)+Trib(n-2)+Trib(n-3): fi: end: # It appears that there is no first-order relation satisfied by the Tribonacci numbers. # Using GuessMPR, I found that # -1+T(n)^3+2*T(n)^2*T(n+1)+T(n)^2*T(n+2)+2*T(n)*T(n+1)^2-2*T(n)*T(n+2)*T(n+1) # -T(n)*T(n+2)^2+2*T(n+1)^3-2*T(n+2)^2*T(n+1)+T(n+2)^3 = 0 # is a second-order, third-degree relation satisfied by the Tribonacci numbers. ##################################################################################### # What is the minimal order and degree, and what is the actual relation, satisfied by # the K-bonacci numbers (the sequence f(n) defined by f(0)=...f(K-2)=0, f(K-1)=1 # f(n)=f(n-1)+...+f(n-K) for n ≥ K), for K=4,5,6. Do you detect a pattern? a4Bon:=proc(n) option remember: if n<3 then return 0: elif n=3 then return 1: else return a4Bon(n-1)+a4Bon(n-2)+a4Bon(n-3)+a4Bon(n-4): fi: end: # I found a third-order relation satisfied by the 4-bonacci numbers: # -(1/2)*a[449]*(-3*x[3]*x[0]*x[1]*x[2]-1-4*x[3]^2*x[0]*x[1]+2*x[3]^2*x[0]*x[2] # -2*x[3]*x[0]^2*x[1]-2*x[3]*x[0]^2*x[2]-3*x[3]*x[0]*x[1]^2-5*x[3]*x[0]*x[2]^2 # -x[1]*x[3]^2*x[2]-2*x[1]^2*x[3]*x[2]-8*x[1]*x[3]*x[2]^2+12*x[0]*x[1]^2*x[2] # +8*x[0]*x[1]*x[2]^2+7*x[0]^2*x[1]*x[2]+2*x[1]*x[3]^3-3*x[1]^3*x[3]+3*x[1]^3*x[2] # +5*x[1]*x[2]^3+7*x[1]^2*x[2]^2+3*x[3]^3*x[2]-2*x[3]^2*x[2]^2-2*x[3]*x[2]^3 # -x[3]^2*x[0]^2+x[3]*x[0]^3+3*x[0]^3*x[1]+2*x[0]^3*x[2]+4*x[0]^2*x[1]^2 # +2*x[0]*x[1]^3+x[0]*x[2]^3+x[3]^3*x[0]+x[0]^4-x[1]^4-x[3]^4+3*x[2]^4)*(-3*x[3]*x[0]*x[1]*x[2] # +1-4*x[3]^2*x[0]*x[1]+2*x[3]^2*x[0]*x[2]-2*x[3]*x[0]^2*x[1]-2*x[3]*x[0]^2*x[2] # -3*x[3]*x[0]*x[1]^2-5*x[3]*x[0]*x[2]^2-x[1]*x[3]^2*x[2]-2*x[1]^2*x[3]*x[2] # -8*x[1]*x[3]*x[2]^2+12*x[0]*x[1]^2*x[2]+8*x[0]*x[1]*x[2]^2+7*x[0]^2*x[1]*x[2] # +2*x[1]*x[3]^3-3*x[1]^3*x[3]+3*x[1]^3*x[2]+5*x[1]*x[2]^3+7*x[1]^2*x[2]^2 # +3*x[3]^3*x[2]-2*x[3]^2*x[2]^2-2*x[3]*x[2]^3-x[3]^2*x[0]^2+x[3]*x[0]^3 # +3*x[0]^3*x[1]+2*x[0]^3*x[2]+4*x[0]^2*x[1]^2+2*x[0]*x[1]^3+x[0]*x[2]^3+x[3]^3*x[0]+x[0]^4-x[1]^4-x[3]^4+3*x[2]^4) # This has degree 8. a5Bon:=proc(n) option remember: if n<4 then return 0: elif n=4 then return 1: else return a5Bon(n-1)+a5Bon(n-2)+a5Bon(n-3)+a5Bon(n-4)+a5Bon(n-5): fi: end: a6Bon:=proc(n) option remember: if n<5 then return 0: elif n=5 then return 1: else return a6Bon(n-1)+a6Bon(n-2)+a6Bon(n-3)+a6Bon(n-4)+a6Bon(n-5)+a6Bon(n-6): fi: end: ##################################################################################### # Old programs ##################################################################################### GenHomPol:=proc(X,d,a,counter) local x,m,counter1,i,P1,X1,P,var: m:=nops(X): counter1:=counter: x:=X[1]: if m=1 then return a[counter1]*x^d,{a[counter1]},counter1+1: else P:=0: var:={}: X1:=[op(2..m,X)]: for i from 0 to d do P1:=GenHomPol(X1,i,a,counter1): P:=expand(P+P1[1]*x^(d-i)): var:=var union P1[2]: counter1:=P1[3]: od: return P,var,counter1: fi: end: GP:=proc(X,d,a) local i,P,var,counter1,G: P:=0: var:={}: counter1:=0: for i from 0 to d do G:=GenHomPol(X,i,a,counter1): P:=P+G[1]: var:=var union G[2]: counter1:=G[3]: od: return P,var: end: # L: a sequence believed to satisfy a polynomial relation # x: symbol for variables x[1], x[2], etc. # ord: order of polynomial relation # deg: degree of polynomial relation # # Output: a polynomial P(x[0],x[1],...,x[order]) such that # for all n, P(L[n],L[n+1],...,L[n+order])=0. GuessPR:=proc(L,x,ord,deg) local X,i,G,P,var,eq,j: X:=[seq(x[i],i=0..ord)]: G:=GP(X,deg,a): P:=G[1]: var:=G[2]: if nops(var)+6+ord>nops(L) then print(`Please make L longer`): return `FAIL`: fi: eq:={seq(subs(seq(x[j]=L[j+i],j=0..ord),P)=0,i=1..nops(var)+6)}: var:=solve(eq,var): return subs(var,P): end: