### Homework 9 - February 17 2012 - Patrick Devlin ### ### I have no preferrence about keeping my work private ### Help:=proc(): print(`From homework 9:`): print(`T(n,x), testProblem1(howFar:=100), SeqtoGF(L,q)`): print(`From earlier stuff:`): print(`GenHomPolG(X, d, a, DegList:=[infinity], co:=0), GPg(X, d, a, DegList:=[infinity], co:=0)`): print(`GuessPR1(L, x, ORDER, DEGREE)`): print(`RecToPR(C, maxDeg, maxOrd, f, n)`): print(`GP(X, d, a), GuessPR(L, x, ORDER, DEGREE)`): print(`GuessRec(L), GuessRec1(L,r)`): print(`GFtoSeq(f,q,m)`): end: ############# # Problem 1 # ############# # Running GuessPR1(L, x, ORDER, DEGREE) for L:=RecToSeq([[-2*x, 1], [1,2*x]], 50) yielded the following recurrence T:=proc(n,x) option remember: if(n <= 1) then return 1: fi: if(n = 2) then return 4*x^2: fi: return ((4*x^2-2)*T(n-1,x) - T(n-2,x) +2): end proc: # This tests the above recurrence against what it should be (it works) testProblem1:=proc(howFar:=100) local L1, L2, i: L1:=RecToSeq([[-2*x,1],[1,2*x]], howFar): L1:=[seq(L1[i]^2, i=1..nops(L1))]: L2:=[seq(T(i,x), i=1..howFar)]: for i from 1 to howFar do if(not testeq(L1[i], L2[i])) then return false: fi: od: return true: end proc: ############# # Problem 2 # ############# # This takes a list and tries to find a C-finite sequence that fits it. # If this is sucessful, it returns a rational generating function for that sequence. # It returns FAIL if this does not work SeqtoGF:=proc(L, q) local C, denom, i, GF, co: C:=GuessRec(L): if(C = FAIL) then return FAIL: fi: denom:=1: for i from 1 to nops(C[1]) do denom:=denom+q^(i) * C[1][i]: od: if (testeq(denom, 0)) then return FAIL: fi: GF:=q/denom: for i from 1 to nops(C[2]) do co:=coeff(series((C[2][i] * q^i / denom) - GF, q, nops(C[2])+i+10), q^i): GF:=GF+co*q^i/denom: od: return GF, C: end proc: ###################### # Stuff from Earlier # ###################### GFtoSeq:=proc(f,q,m) local i, S: S:=series(f,q,m+5): return [seq(coeff(S,q^i), i=1..m)]: end proc: 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: 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 20 for good measure L:=RecToSeq(C,(maxNumberOfCoeffs * 2 + 20)): # 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: 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: ########### # Class 9 # ########### # GenHomPolG(X,d,a,DegList, co): # ## Inputs: # (i) a list of variables X e.g. [x,y,z] (let m=nops(X), the number of variables) # (ii) a pos. integer d, for the degree # (iii) a symbol a by which the coeff. (undetermined) of the pol. are expressed # (iv) a list of degrees for the individual variables (corresponding to X) # (v) co, starting counter # ## Outputs: # (i) A generic HOMOG. polynomial in the variables X of degree d # with coeff. expressed as a[co],a[co+1], ..., a[co+AsNeeded] # and with the restriction that the degree of X[i] <=DegList[i] # (ii) the set of coefficients # (iii) the new value of the counter # # For example, GenHomPolG([x,y],3,a,[2,2],0) --outputs--> a[0]*x^2*y + a[1]*x*y^2, {a[0],a[1]}, 2 GenHomPolG:=proc(X,d,a,DegList:=[infinity],co:=0) local P, var, co1, m, i, x, X1, P1, newDegList: option remember: m:=nops(X): newDegList:=DegList: if(m>nops(DegList)) then newDegList:=[op(DegList), infinity\$m]: fi: # or maybe we wanna return FAIL instead? if(convert(newDegList,`+`) < d) then return (0, {}, co): fi: x:=X[1]: if m=1 then RETURN(a[co]*x^d, {a[co]},co+1): # We don't need to wonder if d > DegList[1] by that above check fi: co1:=co: X1:=[op(2..m,X)]: newDegList:=[op(2..m, newDegList)]: P:=0: var:={}: for i from 0 to min(d,DegList[1]) do P1:=GenHomPolG(X1, d-i, a, newDegList, co1): P:=expand(P+P1[1]*x^i): var:=var union P1[2]: co1:=P1[3]: od: return (P,var,co1): end proc: #GPg(X,d,a, DegList:=[infinity]): ## Inputs: # (i) a list of variables X e.g. [x,y,z] (let m=nops(X), the number of variables) # (ii) a pos. integer d, for the degree # (iii) a symbol a by which the coeff. (undetermined) of the pol. are expressed # (iv) an optional parameter for DegList, which is a list of maximum degrees for the variables # (v) an optional parameter for where the index of the coefficients should start # ## Outputs: # (i) A generic polynomial in the variables X of total degree d with coeff. expressed as a[co],ac[co+1], ..., a[co+AsNeeded] # (ii) the set of coefficients # (iii) the new value of the counter # # For example, GenPol([x],1,a,0) --outputs--> a[0]+a[1]*x,{a[0],a[1]},2 GPg:=proc(X,d,a, DegList:=[infinity], co:=0) local P, var,co1,d1,P1: co1:=co: P:=0: var:={}: for d1 from 0 to d do P1:=GenHomPolG(X,d1,a,DegList,co1): var:=var union P1[2]: co1:=P1[3]: P:=P+P1[1]: od: return (P,var): end proc: # GuessPR1(L,x,ORDER,DEGREE) ## Inputs # (i) a sequene L of numbers # (ii) a SYMBOL x (for the indexed variables x[1],x[2], ...) # (iii) a pos. integer, ORDER # (iv) a pos. integer Degree ## Output # A polynomial, let's call it, P(x[0],x[1], ..., x[ORDER]) such # that for all n P(L[n],L[n+1], ..., L[n+ORDER])=0 # and the degree of X[ORDER] is 1 GuessPR1:=proc(L,x,ORDER,DEGREE, offset:=0) local X,i,P,var,eq,a: X:=[seq(x[i+offset],i=0..ORDER)]: P:=GPg(X,DEGREE,a,[DEGREE\$ORDER, 1]): # Note, this is the only line that was changed from GuessPR [below] var:=P[2]: if nops(var)+6>nops(L) then print(`Please make the list longer, we need`, nops(var)+6 , `terms`): RETURN(FAIL): fi: P:=P[1]: eq:={ seq( subs( {seq(x[i+offset]=L[n+i],i=0..ORDER)}, P)=0, n=1.. nops(var)+6 ) }: var:=solve(eq,var): factor(subs(var,P)): end: 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:=[GPg(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-1))], 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 Class 6 # ###################### # 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: