#John Kim #Use whatever you like #read "C:/Users/John Y. Kim/Documents/Maple/hw16.txt": #old stuff from C13.txt Help13:=proc(): print(` LIF(P,u,N), GuessH1(L,ORD,DEG,n,N)`): end: #Lagrange Inversion applied to tree-enumeration #of ordered #Starting the holonomic ansatz (a.k.a P-recursive) #LIF: Inputs an expression P of the variable #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 eq. #u(t)=t*P(u(t)). For example, #LIF(1+u^2,u,10); should return something with #Catalan numbers LIF:=proc(P,u,N) local n: [seq(coeff(taylor(P^n, u=0,n),u,n-1)/n,n=1..N)]: end: #GuessH1(L,ORD,DEG,n,N): Inputs a sequence of numbers L #and pos. integer ORD and non-neg. integer DEG and letters #n and N and outputs an linear recurrence operator with #poly. coeffs. (conj.) annihilating the sequence L #nops(L)>=(1+ORD)*(1+DEG)+ORD+6 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(`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+DEG)*(1+ORD)+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 print(ope, `didn't work out`): RETURN(FAIL): else RETURN(ope): fi: end: #end old stuff #start new stuff #GuessH(L,n,N): Inputs a sequence of numbers L # and symbols n and N #and finds a minimal (in the sense of ord+deg) # linear recurrence operator with #poly. coeffs. (conj.) annihilating the sequence L with #nops(L)>=(1+ORD)*(1+DEG)+ORD+6 GuessH:=proc(L,n,N) local K,ope,DEG,ORD: 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: #SpellOut(ope,n,N,f) #writes the fact ope(n,N)f(n)=0 in common languae SpellOut:=proc(ope,n,N,f) local i: add(coeff(ope,N,i)*f(n+i),i=0..degree(ope,N))=0: end: #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 outputs the first M terms of #the sequence. For example, #HtoSeq(n+1-N,n,N,[1],4); should yield #[1,2,6,24] #N^(-ORD)*ope(N,n) HtoSeq:=proc(ope,n,N,Ini,M) local L,n1,ORD,BenMiles,kirsten,i: ORD:=degree(ope,N): if nops(Ini)<>ORD then RETURN(FAIL): fi: BenMiles:=expand(subs(n=n-ORD,ope)/N^ORD): L:=Ini: for n1 from nops(Ini)+1 to M do kirsten:=subs(n=n1,coeff(BenMiles,N,0)): if kirsten=0 then print(`blows up at`, n1): print(`so far we have`, L): RETURN(FAIL): fi: L:=[op(L), -add(subs(n=n1,coeff(BenMiles,N,-i))*L[n1-i],i=1..ORD)/kirsten]: od: L: end: #EmpiricalZeilberger(F,n,k,N,M): #inputs, in addition to the positive integer M, 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 be expressible as products/quotients of binomial coefficients if 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 #and EmpiricalZeilberger(binomial(n,k)^2,n,k,N,30); #(n+1)*N-(4*n+2) EmpiricalZeilberger:=proc(F,n,k,N,M) local L,n1,k1: L:=[seq(sum(subs({n=n1,k=k1},F),k1=0..n1),n1=1..M)]: GuessH(L,n,N): end: Help16:=proc(): print(`LIFtoS(F,P,n), NRope(F,P,n,N,K)`): print(`Rope(F,P,n,N) `): end: #LIFtoS(F,P,n): #Inputs an expression F in the variable P. # and a pos. integern #outputs the list of the first n coefficients in the formal power #series that is the UNIQUE solution of the FUNCTIONAL equatation #P(x)=x*F(P(x)) #For example, to get the the enumerating sequence for the #number of complete binary trees with i vertices, up to 10 #do LIFtoS(1+P^2,P,x,10); #P(x)=x*(1+P(x)^2) LIFtoS:=proc(F,P,n) local x,Y,i,p: Y:=coeff(F,P,0)*x: for i from 1 to n do Y:=x*subs(P=Y,F): Y:=add(coeff(Y,x,p)*x^p,p=0..n): od: [seq(coeff(Y,x,i),i=1..n)]: end: ############stuff from AZ.txt############## #AZ.txt: The Almkvist-Zeilberger algorithm #taken from #http://www.math.rutgers.edu/~zeilberg/EKHAD HelpAZ:=proc(): print(` AZd(A,y,n,N) `): end: pashetd:=proc(p,N) local gu,p1,mekh,i: p1:=normal(p): mekh:=denom(p1): p1:=numer(p1): p1:=expand(p1): gu:=0: for i from 0 to degree(p1,N) do gu:=gu+factor(coeff(p1,N,i))*N^i: od: RETURN(gu,mekh): end: goremd:=proc(N,ORDER) local i, gu: gu:=0: for i from 0 to ORDER do gu:=gu+(bpc[i])*N^i: od: RETURN(gu): end: duisd:= proc(A,ORDER,y,n,N) local gorem, conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu, mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,shad,KAMA,i1,apc,apd: KAMA:=10: gorem:=goremd(N,ORDER): conj:=gorem: yakhas:=0: for i1 from 0 to degree(conj,N) do yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A): od: lu1a:=yakhas: LU1:=numer(yakhas): yakhas:=1/denom(yakhas): yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A): P:=LU1: Q:=numer(yakhas): R:=denom(yakhas): j1:=0: while j1 <= KAMA do g:=gcd(R,Q-j1*diff(R,y)): if g <> 1 then Q2:=(Q-j1*diff(R,y))/g: R:=normal(R/g): Q:=normal(Q2+j1*diff(R,y)): P:=P*g^j1: j1:=-1: fi: j1:=j1+1: od: P:=expand(P): R:=expand(R): Q:=expand(Q): Q1:=Q+diff(R,y): Q1:=expand(Q1): l1:=degree(R,y): l2:=degree(Q1,y): meka1:=coeff(R,y,l1): mekb1:=coeff(Q1,y,l2): if l1-1 <>l2 then k1:=degree(P,y)-max(l1-1,l2): else mumu:= -mekb1/meka1: if type(mumu,integer) and mumu > 0 then k1:=max(mumu, degree(P,y)-l2): else k1:=degree(P,y)-l2: fi: fi: fu:=0: if k1 < 0 then RETURN(0): fi: for ia1 from 0 to k1 do fu:=fu+apc[ia1]*y^ia1: od: gugu:=Q1*fu+R*diff(fu,y)-P: gugu:=expand(gugu): degg:=degree(gugu,y): for ia1 from 0 to degg do eq[ia1+1]:=coeff(gugu,y,ia1)=0: od: for ia1 from 0 to k1 do va1[ia1+1]:=apc[ia1]: od: for ia1 from 0 to ORDER do va2[ia1+1]:=bpc[ia1]: od: eq:=convert(eq,set): va1:=convert(va1,set): va2:=convert(va2,set): va1:=va1 union va2 : va1:=solve(eq,va1): fu:=subs(va1,fu): gorem:=subs(va1,gorem): if fu=0 and gorem=0 then RETURN(0): fi: for ia1 from 0 to k1 do gorem:=subs(apc[ia1]=1,gorem): od: for ia1 from 0 to ORDER do gorem:=subs(bpc[ia1]=1,gorem): od: fu:=normal(fu): shad:=pashetd(gorem,N): S:=lu1a*R*fu*shad[2]/P: S:=subs(va1,S): S:=normal(S): S:=factor(S): for ia1 from 0 to k1 do S:=subs(apc[ia1]=1,S): od: for ia1 from 0 to ORDER do S:=subs(bpc[ia1]=1,S): od: RETURN(shad[1],S): end: #AZd(A,y,n,N): inputs a hyperexponential expression A of the #contintuous variable y and the discrete variable n and #outputs a linear recurrence operator annihilating the #integral of A(y,n) along an interval where #A vanishes at the endpoints (including infinite intervals) #or a closed circle around the orign. For example try #AZd((1+y+y^2)^n/y^(n+1),y,n,N); AZd:= proc(A,y,n,N) local ORDER,gu,KAMA: KAMA:=200: for ORDER from 0 to KAMA do gu:=duisd(A,ORDER,y,n,N): if gu<>0 then RETURN(gu): fi: od: FAIL: end: #####End of stuff from AZ.txt #NRope(F,P,n,N,K): Inputs an expression F in variable P and #symbols n and N and a pos. integer K for #guessing, outputs a non-rogorously guessed #linear recurrence operator in the shift operator N # #annihilating the sequence of coeffs. of the unique f.p.s. #solution of #P(x)=x*F(P(x)).For example, #NRope(1+P+P^2,P,n,N,30); NRope:=proc(F,P,n,N,K) local L: L:=LIFtoS(F,P,K): GuessH(L,n,N): end: #Rope(F,P,n,N): Inputs an expression F in variable P and #symbols n and N and a pos. integer K for #guessing, outputs a rigorously derived (using the #not-as-famous-as-it-should-be Almkvist-Zeilberger algorithm #AND Lagrange-Inversion, a #linear recurrence operator in the shift operator N # #annihilating the sequence of coeffs. of the unique f.p.s. #solution of #P(x)=x*F(P(x)).For example, #Rope(1+P+P^2,P,n,N,30); Rope:=proc(F,P,n,N): #LIF: the n-th term is the coeff. of P^(n-1) in #F(P)^n divided by n AZd((F/P)^n/n,P,n,N)[1]: end: ####################### ######START HW 16###### ####################### #LargestAB(L): inputs a list of complex numbers, and finds the one with the largest absolute value. #(Use the evalf and evalc Maple commands) For example, #LargestAB([5,-3,-4,6,1]); #should give 6 and #LargestAB([5,-3,-4,6,-7,1]); #should give -7. LargestAB:=proc(L) local x,norm,max,maxNorm,numMax: if nops(L)=0 then RETURN(FAIL): fi: norm:=0: max:=0: maxNorm:=0: for x in L do: norm:=evalf(abs(evalc(x))): if norm > maxNorm then maxNorm:=norm: max:=x: fi: od: numMax:=0: for x in L do: norm:=evalf(abs(evalc(x))): if norm = maxNorm then numMax:=numMax+1: fi: od: if numMax>1 then RETURN(FAIL): fi: max: end: #FindMu(ope,n,N): inputs a linear recurrence operator with polynomial coefficients, ope(n,N), #and outputs the root of the largest absolute value of lcoeff(ope,n) (a certain polynomial in N), #if it is positive, or returns FAIL if it is negative, or if there are ties (i.e. both a and -a are roots), #or if the root(s) with maximal absolute value are complex. For example, #FindMu((n+2)*N^2-(3*n+3)*N+(2*n+5),n,N); #should return: 2. FindMu:=proc(ope,n,N) local mu, ope1,sol: ope1:=lcoeff(ope,n): sol:=[solve(ope1,N)]: mu:=LargestAB(sol): if mu=FAIL or evalf(mu)<0 then RETURN(FAIL): fi: mu: end: #FindTheta(ope,n,N,Ini): inputs a linear recurrence operator ope, phrased in terms of n and the shift N, #and Ini, a list of length degree(ope,N) indicating the initial values of the sequence, and outputs a #guessed RATIONAL number such that the sequence defined by #a(m):=HtoSeq(ope,n,N,Ini,m)[m] is given by: #a(m) IS ASYMPOTIC TO C*mu^m * m^theta #log(a(m+1)/(mu*a(m))) = theta*log((m+1)/m) #for some constant C, and the above mu that you found, using FindMu(ope,n,N). FindTheta:=proc(ope,n,N,Ini) local mu,a,thetaGuess,theta: mu:=FindMu(ope,n,N): if mu=FAIL then RETURN(mu): fi: a:=HtoSeq(ope,n,N,Ini,10000): thetaGuess:=evalf(log(a[10000]/(mu*a[9999]))/log(10000/9999)): convert(thetaGuess,confrac,'cvgts'): print(cvgts): theta:=cvgts[3]: theta: end: #FindC(ope,n,N,Ini): estimates the C in #a(m) IS ASYMPTOTIC TO C*mu^m* m^theta FindC:=proc(ope,n,N,Ini) local mu,theta,a: mu:=FindMu(ope,n,N): theta:=FindTheta(ope,n,N,Ini): a:=HtoSeq(ope,n,N,Ini,10000): identify(evalf(a[10000]/mu^10000/10000^theta)): end: #Asy(ope,n,N,Ini): inputs a linear recurrence operator ope #in n and the shift operator N, and the initial values #(of length the order of ope, i.e. degree(ope,N) and #outputs (parially non-rig.) the asymptic expression #for the sequence in the form #C*mu^n*n^theta Asy:=proc(ope,n,N,Ini) local mu,theta,C: mu:=FindMu(ope,n,N): theta:=FindTheta(ope,n,N,Ini): C:=FindC(ope,n,N,Ini): C*mu^n*n^theta: end: #AsyBS(F,n,k,L): inputs a binomial coefficient expression phrased in terms of n and k, #and a pos. integer L (for guessing, made large enough) and conjectures asymptotics for #a(n):=add(F(n,k),k=0..n). AsyBS:=proc(F,n,k,L) local ope,ORD,Ini,m,j: ope:=EmpiricalZeilberger(F,n,k,N,L): ORD:=degree(ope,N): Ini:=[seq(evalf(add(subs({n=m,k=j},F),j=0..m)) ,m=1..ORD)]: Asy(ope,n,N,Ini): end: #AsyBS(binomial(n,k)^2,n,k,50); #.56418253122204200526*4^n/sqrt(n) #AsyBS(binomial(n,k)^3,n,k,50); #.36754034533076489914*8^n/n #AsyBS(binomial(n,k)^4,n,k,50); #.25396025799131365077*16^n/n^(3/2) #AsyBS(binomial(n+k,k)^2*binomial(n,k)^2,n,k,50); #.22003455751284834012*(17+12*sqrt(2))^n/n^(3/2)