#Feb. 18, 2016, trying to get famous like Roger Apery Help:=proc(): print(`Z(F,k,n), PRS(R,n,N0) , DS(F,k,n,N0), CS(F,k,n,N0), RAseq(F,k,n,N0)`): print(` EmpDel(F,k,n,N0,w) `): end: #Z(F,k,n): inputs an expression inolving binomial coefficients in k and n #outputs a P-recursive description of the sequence #Sum(F(k,n),k=0..n); Using the Zeilberger algorithm. Z:=proc(F,k,n) local N,ope, Ord,rec,i,n0,k0,ini: ope:=SumTools[Hypergeometric][Zeilberger](F,n,k,N)[1]: Ord:=degree(ope,N): ope:=expand(subs(n=n-Ord,ope)/N^Ord): ope:=add(factor(coeff(ope,N,-i))*N^(-i),i=0..Ord): ope:=ope/coeff(ope,N,0): rec:=[-seq(coeff(ope,N,-i),i=1..Ord)]: ini:=[seq( eval(add( subs({n=n0,k=k0},F),k0=0..n0)) ,n0=0..Ord-1)]: [ini,rec]: end: #PRS(R,n,N): inputs a P-recursive sequence [ini,rec(n)] where ini is a list #of numbers of length k say, starting at the 0-th term #and rec is a list of coefficients of the #P-recursive sequence defined UNIQELY by #x(n)=ini[n], for n=1..k #and rec is a list of numbers [c1(n),c2(n),...,ck(n)] #and for n>=k defined recursively by #x(n)=c1(n)*x(n-1)+c2(n)*x(n-2)+...+ck(n)*x(n-k). #For example the sequence 2^n is coded as #[[1],[2]], and the sequence of Fib. numbers are given by #[[0,1],[1,1]] #It outputs the list consisting of the first N terms of that sequence starting at n=0 PRS:=proc(R,n,N) local ini,rec,Se,k,n0,nt,i: ini:=R[1]: rec:=R[2]: k:=nops(ini): if nops(rec)<>k then RETURN(FAIL): fi: Se:=ini: for n0 from k to N do #nt:=c1(n0)*Se[n0-1]+c2(n0)*Se[n0-2]+...+ ck(n0)*Se[n0-k] nt:=add(subs(n=n0,rec[i])*Se[n0-i+1],i=1..k): Se:=[op(Se),nt]: od: Se: end: #DS(F,k,n,N0): The first N0+1 terms of the binomial coefficients sum #add(F(k,n),k=0..n); computed directly DS:=proc(F,k,n,N0) local k0,n0: [seq( eval(add( subs({n=n0,k=k0},F),k0=0..n0)) ,n0=0..N0)] end: #CS(F,k,n,N0): The first N0+1 terms of the binomial coefficients sum #add(F(k,n),k=0..n); computed cleverly using the amazing Z algorithm CS:=proc(F,k,n,N0) local R: R:=Z(F,k,n): PRS(R,n,N0): end: #RAseq(F,k,n,N0): inputs a binomial coeffient expression in the #symbols k, and n, and outputs the Apery-like sequence (of length N0+1) #of quotients of the a_n/b_n (see van der Poorten's famous paper, "A proof that Euler missed") #(1979, available (freel!) from google scholar) RAseq:=proc(F,k,n,N0) local R,bn, an,ord1,i: R:=Z(F,k,n): ord1:=nops(R[1]): if ord1<2 then RETURN(FAIL): fi: bn:=PRS(R,n,N0): an:=PRS([[0$(ord1-1),1], R[2]],n,N0): [seq(an[i]/bn[i],i=1..nops(an))]: end: #EmpDel(F,k,n,N0,w): the empirical delta of the Apery-like sequence #of rational numbers converging to their limit for terms between N0-w to N0 EmpDel:=proc(F,k,n,N0,w) local c,L: if (N0<40 or w>N0/2) then RETURN(FAIL): fi: L:=RAseq(F,k,n,4*N0): c:=L[nops(L)]: evalf([seq(-log(abs(L[i]-c))/log(denom(L[i])),i=N0-w..N0)],10): end: