#Yusra Naqvi: HW #16 #OK to post ################################################################################ #1 #LargestAB(L): inputs a list of complex numbers and finds the one with the #largest absolute value LargestAB:=proc(L) local i,max: max:=L[1]: for i from 2 to nops(L) do if evalf(abs(L[i]))>evalf(abs(max)) then max:=L[i]: fi: od: max: end: ################################################################################ #2 #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 FindMu:=proc(ope,n,N) local L,mu: L:=[solve(lcoeff(ope,n)=0,N)]: mu:=LargestAB(L): if evalf(Im(evalc(mu)))<>0 or evalf(Re(evalc(mu)))<0 or -mu in L then return(FAIL): else return(mu): fi: end: ################################################################################ #3 #FindTheta(ope,n,N,Ini): inputs a linear recurrence operator ope(n,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 FindTheta:=proc(ope,n,N,Ini) local mu,a,cvgnts: mu:=FindMu(ope,n,N): if mu=FAIL then return(mu): fi: a:=HtoSeq(ope,n,N,Ini,1001): if a=FAIL then return(a): fi: #print(evalf(log(a[1001]/a[1000]/mu)/log(1001/1000))): convert(log(a[1001]/a[1000]/mu)/log(1001/1000),confrac,cvgnts): cvgnts[1]: 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, 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 return(FAIL): fi: L:=[op(L),-add(subs(n=n1,coeff(ope1,N,-i))*L[n1-i],i=1..ord)/k]: od: L: end: ################################################################################ #4 #FindC(ope,n,N,Ini): finds C such that a(m) is asymptotic to C*mu^m*m^theta FindC:=proc(ope,n,N,Ini) local a,mu,theta: option remember: a:=HtoSeq(ope,n,N,Ini,2001)[2000]: mu:=FindMu(ope,n,N): theta:=FindTheta(ope,n,N,Ini): identify(evalf(a/mu^2000/2000^theta)): end: ################################################################################ #5 #Asy(ope,n,N,Ini): finds the expression C*mu^m*m^theta, such that #a(m) is asymptotic to C*mu^m*m^theta Asy:=proc(ope,n,N,Ini) local m,C,mu,theta: option remember: mu:=FindMu(ope,n,N): theta:=FindTheta(ope,n,N,Ini): C:=FindC(ope,n,N,Ini): print(C*mu^m*m^theta): end: ################################################################################ #6 #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 N,ope,Ini,i: ope:=EmpiricalZeilberger(F,n,k,N,L): 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: Asy(ope,n,N,Ini): end: ##### #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: ##### We find: #AsyBS(binomial(n,k)^2,n,k,50); #25.22974831*4^m/m #AsyBS(binomial(n,k)^3,n,k,50); #.3674913416*8^m/m #AsyBS(binomial(n,k)^4,n,k,50); #11.35489279*16^m/m^2 #AsyBS(binomial(n+k,k)^2*binomial(n,k)^2,n,k,50); #9.838602178*(17+12*sqrt(2))^m/m^2 ################################################################################