### Homework 13 - Patrick Devlin - March 4 2012 ### ## I have no preference about keeping my work private ## Help:=proc(): print(`GuessH(L, n, N), OTN(S,N), AllOTN(S,N)`): print(`From class 13:`): print(`LIF(p,u,N)`): print(`GuessHGivenDegAndOrd(L,DEG,ORD,n,N,printErrors:=false)`): print(`-- Type HelpAll() to see all functions --`): end: HelpAll:=proc(): print(`From class 12:`): print(`AlgToSeq(F, P, x, K)\n`): print(`From class 10:`): print(`GuessAlg(L, P, x, MaxDegP:=(nops(L)/2 - 10)), GuessAlgHelp(L, P, x, degP, maxDegX:=infinity)`): print(`GuessAlgExactDegrees(L, P, x, degP, degX)\n`): print(`From class 9:`): print(`GenHomPolG(X, d, a, DegList:=[infinity], co:=0), GPg(X, d, a, DegList:=[infinity], co:=0)`): print(`GuessPR1(L, x, ORDER, DEGREE)`): print(`From previous classes:`): print(`GenPol(X, d, a, co), GenHOMPol(X, d, a, co)`): print(`GP(X, d, a), GuessPR(L, x, ORDER, DEGREE)`): end: ############# # Problem 1 # ############# GuessH:=proc(L, n, N) local orderPlusDeg, order, temp: for orderPlusDeg from 0 to nops(L) do for order from 0 to orderPlusDeg do temp:=GuessHGivenDegAndOrd(L, order, orderPlusDeg-order, n, N): if(temp<>FAIL) then return temp: fi: od: od: return FAIL: end proc: ############# # Problem 2 # ############# OTN:=proc(S,N) local u: return LIF(add(u^s, s in S), u, N): end proc: ############# # Problem 3 # ############# with(combinat): AllOTN:=proc(S,N) local K, L, S1: K:=subsets(S): L:={}: while not K[finished] do S1:=K[nextvalue](): L:=L union {[S1, OTN( {0} union S1, N)]}: od: return L: end proc: ############# # Problem 4 # ############# # The subsets not in Sloane are: # # {1, 5}, {2, 4}, {2, 5}, {3, 4}, {3, 5}, {4, 5}, {1, 2, 4}, {1, 2, 5}, # {1, 3, 5}, {1, 4, 5}, {2, 3, 4}, {2, 3, 5}, {2, 4, 5}, {3, 4, 5}, # {1, 2, 3, 5}, {1, 2, 4, 5}, {1, 3, 4, 5}, {2, 3, 4, 5}. # # Of these, {1, 5} is the MOST interesting, and even still, I personally don't # think it is "interesting enough" to be in OEIS. #################### # From other stuff # #################### LIF:=proc(p, u, N) local n, s: s:=taylor(p,u,N+2): return [seq(coeftayl(s^n, u=0, n-1)/n,n=1..N)]: # not even checking if p is in fact a function in u holomorphic in the neighborhood u=0. end proc: ## This takes a list, L, and tries to output a linear recurrence operator with polynomial coefficients # of order exactly ORD and maximum degree exactly DEG that annihilates L. # It's in terms of the shift operator N and polynomials n. # ## It will output something like (a_0 (n) + a_1 (n) N + a_2 (n) N^2 + ... + a_ORD (n) N^ORD) L = 0 # or FAIL if it fails. # If the polynomials a_0 (n) = a_1 (n) = ... = a_ORD (n) = 0 are the only solution, then it will return FAIL GuessHGivenDegAndOrd:=proc(L, ORD, DEG, n, N, printErrors:=false) local ope, a, i, j, vars, eqns, n1, solvedVars: if(nops(L) < (1+ORD)*(1+DEG)+ORD+6) then if(printErrors) then print(`List is too short! Need at least `, (1+ORD)*(1+DEG)+ORD+6, ` terms in L`): fi: return FAIL: fi: ope:=add(add(a[i,j]*n^i * N^j, i=0..DEG),j=0..ORD): vars:={seq(seq(a[i,j], i=0..DEG),j=0..ORD)}: eqns:=[seq( add(eval(coeff(ope,N,i)*L[n1+i],{n=n1,N=1}), i=0..ORD) = 0, n1=1..min((1+DEG)*(1+ORD)+6, (nops(L)-ORD)))]: #this "min" is actually just for good measure. We already checked this at first. # We don't solve for ALL the equations cause that would be solving perhaps a LOT of stuff if L is large. solvedVars:=solve(eqns, vars): ope:=subs(solvedVars, ope): if(testeq(ope,0)) then return FAIL: fi: eqns:=[seq( # after passing first "fast" test, we check ALL the data to make sure it worked add(eval(coeff(ope,N,i)*L[n1+i],{n=n1,N=1}), i=0..ORD) = 0, n1=1..(nops(L)-ORD))]: solvedVars:=solve(eqns, vars): ope:=subs(solvedVars, ope): if(testeq(ope,0)) then return FAIL fi: ope:=subs({seq(v=1, v in var)}, ope): ope:=add(factor(coeff(ope,N,j))*N^j, j=0..degree(ope,N)): return ope: end proc: ############ # Class 12 # ############ ### AlgToSeq(F,P,x,K) ## Inputs: # (i) an algebraic expression F in variables P and x, where P is a the gen func. of a mystery sequence in the variable x # (ii) P, the generating function # (iii) an indeterminant x that is the dummy variable in the gen function P # (iv) a positive integer K ## Outputs: # A sequence of length K+1 denoting the coefficient of x^0 , x^1 , ... , x^K in the Maclaurin expansion of P # - If this fails, it outputs FAIL. - # ## Examples: # AlgToSeq(P*(1-x)-1, P, x, 4) --outputs--> [1,1,1,1,1] # AlgToSeq(P-1-x*P^2, P, x, 4) --outputs--> [1,1,2,5,14] AlgToSeq:=proc(F,P,x,K,start:=1) local H, G, i, F1, f: H:=coeff(F,P,1): if(testeq(0, eval(H,x=0))) then return FAIL: fi: # Need to check this or the method will fail # Think...! The expression F says... P = (P*H-F)/H. So set this F1:=(P*H-F)/H: # Then F1 = -(a0 (x) + a2(x)*P^2 + ...) / H(x). # So if we iterate this, the P's ought to go further and further away, and we'll be left with F1(x) = STUFF(x) + P^BIG * OTHERSTUFF(P,x) f:=start: for i from 1 to K do f:=expand(eval(F1, P=f)): f:=series(f,x=0,2*K+10): # truncate it for speed od: f:=series(f,x=0,K+1): return [seq(coeff(f,x,i), i=0..K)]: end proc: ############ # Class 10 # ############ # GuessAlg(L, P, x, MaxDegP:=(nops(L)/2 - 10)): Takes an input sequence of numbers L and # tries to guess and algebraic equation of the form: # F(P,x), where P is the symbol for the generating function of L, [as best as we can tell!], # and x is an indeterminant (corresponding to generating function for L). GuessAlg:=proc(L, P, x, MaxDegP:=(nops(L)/2-10), MaxDegX:=infinity) local g, DegP: for DegP from 1 to MaxDegP do g:=GuessAlgHelp(L, P, x, DegP, MaxDegX): if(g <> FAIL) then return g: fi: od: return FAIL: end proc: GuessAlgHelp:=proc(L, P, x, DegP, MaxDegX:=infinity) local g, DegX: for DegX from 0 while( ((1+DegX)*(1+DegP) < nops(L) - 6) and (DegX <= MaxDegX)) do g:=GuessAlgExactDegrees(L, P, x, DegP, DegX): if(g <> FAIL) then return g: fi: od: return FAIL: end proc: ### ## Given a list L, formal symbols, P, and x, and two integers, DegP, DegX, ## it tries to return a polynomial of the form ## F(P, x) [not equal 0], with max degree of P equal to DegP, and max degree of x equal to DegX ## such that if genP(x) is the generating function of L, then F(genP, x) is close enough to 0. GuessAlgExactDegrees:=proc(L, P, x, DegP, DegX) local F, a, i, j, vars, genP, FofP, eqns, solvedVars, v: F:=add(add(a[i,j]*x^i * P^j, i=0..DegX), j=0..DegP): vars:={seq(seq(a[i,j],i=0..DegX), j=0..DegP)}: genP:=add(L[i] * x^(i-1), i=1..nops(L)): FofP:=series(eval(F, P=genP), x, nops(L)+10): eqns:={ seq( coeff(FofP, x, i) = 0, i=0..nops(L)-1)}: solvedVars:=solve(eqns, vars): F:=eval(F, solvedVars): if (testeq(F,0)) then return FAIL: fi: # This next line is "optional". Dr. Z. put it in to get rid of "stupid looking" free parameters F:=eval(F, {seq(v = 1, v in vars)}): return add(factor(coeff(F,P, i))*P^i, i=0..degree(F,P)): # this is just F cleaned up a little 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) local X,i,P,var,eq,a: X:=[seq(x[i],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]=L[n+i],i=0..ORDER)}, P)=0, n=1.. nops(var)+6 ) }: var:=solve(eq,var): factor(subs(var,P)): end: #GuessPR(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 GuessPR:=proc(L,x,ORDER,DEGREE) local X,i,P,var,eq,a: X:=[seq(x[i],i=0..ORDER)]: P:=GP(X,DEGREE,a): 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( {x[0]=L[n],x[1]=L[n+1], ..., x[ORDER]=L[n+ORDER]}, # P)=0, n=1..min(nops(L)-ORDER, nops(var)+6) }: eq:={ seq( subs( {seq(x[i]=L[n+i],i=0..ORDER)}, P)=0, n=1.. nops(var)+6 ) }: var:=solve(eq,var): factor(subs(var,P)): end: