#Homework by Ben Miles for Experimental math assigned 2/13 OK to post GuessPR:=proc(L,x,ORDER,DEGREE) local X,i,P,var,var1,eq,a,Q: X:=[seq(x[i],i=0..ORDER)]: P:=GP(X,DEGREE,a): var:=P[2]: if nops(var)+6>nops(L) then RETURN(FAIL): fi: P:=P[1]: eq:={ seq( subs( {seq(x[i]=L[n+i],i=0..ORDER)},P)=0, n=1.. nops(var)+6 ) }: var:=solve(eq,var): var1:=[]; for i from 1 to nops(var) do Q:=a[i-1]; if type(subs(var[i],Q),rational)<>true then var1:=[op(var1),a[i-1]=1]: else var1:=[op(var1),var[i]]: fi: od: P:=factor(subs(var,P)): P:=subs(var1,P): end: GP:=proc(X,d,a) local P, var,co1,d1,P1: co1:=0: P:=0: var:={}: for d1 from 0 to d do P1:=GenHOMPol(X,d1,a,co1): var:=var union P1[2] : co1:=P1[3]: P:=P+P1[1]: od: P,var: end: GenHOMPol:=proc(X,d,a,co) local P, var, co1,m,i,x,X1,P1: option remember: m:=nops(X): x:=X[1]: if m=1 then RETURN(a[co]*x^d, {a[co]},co+1): fi: co1:=co: X1:=[op(2..m,X)]: P:=0: var:={}: for i from 0 to d do P1:=GenHOMPol(X1,d-i,a,co1): P:=expand(P+P1[1]*x^i): var:=var union P1[2]: co1:=P1[3]: od: P,var,co1: end: RecToSeq := proc(Seq, n) local ini, co, out, i, j: ini := Seq[2]: co := Seq[1]: out := ini: for i from nops(ini)+1 to n do: out := [op(out), add(-co[j]*out[i-j], j=1..nops(co))]: od: return out: end: #New stuff #1. 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,ORDER) local D,P: for D from 0 do: P:=GuessPR(L,x,ORDER,D): if P= FAIL then return(FAIL) elif P<>0 then return(P) fi: od: end: #2. RecToPR(C,DEG,ORD,f,n) 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. RecToPR:=proc(C,DEG,ORD,f,n) local L,k,x,P,i: k:=((DEG+1))^((ORD+1))+10; L:=RecToSeq(C,k); P:=GuessPR(L,x,ORD,DEG); for i from 0 to ORD do: P:=subs(x[i]=f(n+i), P); od: P=0; end: #3. tribonacci tribonacc:=proc(m,k) local i,j,P1,P,L: L=[]: for i from 1 to m do for j from 1 to k do P1:= 0 = 0: P:=RecToPR([[-1,-1,-1],[0,0,1]],j,i,f,n): if P<>P1 then print(P,j,i): else print(j,i): fi: od: od: end: #It appears as though the lowest is RecToPR([[-1,-1,-1],[0,0,1]],3,2,f,n): #4.kbonacci kbonacc:=proc(a,b,k)local i,j,P1,P,L: L=[]: for i from a to a do for j from 1 to b do P1:= 0 = 0: P:=RecToPR([[-1\$k],[0\$(k-1),1]],j,i,f,n): if P<>P1 then return(P,j,i): else print(j,i): fi: od: od: end: #for k=4 it seems to be RecToPR([[-1,-1,-1,-1],[0,0,0,1]],8,3,f,n) # for k=5 it seems to be RecToPR([[-1,-1,-1,-1,-1],[0,0,0,0,1]],5,4,f,n)