### Pat Devlin - Homework 10 - February 25 2012 -- I have no preferrence about keeping my work private ### Help:=proc(): print(`DiagToAlg(R, x, y, P, t, MaxDegP), AlgToSeq(F,P,x,K)`): 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)`): print(`Type HelpAll() to see all functions`): end: HelpAll:=proc(): print(`From previous classes:`): print(`GenHomPolG(X, d, a, DegList:=[infinity], co:=0), GPg(X, d, a, DegList:=[infinity], co:=0)`): print(`GuessPR1(L, x, ORDER, DEGREE)`): print(`GenPol(X, d, a, co), GenHOMPol(X, d, a, co)`): print(`GP(X, d, a), GuessPR(L, x, ORDER, DEGREE)`): end: ############# # Problem 2 # ############# ### DiagToAlg(R, x, y, P, t, MaxDegP) # Inputs a function R of the variables x and y, and # tries to guess an algebraic equation F(P, t) with deg(F,P) <= MaxDegP # such that if R = sum_{i,j>=0} a[i,j] x^i y^j, and Di(t):= sum_{i>=0} a[i,i] t^i # then F(Di,t) = 0. # Else it returns fail. DiagToAlg:=proc(R,x,y,P,t,MaxDegP, maxDegX:=infinity): return GuessAlg([seq(coeftayl(R, [x,y]=[0,0], [i,i]), i=0..MaxDegP*10+15)], P, t, MaxDegP, maxDegX): end proc: ## Testing this on the desired inputs yields.. # (1) DiagToAlg(1/((1-x)*(1-y)),x,y,P,t,6) = 1 + (t-1)*P # (2) DiagToAlg(1/(1-x-y),x,y,P,t,6) = 1/4 + (t-1/4)*P^2 # (3) DiagToAlg(1/(1-x-y+3*x*y),x,y,P,t,6) = -1/9 + (t^2 + 2/9*t + 1/9)*P^2 # (4) DiagToAlg(1/(1-x-x^2-y-y^2),x,y,P,t,10) = 1/128 * ( 1 + 2*(4*t-3)*(4*t^2 - 8*t + 1)*P^2 + (8*t + 5)*(4*t^2 - 8*t +1)^2 P^4 ) ############# # Problem 3 # ############# AlgToSeq:=proc(F,P,x,K) local H, G, iteratedResult, iter: H:=coeff(F,P,1): G:=H*P-F: if(testeq(H,0)) then return FAIL: fi: iteratedResult:=1: for iter from 1 to K do try: iteratedResult:=eval(G/H, P=iteratedResult): try: eval(iteratedResult, x=0): catch: return FAIL: end try: iteratedResult:=series(expand(iteratedResult),x,K+5): catch: return FAIL: end try: od: return [seq(coeff(iteratedResult,x,i),i=0..K-1)]: end proc: ###################### # Stuff from earlier # ###################### # 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: