#Yusra Naqvi: HW #14 #OK to post ################################################################################ #1 ### #GuessH1(L,ord,deg,n,N): inputs a sequence L, non-negative integers ord and deg, #and symbols n and N, and outputs a linear recurrence operator with polynomial #coefficients conjecturally annihilating the sequence L GuessH1:=proc(L,ord,deg,n,N) local ope,a,i,j,var,eq,n1,var1: if nops(L)<(1+ord)*(1+deg)+ord+6 then print(cat(`We need at least `,(1+ord)*(1+deg)+ord+6,` terms in L.`)): return(FAIL): fi: ope:=add(add(a[i,j]*n^i*N^j,i=0..deg),j=0..ord): var:={seq(seq(a[i,j],i=0..deg),j=0..ord)}: eq:={seq( add(subs(n=n1,coeff(ope,N,j))*L[n1+j],j=0..ord), n1=1..(1+ord)*(1+deg)+6)}: var1:=solve(eq,var): ope:=subs(var1,ope): if 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)): eq:=expand({seq( add(subs(n=n1,coeff(ope,N,j))*L[n1+j],j=0..ord), n1=1..nops(L)-ord)}): if eq<>{0} then return(FAIL): else ope: fi: end: ### #GuessH(L,n,N): inputs a sequence L, and symbols n and N, and #outputs a minimal (in the sense of order+degree) linear recurrence operator #with polynomial coefficients conjecturally annihilating the sequence L #with nops(L)>=(1+ord)*(1+deg)+ord+6 GuessH:=proc(L,n,N) local K,ope,ord,deg: ord:=1: deg:=0: for K from 1 while nops(L)>=(1+ord)*(1+deg)+ord+6 do for ord from 1 to K do deg:=K-ord: ope:=GuessH1(L,ord,deg,n,N): if ope<>FAIL then return(ope): fi: od: od: FAIL: end: ### #EmpiricalZeilberger(F,n,k,N,M): inputs a positive integer M, a symbol N, and #an expression F in the letters n and k that is a summand #of a binomial coefficient sum of the form #binomial(n,k)^r*(ProductsOrQuotientOfTermsOfTheForm (a*n+b*k+c)!) #(a,b integers, c number) #(that are expressible as products/quotients of binomial coefficients of the #form binomial(a*n+b*k,a1*n+b1(k)) #and outputs a conjenctured recurrence satisfied by the sequence #a(n1):=Sum(F(n1,k),k=0..n1); #For example, #EmpiricalZeilberger(binomial(n,k)*2^k,n,k,N,30); #should return #N-3 EmpiricalZeilberger:=proc(F,n,k,N,M) local L,i,j: L:=[]: for i from 1 to M do L:=[op(L),expand(add(subs({n=i,k=j},F),j=0..i))]; od: GuessH(L,n,N): end: ################################################################################ #2 #F=binomial(n,k)^4; #-(4*(4*n+5))*(4*n+3)*(1+n)-(2*(2*n+3))*(3*n^2+9*n+7)*N+(n+2)^3*N^2 #F=binomial(n,k)^2*binomial(n+k,k)^2 #(1+n)^3-(2*n+3)*(17*n^2+51*n+39)*N+(n+2)^3*N^2 #F=binomial(n,k)^3*binomial(n+k,k)^2 #Too long to list here: has degree 9 and order 7. #Coefficients involve large numbers. ################################################################################ #3 #This gives the same output in all the suggested example cases except for: #F=binomial(n,k)^3*binomial(n+k,k)^2 #where the output in this case has degree 16 and order 5. #This is also much faster! ################################################################################ #4 ### #HtoSeq(ope,n,N,ini,M): inputs a linear recurrence operator ope #phrased in terms of the discrete variable n and its corresponding #shift operator N, and ord=degree(ope,N) values for the sequence, and #outputs the first M terms of the sequence. #For example: #HtoSeq(n+1-N,n,N,[1],4); #should yield #[1,2,6,24] HtoSeq:=proc(ope,n,N,ini,M) local ord,ope1,L,n1,k,i: ord:=degree(ope,N): if nops(ini)<>ord then return(FAIL): fi: ope1:=expand(subs(n=n-ord,ope)/N^ord): L:=ini: for n1 from nops(ini)+1 to M do k:=subs(n=n1,coeff(ope1,N,0)): if k=0 then print(cat(`Blows up at `,n1,`.`)): print(cat(`So far, we have `,L)): return(FAIL): fi: L:=[op(L),-add(subs(n=n1,coeff(ope1,N,-i))*L[n1-i],i=1..ord)/k]: od: L: end: ### #FirstMterms(F,n,k,M1,M2): finds the first M1 terms of the sequence #a(n1):=Sum(F(n1,k),k=0..n1); #and then guesses a recurrence and uses it to output #the first M2 terms of the same sequence. FirstMterms:=proc(F,n,k,M1,M2) local N,ope,ini,i: ope:=EmpiricalZeilberger(F,n,k,N,M1): ini:=[]: for i from 1 to degree(ope,N) do ini:=[op(ini),expand(add(subs({n=i,k=j},F),j=0..i))]; od: HtoSeq(ope,n,N,ini,M2): end: ### #This gives us: #FirstMterms(binomial(n,k),n,k,20,1000)[1000]; #107150860718626732094842504906000181056140481170553360744375038837035105112493 #612249319837881569585812759467291755314682518714528569231404359845775746985748 #039345677748242309854210746050623711418779541821530464749835819412673987675591 #65543946077062914571196477686542167660429831652624386837205668069376 #FirstMterms(binomial(n,k)^2,n,k,20,1000)[1000]; #204815162698948971433516250298082504439642488798139703382038263767174818620208 #375582893299418261020620146476631999802369241548179800452479201804754976926157 #856301289663432064714851152395251651227768588611539546256147907378668464154444 #533617613770073855673814589630071306510455959514479888746206368718514551828551 #173166276253663773084682932255389049743859481431755030783796444370810085163724 #827462791417016619883764840843541430817785947037746565188475514680749694674923 #803033101818723298009668567458560252549910118113525353465888794196665367490451 #1306110096311906270342502293155911108976733963991149120 #FirstMterms(binomial(n,k)^4,n,k,20,1000)[1000]; #105810444912123795746231188557055866350850558228633681091593665418370891494186 #275699363465887734946955573051849447639465543296792349845753260795072348672281 #435244339445003035470626981308227036571425935313621568747047661757164581052168 #975616577819084529082449013216134436958891223182786100322813203762145710894630 #450188413914920375083595574265703359026647551016680914812050327692150253932517 #553517682754328737817391626069875857890121286766293006090656371389265055516027 #461061014824692725429540752107600345915254195328546681146533848930777668240326 #735505658783193375415535837739494931868932224703232279412455484853533268921741 #004604852764742330532840806455047020202486696009448291836958936259156477010359 #581246025046926254548470278218976807808650103588050982636548700042842293937137 #848051390625872102175960282257212742133546253368076535931728158394219114048806 #842179920660399068531367371391904541836071070807649588481957206221580439672161 #712514061734033231350845783017471436709709958415079454130204505015942750969106 #879666330319919439079121140789444100878283595892393949012927409847571897262243 #425940707180177270366024748975695858386958047762220607527688791414832542037054 #612552185462380343555329821248 ################################################################################ #5 #qGuessH(L,n,N,q): Inputs a sequence L, of polynomials or rational functions #in q and tries to conjecture a recurrence operator ope(q^n,q,N) #annihilating L. qGuessH:=proc(L,n,N,q) local K,ope,ord,deg: ord:=1: deg:=0: for K from 1 while nops(L)>=(1+ord)*(1+deg)+ord+6 do for ord from 1 to K do deg:=K-ord: ope:=qGuessH1(L,ord,deg,n,N,q): if ope<>FAIL then return(ope): fi: od: od: FAIL: end: ### #qGuessH(L,n,N,q): Inputs a sequence L, of polynomials or rational functions #in q and tries to conjecture a recurrence operator ope(q^n,q,N) with #order ord and degree deg, annihilating L. qGuessH1:=proc(L,ord,deg,n,N,q) local ope,a,i,j,var,eq,n1,var1: if nops(L)<(1+ord)*(1+deg)+ord+6 then return(FAIL): fi: ope:=add(add(a[i,j]*q^(n*i)*N^j,i=0..deg),j=0..ord): var:={seq(seq(a[i,j],i=0..deg),j=0..ord)}: eq:={seq( add(subs(n=n1,coeff(ope,N,j))*L[n1+j],j=0..ord), n1=1..(1+ord)*(1+deg)+6)}: var1:=solve(eq,var): ope:=subs(var1,ope): if 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)): eq:=expand({seq( add(subs(n=n1,coeff(ope,N,j))*L[n1+j],j=0..ord), n1=1..nops(L)-ord)}): if eq<>{0} then return(FAIL): else ope: fi: end: ################################################################################ #6 ### #LIF(P,u,N): inputs an expression P(u) that possesses a Maclaurin expansion in u #and outputs the first N terms of the series expansion of the unique f.p.s u(t) #satisfying the functional equation u(t)=t*P(u(t)). #For example: #LIF(1+u^2,u,12); #should return #[1,0,1,0,2,0,5,0,14,0,42,0] LIF:=proc(P,u,N) local n: [seq(coeff(taylor(P^n,u=0,n),u,n-1)/n,n=1..N)]: end: ### #OTN(S,N): inputs a set of positive integers, S, and outputs the list L of #length N of non-neg. integers such that L[n] is the number of ordered trees #with n vertices where every vertex must have a number of children that is a #member of S or 0. #Note that this differs from the statement of problem in that we always take 0 #to be one of the possible numbers of children in order to get finite trees. OTN:=proc(S,N) local S1,P,u,s: S1:=S union {0}: P:=0: for s in S1 do P:=P+u^s: od: LIF(P,u,N): end: ### #TreeH(S,N1,n,N): inputs a set of non-negative integers, S, and conjectures #a linear recurrence operator, ope(n,N) annihilating the enumerating sequence #for the number of ordered trees with n vertices where every vertex must have #a number of children that is a member of S union {0}. TreeH:=proc(S,N1,n,N) local L: L:=OTN(S,N1): GuessH(L,n,N): end: ### #We get: #TreeH({1,2},40,n,N); #-3*n+(-3-2*n)*N+(3+n)*N^2 ################################################################################ #7 #DiagRec(F,x,y,z,M,n,N): inputs a rational function F of the variables x,y,z #(whose denominator does not vanish at (x,y,z)=(0,0,0)), and a positive integer #M (indicating how many terms to use for guessing), and symbols n and N and #outputs a linear recurrence operator with polynomial coefficients ope(n,N) #annihilating the sequence a(n):=coeff of x^n*y^n*z^n in the Taylor expansion #of F around the origin. #For example, #DiagRec(1/(1-x-y-z),x,y,z,40,n,N); #should return #(n+1)^2*N-3*(3*n+1)*(3*n+2) DiagRec:=proc(F,x,y,z,M,n,N) local L: L:=[seq(coeftayl(F,[x,y,z] = [0,0,0], [i,i,i]),i=1..M)]: print(L): GuessH(L,n,N): end: ################################################################################ #8 #F=1/(1-x-y-z+4*x*y*z); #-8*(n+1)^2+(-16-21*n-7*n^2)*N+(n+2)^2*N^2 #Sequence A172 in OEIS #The recurrence is also given. #F=1/(1-x-y-z+x*y+x*z+y*z); #(3*(3*n+4))*(3*n+2)+(n+2)^2*N^2 #If you remove the alternating 0's in the sequence, and take absolute values, #this is Sequence A6480 in OEIS. ################################################################################