######################################################################
##AperyAcc: Save this file as AperyAcc. To use it, stay in the       #   
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read AperyAcc : <Enter>                            #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg@math.rutgers.edu.                                          # 
######################################################################

Digits:=100:

#VERSION of July 8, 2020
#Created: Jan. 2, 2002
#This version: Jan. 2, 2001
#AperyAcc: A Maple package to study Apery's Original Method
#As described in:"Interpolation de Fractions Continues et
# Irrationalite de Certaines Constantes":
#Bull. de la section des sci. du C.T.H.S. no. 3, p. 37-53
#Bib. Nat. Paris, Math Reviews 83b:10039
#Please report bugs to zeilberg@math.rutgers.edu




print(` AperyAcc: A Maple package to study Apery's Original Method`):
print(`As described in his masterpiece:`):
print(`"Interpolation de Fractions Continues et Irrationalite de `):
print(` Certaines Constantes":`):
print(`Bull. de la section des sci. du C.T.H.S. no. 3, p. 37-53`):
print(` (Math Reviews 83b:10039) `):
print(``):
print(`It is one of the Maple packages accompanying the article`):
print(`COMPUTERIZED DECONSTRUCTION by Doron Zeilberger`):
print(``):
print(` Please report bugs to zeilberg@math.rutgers.edu `):
print(`Written by Doron Zeilberger, zeilberg@math.rutgers.edu`):

print(`Created: Jan. 2, 2002.`):
print(`This version: Jan. 2, 2002 `):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name); .` ):
print(`For a list of all procedures, type: ezra1();  .`):
 print(``):
ezra1:=proc()
if args=NULL then
print(` AperyAcc: A Maple package to study Apery's Original Method`):
 print(`Contains the following procedures: AccRec, AccRecN, Aitken `):
 print(` Aluf, Apery, Aperyh,Aperyh1, AperyhOper`):
 print(`AperyNes, AperySeq , Appx`):
 print(` ApplyDiffOpe, ApplyIni, ApplyOpe, ApplyRec, Appx, `):
  print(` BestAcc, BestAcc1, CheckC, CheckN `):
 print(` CheckOpehz2, CODV, CODVpol, DiagPQ, DiagLogt, Diagz2, Diagz3,`):
 print(` delt,delta, Expl, FindaB, FindAcc, FindB, FindB1, FindOpe`):
 print(` FindOpeNew, FindOpePol, GF1Q, GFLogtqh,`):
 print(`GFz2qh, GFz3qh, GuessPol, Halevay, hthRow , ilcmP, ImParSum, `):
 print(`IniOpe, IsNumericSol, IsRatPol, Lucky, MainDiagPQ, `):
 print(` MakePol, Mult, OpeDh, OpehLogt `):
 print(` Opehz2, Opehz3, ParSum,`):
 print(` pBfromC, PhLogt, Phz2, Phz3, pnh, PolNiceBasis`):
 print(` PQtoOpe, qnh, Qnh0, RatImp, RatImpD, RatImpCfD, RatImpD1 `):
 print(` RecToDiff, ResultOpe, ResultOpe1, Rnh0, `):
 print(`SeqFromRec, SeqQuotFromRec, ShiftPQrs, Shiftr, Sidra, SidraCf `):
 print(` yafe , Ztrans`):
 print(`For help with any one of them type: ezra1(ProcedureName);  `):
 else
 ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then
 print(`For a list of all the procedures type ezra1();`):
 print(`The main procedures are : `):
 print(` AccRec, Apery, Aperyh,  AperyNes, AperySeq, Appx, AppxDel, RatImp`):
fi:


if nops([args])=1 and op(1,[args])=AccRec then
print(`AccRec(pol,x,n,N): Given a polynomial pol`):
print(`in the variable n, it tries to uses Apery's method to `):
print(`find a convergence-acceleration recurrence (using the diagonal)`):
print(`followed by the initial conditions [[a0,a1],[b0,b1]]`):
print(`such that the constant sum(x^n/pol(n),n=1..infinity) is`):
print(`approximated by the quotients of solution of the recurrence`):
print(`with the given initial conditions for the numerator and`):
print(`denominator. For example, for a fast recurrence for Catalan's:`):
print(`Constant, found recently by Wadim Zudilin in a semi-human way`):
print(`type: AccRec(-(2*n-1)^2,-1,n,N);`):
fi:

if nops([args])=1 and op(1,[args])=AccRecN then
print(`AccRecN(pol,n,x0,N,g): Like AccRec(pol,x0,n,N), but with`):
print(` a normalizing function g entered `):
print(`It returns: what AccRec returns plus the normalized diagonal`):
print(` sequences (the pnh(n,n) multiplied by ilcm(pol(1),..,pol(n)) `):
print(` For example, try: `):
print(` AccRecN(n^2,1,n,N,2^n/n!^2); AccRecN(n^3,1,n,N,1/n!^3); `):
fi:

if nops([args])=1 and op(1,[args])=Aitken then
print(` Aitken(List): The Aitken transform of the sequence List `):
fi:

if nops([args])=1 and op(1,[args])=Aluf then
print(`Aluf(pol,Rat,n,L): the best appx. (according to delta)`):
print(`amongst sum(1/pol(i),i=1..n0)+Rat(n0) to `):
print(` a:=sum(1/pol(i),i=1..infinity) for n0=2..L`):
print(` For example, try: Aluf(n^2,1/(n+1/2),n,5); `):
fi:

if nops([args])=1 and op(1,[args])=Apery then
print(`Apery(pol,x,n,R):Given a polynomial pol in the variable n`):
print(`with polynomial coefficient, uses Apery's ORIGINAL method`):
print(`for TRYING to prove irrationality for the constant`):
print(`sum(x^i/pol(i),i=1..infinity)`):
print(`as described in his beautiful C.T.H.S. paper. It continues`):
print(`to accelerate up to R rows, if possible, and returns`):
print(`the list of accelerating B's, followed by the list of P's`):
print(`Followed by the list of Q's`):
print(`that feature in the recurrence for row h`):
print(`followed by the explicit expression for row h, h=0,1,2,..`):
print(` that is "easy" solution in row i divided by pol(1)pol(2)...pol(n) .`):
print(` If it fails, then it returns 0 `):
print(`For example, try Apery(n^2,1,n,5);, Apery(n^3,1,n,5);`):
fi:

if nops([args])=1 and op(1,[args])=Aperyh then
print(`Aperyh(pol,x,n,h):Given a polynomial pol in the variable n`):
print(`uses Apery's ORIGINAL method`):
print(`To guess a convergence-acceleration scheme for`):
print(`the constant sum(x^n/pol(n),n=1..infinity)`):
print(`Here h is a variable.`):
print(`It returns B(n,h),P(n,h),Q(n.h) if successful`):
print(` otherwise it returns 0. For example, try `):
print(` Aperyh(n^2,1,n,h); and Aperyh(n^3,1,n,h); `):
fi:

if nops([args])=1 and op(1,[args])=Aperyh1 then
print(`Aperyh1(pol,x,n,R,h):Given a polynomial pol in the variable n`):
print(`with polynomial coefficient, uses Apery's ORIGINAL method`):
print(`and a variable h tries, if possible, to fit the`):
print(`accelerating B(n)'s and the coeffs.(P,Q) of the recurrence`):
print(`for row h by polynomials of degree h`):
print(`by fitting the data from Apery(pol,n,R).`):
print(`It returns B(n,h),P(n,h),Q(n.h)if successful`):
print(` otherwise it returns 0. For example, try `):
print(` Aperyh1(n^2,1,n,8,h); and Aperyh1(n^3,1,n,8,h); `):
fi:

if nops([args])=1 and op(1,[args])=AperyNes then
print(`AperyNes(pol,x,n,N,h) : Given a polynomial pol, in n`):
print(`finds, if possible, an Apery-style miracle pair`):
print(`accelarating operator: Acc0=N+B(n+1,h+1) and a row-operator`):
print(`ope:=N^2+P(n,h)*N+pol(n+1)^2 such that`):
print(`ResultOpe(ope,[1,0],[0,1],Acc,n,N) is ope with`):
print(`h replaced by h+1, and B is a poly of (n,h) of total degree`):
print(`degree(pol,n), and P is a (not necessarily homog.) `):
print(`polynomial of degree degree(pol,n)`):
print(`For example, try:AperyNes(n^2,1,n,N,h) `):
fi:

if nops([args])=1 and op(1,[args])=AperyhOper then
print(`AperyhOper(pol,x,n,h,N):Like Aperyh(pol,x,n,h,N) but`):
print(`returns the accelerating operator followed by `):
print(`the recurrence satistied in the n-direction`):
print(`For example, try: AperyhOper(n^2,1,n,h,N);`):
fi:

if nops([args])=1 and op(1,[args])=AperySeq then
print(`AperySeq(pol,x,n,L): `):
print(`Returns the sequence of diagonal approximants`):
print(`up to the L-th term  `):
print(`in Apery's scheme for sum(x^n/pol(n),n=1..infinity)`):
print(`if it fails, it returns 0`):
print(` For example, try `):
print(` AperySeq(n^2,1,n,40); and AperySeq(n^3,1,n,40); `):

fi:

if nops([args])=1 and op(1,[args])=ApplyDiffOpe then
print(`ApplyDiffOpe(ope,x,D1,F): applies the differential operator`):
print(`ope(x,D1) to F(x)`):
fi:


if nops([args])=1 and op(1,[args])=ApplyIni then
print(`ApplyIni(ope,Q,n,N,Lis1): Given a recurrence operator `):
print(`ope(n,N) and another one Q(n,N)`):
print(`in the variable n and its corresponding shift N`):
print(`and a list Lis1 of (size degree(Q,N)) of initial values`):
print(`finds the initial values for the sequence ope(n,N) applied to`):
print(`the sequence whose inital values are given by Lis1 `):
print(`and annihilated by Q(n,N) For example, try:`):
print(`ApplyIni(N+2,N-1,n,N,[1]);`):
fi:

if nops([args])=1 and op(1,[args])=ApplyOpe then
print(`ApplyOpe(P,Q,n,N): Given two recurrence operators P(n,N) and Q(n,N)`):
print(`finds the operator R(n,N) of order degree(P,N)-1 such that`):
print(`if f is a solution of P(n,N)f=0 then Q(n,N)f=R(n,N)f`):
print(`For example, try ApplyOpe(N^2-N-1,N^2-N-1,n,N)`):
fi:

if nops([args])=1 and op(1,[args])=ApplyRec then
print(`ApplyRec(ope,n,N,Lis): applies the recurrence ope(N,n) to the `):
print(`list Lis, For example, try: ApplyRec(N^2-N-1,n,N,[1,1,2,3,5,8]);`):
fi:

if nops([args])=1 and op(1,[args])=Appx then
print(`Appx(pol,x,n,L): Gives the Lth diag. approximant`):
print(`Apery's recurrence for x/pol(1)+x^2/pol(2)+x^3/pol(3)+... `):
print(`followed by the error. For example: try Appx(n^2,1,n,20);`):
print(`to get a good appx. for Pi^2/6, or to get an approximation for`):
print(` Catalan's constant, type: Appx(-(2*n-1)^2,-1,n,20);`):
fi:

if nops([args])=1 and op(1,[args])=AppxDel then
print(`AppxDel(pol,x,n,L): Gives the Lth diag. approximant`):
print(`Apery's recurrence for x/pol(1)+x^2/pol(2)+x^3/pol(3)+... `):
print(`followed by the error and the delta. For example: `):
print(`try AppxDel(n^2,1,n,20);`):
fi:

if nops([args])=1 and op(1,[args])=BestAcc then
print(`BestAcc(ope,IniT,IniB,n,N,ListDegs): Given an operator ope(n,N).`):
print(`and initial values IniT, IniB (for top and bottom)`):
print(`and a list of degrees ListDegs tries to find`):
print(`the best accelerating operator ope1=P0(n)+P1(n)*N+P2(n)*N^2+...`):
print(`with ListDegs=[degree(P0(n),n),degree(P1(n),n),..]`):
print(`For example, try: BestAcc(IniOpe(n^2,1,n,N),n,N,[2,0]);`):
fi:

if nops([args])=1 and op(1,[args])=BestAcc1 then
print(`BestAcc1(ope,IniT,IniB,n,N,ListDegs,L): Given an operator ope(n,N).`):
print(`and initial values IniT, IniB (for top and bottom)`):
print(`and a list of degrees ListDegs tries to find`):
print(`the best accelerating operator (within L)`):
print(`ope1=P0(n)+P1(n)*N+P2(n)*N^2+...`):
print(`with ListDegs=[degree(P0(n),n),degree(P1(n),n),..]`):
print(`For example, try: BestAcc1(IniOpe(n^2,1,n,N),n,N,[2,0],0);`):
fi:

if nops([args])=1 and op(1,[args])=CheckC then
print(`CheckC(C,n,n0): Given an expression C(n) and an integers n0`):
print(` computes `):
print(` sum((C(i+1)-C(i))/(C(i+1)*C(i)+C(i+1)+C(i)),i=n0+1..infinity)- `):
print( ` 1/(C(n0+1)). It should be small `):
print(` For example, try CheckC(2*n*(n-1),n,10); `):
fi:

if nops([args])=1 and op(1,[args])=CheckN then
print(`CheckN(pol,x0,n,h,g,L): Given a polynomial pol, in the variable`):
print(` n, a number x0,`):
print(` and another variable h, and a normalization function g `):
print(` checks whether q(n,h)*g(n,h) are integers and `):
print(` lcm(pol(1), ..., pol(n))*p(n,h)*g(n,h) are integers `):
print(` For example, try: CheckN(n^3,1,n,h,1/h!^3,20); `):
fi:

if nops([args])=1 and op(1,[args])=CheckOpehz2 then
print(`CheckOpehz2(h): Checks The recurrence operator annihilated by`):
print(`the h-th row for Apery's scheme for zeta(2)`):
print(`For example, try: CheckOpehz2(5)`):
fi:

if nops([args])=1 and op(1,[args])=CODV then
print(`CODV(ope,n,m,N,M,g): Given a partial recurrence operator`):
print(`(with poly. coeffs.), ope in n,m, and the shift operators`):
print(`N and M, and a closed-form expression g (in n and m)`):
print(`outputs the operator annihilating g(n,m)F(n,m)`):
print(`where F(n,m) is a function annihilated by ope`):
print(`For example, try CODV(N-1,n,m,N,M,2^n)`):
fi:

if nops([args])=1 and op(1,[args])=CODVpol then
print(`CODVpol(ope,n,N,pol): Given a partial recurrence operator`):
print(`(with poly. coeffs.), ope in n and the shift operator`):
print(`N, and a polynomial pol(in n) `):
print(`outputs the operator annihilating F(n)/(pol(1)...pol(n))`):
print(`where F(n) is a function annihilated by ope`):
print(`For example, try CODVpol(N-2*(n+1),n,N,n)`):
fi:

if nops([args])=1 and op(1,[args])=DiagPQ then
print(`DiagPQ(P,Q,n,m,N,M): Input: recurrence operators P(n,m,N)`):
print(`and Q(n,m,N) (n,m, are variables and N and M are the shift`):
print(`operators in n,m, resp.); Output: an operator of the form`):
print(`R(n,m,N*M) that annihilates any function that is annihilated `):
print(`by P(n,m,N) and M-Q(n,m,N) For example, try:`):
print(`DiagPQ(Opehz2(n,N,h) ,Phz2(n,N,h),n,h,N,H);`):
fi:

if nops([args])=1 and op(1,[args])=DiagLogt then
print(`DiagLogt(N,n,t): the recurrence satisfied by the diagonal`):
print(`of Apery's accelerating scheme for  Log(t+1)` ):
fi:

if nops([args])=1 and op(1,[args])=Diagz2 then
print(`Diagz2(N,n): the recurrence satisfied by the diagonal`):
print(`of Apery's accelerating scheme for  Zeta(2)`):
fi:

if nops([args])=1 and op(1,[args])=Diagz3 then
print(`Diagz3(N,n): the recurrence satisfied by the diagonal`):
print(`of Apery's accelerating scheme for  Zeta(3)`):
fi:

if nops([args])=1 and op(1,[args])=delt then
print(`delt(a,c): Given a rational number a and a constant`):
print(`c, finds the delta such that |a-c|=1/denom(a)^(1+delta)`):
print(`For example, try delta(22/7,evalf(Pi)):`):
fi:
 
if nops([args])=1 and op(1,[args])=Expl then
print(`Expl(P,N): Given a recurrence operator P(n,N) (or order I)`):
print(` outputs it in explicit form with f(n+I) expressed in terms`):
print(`of (n),...,f(n+I-1) For example, try Expl(N^2-N-1,N)`):
fi:

if nops([args])=1 and op(1,[args])=FindaB then
print(`FindaB(Q1,Q,pol,n): Given polynomials Q1(n),Q(n),`):
print(`pol(n) tries to find a constant a and a polynomial`):
print(`B(n) such that a*Q1(n)-Q(n+1)*pol(n+1)-Q(n)*B(n)=0`):
print(`If unsuccesful, returns 0`):
print(`For example, try: FindaB(2*n^3+3*n^2+3*n+1,n,n^2,n);`):
fi:

if nops([args])=1 and op(1,[args])=FindAcc then
print(`FindAcc(pol,n): Uses Apery's method to find a rational`):
print(`function R(n) such that sum(1/pol(i),i=1..n)+R(n)`):
print(`is a better appx. to sum(1/pol(i),i=1..infinity) then`):
print(`the the partial sum itself`):
print(`For example, try FindAcc(n^4+7/16,n);`):
fi:

if nops([args])=1 and op(1,[args])=FindB then
print(`FindB(P,Q,n): Uses Apery's method to find a polynomial`):
print(`B such that B(n)*(Q(n)+B(n+1))-P(n) is of  smallest degree`):
print(`possible. For example: try FindB(-n^4,2*n^2+2*n+1,n);`):
print(`and FindB1(-n^4,2*n^2+2*n+3,n);`):
print(`It exists {} if it can't find any of degree lower than P`):
fi:

if nops([args])=1 and op(1,[args])=FindB1 then
print(`FindB1(P,Q,n,L): Uses Apery's method to find a polynomial`):
print(`B such that B(n)*(Q(n)+B(n+1))-P(n) is of  degree L`):
print(`possible. For example: try FindB1(-n^4,2*n^2+2*n+1,n,0);`):
print(`and FindB1(-n^4,2*n^2+2*n+3,n,0);`):
fi:


if nops([args])=1 and op(1,[args])=FindOpe then
print(`FindOpe(Q,d,n,N): tries to find a second-order operator of the`):
print(`form ope(N,n)=P2*N^2+P1(n)*N+P0(n) where P0,P1,P2 are`):
print(`polynomials of degree d`):
print(`annihilating the polynomial Q(n)`):
print(`For example, try FindOpe(n,0,n,N);`):
fi:

if nops([args])=1 and op(1,[args])=FindOpeNew then
print(`FindOpeNew(Q,d,n,N,pol): tries to find a second-order`):
print(` operator of the`):
print(`form ope(N,n)=pol(n+2)*P2*N^2+P1(n)*N+pol(n+1)*P0(n) `):
print(`where ope(N,n) is`):
print(`of degree d in n,`):
print(`annihilating the polynomial Q(n)`):
fi:

if nops([args])=1 and op(1,[args])=FindOpePol then
print(`FindOpePol(Q,pol,n,N): tries to find an operator of the`):
print(`form ope(N,n)=pol(n+1)*pol(n+2)*N^2+pol(n+1)*P1(n)*N+P2(n)`):
print(`annihilating the polynomial Q(n)`):
print(`For example, try FindOpePol(1,n,n,N);`):
fi:

if nops([args])=1 and op(1,[args])=ImParSum then
print(`ImParSum(pol,n,n0): The improved partial sum`):
print(`e.g. ImParSum(n^3,n,10)`):
fi:



if nops([args])=1 and op(1,[args])=GF1Q then
print(`GF1Q(pol,x0,n,x,L): the list of the generating functions`):
print(`sum(q(n,h)*x^n,n=1..infinity), h=0..L`):
print(`Where q(n,h) are the denominators in Apery's scheme`):
print(`for approximating sum(x0^n/pol(n),n=1..infinity)`):
print(` for example, try`):
print(`GF1Q(n^2,1,n,x,4): `):
fi:

if nops([args])=1 and op(1,[args])=GF1R then
print(`GF1R(pol,x0,n,x,L): the list of the generating function`):
print(`of the first L+1 rows in Apery's scheme for r(n,h) where`):
print(` for example, try`):
print(`GF1R(n^2,1,n,x,4): `):
fi:

if nops([args])=1 and op(1,[args])=GFLogtqh then
print(`GFLogtqh(x,h,t): The generating function`):
print(`(1/h!)sum(qnh(n,h)/n!*x^n,n=0..infinity)`):
print(` implied by PhLogt `):
print(`For example, try: GFLogtqh(x,3,t)`):
fi:

if nops([args])=1 and op(1,[args])=GFz2qh then
print(`GFz2qh(x,h): The generating function`):
print(`(1/h!^2)sum(qnh(n,h)/n!^2*x^n,n=0..infinity)`):
print(` implied by Phz2  `):
print(`For example, try: GFz2qh(x,3)`):
fi:

if nops([args])=1 and op(1,[args])=GFz3qh then
print(`GFz3qh(x,h): The generating function`):
print(`(1/h!^3)sum(qnh(n,h)/n!^3*x^n,n=0..infinity)`):
print(` implied by Phz3 `):
print(`For example, try: GFz3qh(x,3)`):
fi:

if nops([args])=1 and op(1,[args])=GuessPol then
print(`GuessPol(List,h): Given a list of expressions, List,`):
print(`guesses a polynomial P,of degree<=nops(List)-3 such that`):
print(`it fits List=[P(0),P(1),P(nops(List)-1)], and returns`):
print(`0 if none exists`):
fi:


if nops([args])=1 and op(1,[args])=Halevay then
print(`Halevay(pol,n,h,N,start1,end1,res1,degn,Maxdegh): Given a polynomial`):
print(`pol, tries to guess second-order operators of degn`):
print(`in n and <=Maxdegh in h for the rows of the output of RatImpD`):
print(`starting at start1, ending at end1, and with resolution res1`):
print(`For example, try:`):
print(`Halevay(n^2,n,h,N,2,6,1,2,2)`):
fi:

if nops([args])=1 and op(1,[args])=hthRow then
print(`hthRow(n,pol,h0,L): The first L+1 terms in the h0^th row`):
print(`of the Apery Scheme for 1/pol(1)+...+1/pol(n)+...`):
fi:

if nops([args])=1 and op(1,[args])=ilcmP then
print(`ilcmP(pol,n,n0): The least common multiple of pol(1), ..., pol(n0)`):
print(`where pol is a polynomial in n `):
print(`For example, do: ilcmP(n^2,n,5); `);
fi:


if nops([args])=1 and op(1,[args])=IniOpe then
print(`IniOpe(pol,x,n,N): The second-order operator,`):
print(`in the variable n and its shift operator N,  satisfied`):
print(`by the numerator and denominator of sum(x^i/pol(i),i=1..n)`):
print(`followed by the initial values for the numerator and deno.`):
print(` respectively `):
print(`For example, try IniOpe(n^2,1,n,N);`):
fi:

if nops([args])=1 and op(1,[args])=IsNumericSol then
print(`IsNumericSol(Sol): Given a solution set, Sol,`):
print(`returns true iff it is numeric;`):
print(`For example, try: IsNumericSol({a=a}); IsNumericSol({a=5});`):
fi:

if nops([args])=1 and op(1,[args])=Lucky then
print(`Lucky(n,x,deg,a,b): Given a variable n, and letters a,b, an integer`):
print(`deg, finds the general monic polynomial pol of dergee `):
print(`deg that allows Apery-acceleration for sum(1/pol(i),i=1..infinity)`):
print(`For example, try: Lucky(n,1,2,a,b);`):
fi:

if nops([args])=1 and op(1,[args])=IsRatPol then
print(`IsRatPol(P,n): Given a polynomial P in n decides whether`):
print(`all its coefficients are rational numbers `):
fi:

if nops([args])=1 and op(1,[args])=MainDiagPQ then
print(`MainDiagPQ(P,Q,n,m,N,L1,L2): Like DiagPQ, but `):
print(`there is no input of the shift M, and instead`):
print(`there is a normalizing sequence L1(n) and L2(m)`):
print(`and the output is for the recurrence`):
print(` along the main diagonal using N `):
fi:

if nops([args])=1 and op(1,[args])=MakePol then
print(`MakePol(sid,x,y,deg): Given a sequence sid, of polynomials`):
print(`in x,y tries to find a sequence a[h], of the same`):
print(`length of sid, such that sid[h]*a[h] is a polynomial of degree`):
print(`deg in h`):
print(`For example, try `):
print(`MakePol([(1+n*m+n^2)/2,(1+n*m+n^2)*3,(1+n*m+n^2)*4],n,m,1);`):
fi:

if nops([args])=1 and op(1,[args])=Mult then
print(`Mult(P,Q,N,n): Given two recurrence operators P(N,n)`):
print(`and Q(N,n), (where n is the discrete variable and N is the`):
print(`shift operator N : Nx(n):=x(n+1), computes the`):
print(`product P(N,n)Q(N,n); For example, try:`):
print(`Mult(N,n^2*N,N,n)`):
fi:

if nops([args])=1 and op(1,[args])=OpeDh then
print(`OpeDh(pol,n,x,D1,h0) The differential operator that yields`):
print(`the generating function for the (h0-1)th row to the h0^th row`):
print(`for Apery's scheme for 1/pol(1)+...+1/pol(n)+...`):
print(`For example, try: OpeDh(n^2,n,x,D1,2);`):
fi:

if nops([args])=1 and op(1,[args])=OpehLogt then
print(`OpehLogt(n,N,h,t): The recurrence operator annihilated by`):
print(`the h-th row for Apery's scheme for log(1+t)`):
fi:

if nops([args])=1 and op(1,[args])=Opehz2 then
print(`Opehz2(n,N,h): The recurrence operator annihilated by`):
print(`the h-th row for Apery's scheme for zeta(2)`):
fi:

if nops([args])=1 and op(1,[args])=Opehz3 then
print(`Opehz3(n,N,h): The recurrence operator annihilated by`):
print(`the h-th row for Apery's scheme for zeta(3)`):
fi:

if nops([args])=1 and op(1,[args])=ParSum then
print(`ParSum(pol,n,n0): sum(1/pol(i),i=1..n0)`):
fi:

if nops([args])=1 and op(1,[args])=pBfromC then
print(`pBfromC(C,n): Gives a generic good pair (p,B) from C`):
print(`For example, try pBfromC(2*n*(n-1),n)`):
fi:

if nops([args])=1 and op(1,[args])=PhLogt then
print(`PhLogt(n,N,h,t): The accelerating operator of Apery for log(1+t) `):
print(`going from the h-th row to the next`):
fi:

if nops([args])=1 and op(1,[args])=Phz2 then
print(`Phz2(n,N,h): The accelerating operator of Apery for Zeta(2) `):
print(`going from the h-th row to the next`):
fi:

if nops([args])=1 and op(1,[args])=Phz3 then
print(`Phz3(n,N,h): The accelerating operator of Apery for Zeta(3) `):
print(`going from the h-th row to the next`):
fi:

if nops([args])=1 and op(1,[args])=pnh then
print(`pnh(n,h,pol,x,B,n0,h0) Given a polynomial in the variable n`):
print(`and an accelerating polynomial B(n,h), computes the numerator`):
print(`in the approximating scheme for the n0-th entry in row h0`):
print(`already normalized by dividing by pol(1)..pol(n)`):
print(`for example, try: pnh(n,h,n^2,1,-n^2+(h+1)*n-(h+1)^2/2,2,2);`):
fi:

if nops([args])=1 and op(1,[args])=PolNiceBasis then
print(`PolNiceBasis(pol,n): Given a polynomial pol in the variable`):
print(` n, expresses it in terms of binomial(n+i,d), (i=0..d) `):
print(` where d is the degree, e.g., try: PolNiceBasis(n^2+n,n); `):
fi:

if nops([args])=1 and op(1,[args])=PQtoOpe then
print(`PQtoOpe(P,Q,N,n): Given functions of n, P(n) and Q(n)`):
print(`constructs the operator: N^2-Q(n+1)*N-P(n+1)`):
fi:

if nops([args])=1 and op(1,[args])=qnh then
print(`qnh(n,h,pol,B,n0,h0) Given a polynomial in the variable n`):
print(`and an accelerating polynomial B(n,h), computes the denominator`):
print(`in the approximating scheme for the n0-th entry in row h0`):
print(`already normalized by dividing by pol(1)..pol(n)`):
print(`for example, try: qnh(n,h,n^2,-n^2+(h+1)*n-(h+1)^2/2,2,2);`):
fi:

if nops([args])=1 and op(1,[args])=Qnh0 then
print(`Qnh0(pol,B,n,h,h0): an exact expression, in terms of the symbol n`):
print(`for qnh in Apery's scheme for h=h0. For example, try:`):
print(` Qnh0(n^2,Aperyh(n^2,1,n,h)[1],n,h,3): `):

fi:
if nops([args])=1 and op(1,[args])=RatImp then
print(`RatImp(pol,n,x,h,h0): The rational functions that`):
print(`improves the convergence of the partial sum of`):
print(`sum(x^i/pol(i),i=1..n) to the infinite sum`):
print(`implied by the h0-th row of Apery's scheme`):
print(` for example, try RatImp(n^2,n,1,h,3); `):
fi:

if nops([args])=1 and op(1,[args])=RatImpCfD then
print(`RatImpCfD(pol,Cf,n,deg): Given a rational function pol(n)`):
print(`and a closed-form function pol(n)`):
print(`Tries to find the rational function`):
print(`R(n), with monic denominator of of degree deg,`):
print(`(Cf(n)*R(n)-Cf(n-1)*R(n-1)-Cf(n)/pol(n))/Cf(n) `):
print(`has numerator of degree`):
print(`as small as possible`):
print(`For example, try: RatImpCfD(n^2,1,n,3);`):
fi:

if nops([args])=1 and op(1,[args])=RatImpD then
print(`RatImpD(pol,n,deg): Tries to find the rational function`):
print(` R(n), with monic denominator of of degree deg, `):
print(` R(n)-R(n-1)+1/pol(n) has numerator of degree `):
print(` as small as possible `):
print(` For example, try: RatImpD(n^2,n,3); `):
fi:

if nops([args])=1 and op(1,[args])=RatImpD1 then
print(`RatImpD1(pol,n,deg,L): Tries to find the rational function`):
print(`R(n), with monic denominator of of degree deg,`):
print(` R(n)-R(n-1)+1/pol(n) has numerator of degree L `):
print(` For example, try: RatImpD1(n^2,n,3,0);`):

fi:

if nops([args])=1 and op(1,[args])=RecToDiff then
print(`RecToDiff(ope,N,n,x,D1): given a recurrence operator`):
print(`ope(N,n) finds the differential operator P(x,D1) (D1=d/dx)`):
print(`annihilating the generating function f(x) of the sequence that is `):
print(`annihilated by ope(N,n). For example, try RecToDiff(N-1,N,n,x,D1)`):
fi:



if nops([args])=1 and op(1,[args])=ResultOpe then
print(` ResultOpe(Ope1,IniTop,IniBot,Ope2,n,N): `):
print(` Given two recurrence operators Ope1(n,N) and `):
print(` Ope2(n,N), and lists IniTop, IniBot, for the`):
print(` initial values `):
print(` finds the operator Ope3(n,N) of the same order as Ope1 such that`):
print(`if f is a solution of Ope1(n,N)f=0 then Ope3(n,N)(Ope2(n,N)f)=0`):
print(`followed by the initial values for the top and bottom`):
print(`of the sequence Ope2(n,N)f`):
print(`For example, try ResultOpe(N^2-N-1,[0,1],[1,0],N+1,n,N) `):
fi:


if nops([args])=1 and op(1,[args])=ResultOpe1 then
print(`ResultOpe1(Ope1,Ope2,n,N): Given two recurrence operators`):
print(` Ope1(n,N) and Ope2(n,N)`):
print(`finds the operator Ope3(n,N) of the same order as Ope1 such that`):
print(`if f is a solution of Ope1(n,N)f=0 then Ope3(n,N)(Ope2(n,N)f)=0. `):
print(`And the coeff. of N^ORDER is a const., it returns 0 otherwise`):
print(`For example, try ResultOpe1(N^2-N-1,N+1,n,N) `):
fi:

if nops([args])=1 and op(1,[args])=Rnh0 then
print(`Rnh0(pol,B,n,h,h0): an exact expression, in terms of the symbol n`):
print(`for r(n,h) in Apery's scheme for h=h0.`):
print(`Where r(n,h) is the polynomial in h such that`):
print(`p(n,h)=(1/pol(1)+...+1/pol(n))*q(n,h)+r(n,h)`):
print(`For example, try: Rnh0(n^2,Aperyh(n^2,1,n,h)[1],n,h,3):`):
fi:

if nops([args])=1 and op(1,[args])=SeqFromRec then
print(`SeqFromRec(ope,n,N,ini,L): Given an (ordinary) recurrence operator`):
print(`ope(n,N) in the variable n and the shift operator N (Nf(n):=f(n+1))`):
print(`and given initial values [ini0,ini1,..,ini_{ORD-1}],computes`):
print(`the first L+1 terms of the sequence a(n) satisfying`):
print(`ope(n,N)a(n)=0 and a[i]=ini[i] for i=0,...,ORD-1 .`):
print(`For example, try: SeqFromRec(N^2-N-1,n,N,[0,1],10);`):
fi:
if nops([args])=1 and op(1,[args])=SeqQuotFromRec then
print(`SeqQuotFromRec(ope,ini1,ini2,n,N,L): `):
print(`Given an (ordinary) recurrence operator`):
print(`ope(n,N) in the variable n and the shift operator N (Nf(n):=f(n+1))`):
print(`and given initial values ini1,ini2 finds`):
print(`the first L+1 terms of the sequence a(n)/b(n) where both `):
print(`a(n) and b(n) are annihilated by ope(n,N)`):
print(`and a[i]=ini1[i] for i=0,...,ORD-1 ; `):
print(`b[i]=ini2[i] for i=0,...,ORD-1`):
print(`For example, try: SeqQuotFromRec(N^2-N-1,[0,1],[1,1],n,N,10);`):
fi:


if nops([args])=1 and op(1,[args])=ShiftPQrs then
print(`ShiftPQrs(P,Q,n,m,N,r,s): If f is a solution of P(n,m,N)f=0, `):
print(`and Mf(m,n)=Q(n,m,N)f finds the operator R(n,m,N)`):
print(`of order degree(P,N)-1 such that`):
print(`N^rM^sf=R(n,m,N)f. For example, try: `):
print(`ShiftPQrs(N^2-N-1,1+n*N,n,m,N,3,2);`):
fi:

if nops([args])=1 and op(1,[args])=Shiftr then
print(`Shiftr(P,n,N,r): If f is a solution of P(n,N)f=0, then`):
print(`the operator Q(n,N) or order degree(P,N)-1 such that`):
print(` N^rf=Q(n,N)f  `):
print(` For example, try: Shiftr(N^2-N-1,n,N,3);`):
fi:

if nops([args])=1 and op(1,[args])=Sidra then
print(`Sidra(pol,n,L): The sequence of best appx. followed by the`):
print(`delta to sum(1/pol(i),i=1..infinity) furnished by the`):
print(`rational-function corrections supplied by `):
print(`RatImpD(pol,n,d2), d2=degree(pol,n)-1..L`):
print(`For example, try: Sidra(n^2,n,10);`):
fi:

if nops([args])=1 and op(1,[args])=SidraCf then
print(`SidraCf(pol,Cf,n,L): The sequence of best appx. followed by the`):
print(`delta to sum(Cf(i)/pol(i),i=1..infinity) furnished by the`):
print(`rational-function corrections supplied by `):
print(`RatImpCfD(pol,Cf,n,d2), d2=degree(pol,n)-1..L`):
print(`For example, try: SidraCf(n^2,1,n,10);`):
fi:


if nops([args])=1 and op(1,[args])=yafe then
print(`yafe(ope,N): given a recurrence  operator in the shift`):
print(`N, outputs it in nice form`):
fi:

if nops([args])=1 and op(1,[args])=Ztrans then
print(`Ztrans(pol,n,x): Given a polynomial pol, in n, finds`):
print(`the rational function sum(pol(i)*x^i,i=0..infinity)`):
fi:

end:


#Expl(P,N): Given a recurrence operator P(n,N) (or order I) outputs
# it in explicit form with f(n+I) expressed in terms
#of (n),...,f(n+I-1) For example, try Expl(N^2-N-1,n,N)
Expl:=proc(P,N)
local mu,I1:
I1:=degree(P,N):
mu:=coeff(P,N,I1):
expand((mu*N^I1-P)/mu):
end:

#Shiftr(P,n,N,r): If f is a solution of P(n,N)f=0, then
#the operator Q(n,N) or order degree(P,N)-1 such that
#N^rf=Q(n,N)f. For example, try: Shiftr(N^2-N-1,n,N,3);
Shiftr:=proc(P,n,N,r)
local I1,gu,mu:
I1:=degree(P,N):

if r<0 then
ERROR(` r>=0 `):
fi:
if r<I1 then
RETURN(N^r):
fi:

if r=I1 then
RETURN(Expl(P,N)):
fi:

gu:=Shiftr(P,n,N,r-1):
gu:=expand(subs(n=n+1,gu)*N):
mu:=coeff(gu,N,I1):
gu:=expand(gu-mu*N^I1+mu*Expl(P,N)):
gu:
end:


#ApplyOpe(P,Q,n,N): Given two recurrence operators P(n,N) and Q(n,N)
#finds the operator R(n,N) of order degree(P,N)-1 such that
#if f is a solution of P(n,N)f=0 then Q(n,N)f=R(n,N)f
#For example, try ApplyOpe(N^2-N-1,N^2-N-1,n,N)
ApplyOpe:=proc(P,Q,n,N)
local gu,i:
gu:=0:

for i from 0 to degree(Q,N) do
gu:=expand(gu+coeff(Q,N,i)*Shiftr(P,n,N,i)):
od:
gu:
end:


#ResultOpe1(Ope1,Ope2,n,N): Given two recurrence operators Ope1(n,N) and 
#Ope2(n,N)
#finds the operator Ope3(n,N) of the same order as Ope1 such that
#if f is a solution of Ope1(n,N)f=0 then Ope3(n,N)(Ope2(n,N)f)=0
#For example, try ResultOpe1(N^2-N-1,N+1,n,N) 
ResultOpe1:=proc(P,Q,n,N)
local ope,i,eq,var,a,I1,gu,var1:
I1:=degree(P,N):
ope:=0:
var:={}:
for i from 0 to I1 do
 ope:=ope+a[i]*N^i:
 var:=var union {a[i]}:
od:
gu:=0:
for i from 0 to I1 do
 gu:=expand(gu+a[i]*ApplyOpe(P,expand(subs(n=n+i,Q)*N^i),n,N)):
od:

eq:={}:
for i from 0 to I1 do
eq:=eq union {coeff(gu,N,i)}:
od:
var1:=solve(eq,var):
ope:=subs(var1,ope):
for i from 1 to nops(var) do
ope:=subs(var[i]=1,ope):
od:
ope:=numer(normal(ope)):
ope:=normal(ope/coeff(ope,N,I1)):
ope:=yafe(expand(ope),N):
if not type(coeff(ope,N,I1),rational) then
#print(`ope is`,ope):
RETURN(0):
fi:
ope/coeff(ope,N,2):
end:


#IniOpe(pol,x,n,N): The second-order operator,
#in the variable n and its shift operator N,  satisfied
#by the numerator and denominator of sum(x^i/pol(i),i=1..n)
#followed by the initial values for the numerator and deno.
#respectively
#For example, try IniOpe(n^2,1,n,N);
IniOpe:=proc(pol,x,n,N) local P,Q:
P:=-x*pol^2: Q:=x*pol+subs(n=n+1,pol):
N^2-subs(n=n+1,Q)*N-subs(n=n+1,P),[0,x],[1,subs(n=1,pol)]:
end:


#PhLogt(n,N,h,t): The accelerating operator of Apery for
#going from the h-th row to the next
PhLogt:=proc(n,N,h,t): N+(n-h)*t:end:

#Phz2(n,N,h): The accelerating operator of Apery for
#going from the h-th row to the next
Phz2:=proc(n,N,h): expand(2*(N-(n+1)^2+(h+1)*(n+1)-(h+1)^2/2)):end:

#Phz3(n,N,h): The accelerating operator of Apery for
#going from the h-th row to the next for Zeta(3)
Phz3:=proc(n,N,h): N-(n+1)^3+2*(h+1)*(n+1)^2-
2*(h+1)^2*(n+1)+(h+1)^3:end:


#OpehLogt(n,N,h,t): The recurrence operator annihilated by
#the h-th row for Apery's scheme for log(1+t)
OpehLogt:=proc(n,N,h,t): N^2-
((n+2)-(n+1)*t+h*(1+t))*N-
(n+1)^2*t:
end:

#Opehz2(n,N,h): The recurrence operator annihilated by
#the h-th row for Apery's scheme for zeta(2)
Opehz2:=proc(n,N,h): N^2-(2*(n+1)^2+2*(n+1)+h^2+h+1)*N+(n+1)^4:end:

#Opehz3(n,N,h): The recurrence operator annihilated by
#the h-th row for Apery's scheme for zeta(3)
Opehz3:=proc(n,N,h): N^2-
(2*(n+1)^3+3*(n+1)^2+(4*h^2+4*h+3)*(n+1)+(2*h^2+2*h+1))*N+
(n+1)^6:end:

#CheckOpehz2(h): Checks The recurrence operator annihilated by
#the h-th row for Apery's scheme for zeta(2)
CheckOpehz2:=proc(h) local n,N:
evalb(
expand(ResultOpe1(Opehz2(n,N,h), Phz2(n,N,h) ,n,N)-
Opehz2(n,N,h+1))=0
     ):
end:


#CODV(ope,n,m,N,M,g): Given a partial recurrence operator
#(with poly. coeffs.), ope in n,m, and the shift operators
#N and M, and a closed-form expression g (in n and m)
#outputs the operator annihilating g(n,m)F(n,m)
#where F(n,m) is a function annihilated by ope
#For example, try CODV(N-1,n,m,N,M,2^n)
CODV:=proc(ope,n,m,N,M,g)
local i,j,ope1,gu,mu,lu:
ope1:=expand(ope):
gu:=0:

for i from ldegree(ope1,N) to degree(ope1,N) do
lu:=coeff(ope1,N,i):
for j from ldegree(ope1,M) to degree(ope1,M) do
mu:=coeff(lu,M,j):
gu:=gu+mu*normal(simplify(g/subs({n=n+i,m=m+j},g)))*N^i*M^j:
od:
od:
gu:
end:

#Mult(P,Q,N,n): Given two recurrence operators P(N,n)
#and Q(N,n), (where n is the discrete variable and N is the
#shift operator N : Nx(n):=x(n+1), computes the
#product P(N,n)Q(N,n); For example, try:
#Mult(N,n^2*N,N,n)
Mult:=proc(P,Q,N,n)
local i,gu:
gu:=0:
for i from ldegree(P,N) to degree(P,N) do
gu:=gu+expand(coeff(P,N,i)*subs(n=n+i,Q)*N^i):
od:
collect(gu,N):
end:

#ShiftPQrs(P,Q,n,m,N,r,s): If f is a solution of P(n,m,N)f=0, 
#and Mf(m,n)=Q(n,m,N)f finds the operator R(n,m,N)
#of order degree(P,N)-1 such that
#N^rM^sf=R(n,m,N)f. For example, try: 
#ShiftPQrs(N^2-N-1,1+n*N,n,m,N,3,2);
ShiftPQrs:=proc(P,Q,n,m,N,r,s)
local gu,i,lu:
if s=0 then
RETURN(Shiftr(P,n,N,r) ):
fi:
gu:=ShiftPQrs(P,Q,n,m,N,r,s-1):
gu:=subs(m=m+1,gu):
lu:=Mult(gu,Q,N,n):

gu:=0:
for i from 0 to degree(lu,N) do
gu:=gu+coeff(lu,N,i)*Shiftr(P,n,N,i):
od:
collect(gu,N):
end:



#DiagPQ(Ope1,Ope2,n,m,N,M): Input: recurrence operators Ope1(n,m,N)
#and Ope2(n,m,N) (n,m, are variables and N and M are the shift
#operators in n,m, resp.); Output: an operator of the form
#Ope3(n,m,N*M) that annihilates any function that is annihilated 
#by Ope1(n,m,N) and M-Ope2(n,m,N) For example, try:
#DiagPQ(Opehz2(n,N,h) ,Phz2(n,N,h),n,h,N,H);
DiagPQ:=proc(P,Q,n,m,N,M)
local i,ope,eq,var,a,Ord,lu,K,var1:
Ord:=degree(P,N):
ope:=0:
var:={}:
for i from 0 to Ord do
 ope:=ope+a[i]*K^i:
 var:=var union {a[i]}:
od:

var1:=var:
lu:=0:
for i from 0 to Ord do
lu:=lu+a[i]*ShiftPQrs(P,Q,n,m,N,i,i):
od:
lu:=expand(lu):

eq:={}:

for i from 0 to Ord do
 eq:=eq union {coeff(lu,N,i)}:
od:

var:=solve(eq,var):

ope:=subs(var,ope):
ope:=normal(ope):
ope:=numer(ope):
ope:=yafe(ope,K):
ope:=subs(K=N*M,ope):

for i from 1 to nops(var1) do
ope:=subs(var1[i]=1,ope):
od:
ope:
end:



#DiagLogt(N,n,t): the recurrence satisfied by the diagonal
#of Apery's accelerating scheme for  log(1+t)
DiagLogt:=proc(N,n,t) local h:
MainDiagPQ(OpehLogt(n,N,h,t) ,PhLogt(n,N,h,t),
n,h,N,1/n!,1/h!): end:

#Diagz2(N,n): the recurrence satisfied by the diagonal
#of Apery's accelerating scheme for  Zeta(2)
Diagz2:=proc(N,n) local h:
MainDiagPQ(Opehz2(n,N,h) ,Phz2(n,N,h),n,h,N,1/n!^2,1/h!^2*2^h): end:

#MainDiagPQ(P,Q,n,m,N,L1,L2): Like DiagPQ, but 
#there is no input of the shift M, and instead
#there is a normalizing sequences L1(n) and L2(m)
#and the output is for the recurrence
#along the main diagonal using N
MainDiagPQ:=proc(P,Q,n,m,N,L1,L2)
local P1,Q1,M,R:
P1:=CODV(P,n,m,N,M,L1):
Q1:=CODV(Q,n,m,N,M,L1):
R:=DiagPQ(P1,Q1,n,m,N,M):
R:=CODV(R,n,m,N,M,L2):
R:=subs({M=1,m=n},R):
R:=yafe(R,N):
R:
end:


#yafe(ope,N): given a recurrence  operator in the shift
#N, outputs it in nice form
yafe:=proc(ope,N)
local ope1,ope2,i:
ope1:=expand(ope):
ope1:=ope1/coeff(ope1,N,degree(ope1,N)):
ope1:=normal(ope1):
ope1:=numer(ope1):
ope2:=0:

for i from ldegree(ope1,N) to degree(ope1,N) do
 ope2:=ope2+factor(coeff(ope1,N,i))*N^i:
od:
ope2:
end:

#Diagz3(N,n): the recurrence satisfied by the diagonal
#of Apery's accelerating scheme for  Zeta(3)
Diagz3:=proc(N,n) local h:
MainDiagPQ(Opehz3(n,N,h) ,Phz3(n,N,h),n,h,N,1/n!^3,1/h!^3): end:


#Atom1(x,i,j,D1): the differential operator f->(x d/dx)^i (x^(j)*f)
Atom1:=proc(x,i,j,D1)
local lu:
if i<0 then
ERROR(`Bad input`):
fi:
if i=0 then
 RETURN(x^j):
fi:
lu:=Atom1(x,i-1,j,D1):
expand(x*D1*lu+x*diff(lu,x)):
end:

#RecToDiff(ope,N,n,x,D1): given a recurrence operator
#ope(N,n) finds the differential operator P(x,D1) (D1=d/dx)
#annihilating
#the generating function f(x) of the sequence that is 
#annihilated by ope(N,n). For example, try RecToDiff(N-1,N,n,x,D1)
RecToDiff:=proc(ope,N,n,x,D1)
local gu,i,j,lu,lu1,lu11:
lu:=expand(ope):
gu:=0:
for j from ldegree(ope,N) to degree(ope,N) do
lu1:=coeff(lu,N,j):

 for i from 0 to degree(lu1,n) do
  lu11:=coeff(lu1,n,i):
  gu:=expand(gu+lu11*Atom1(x,i,-j,D1)):
 od:
od:
gu:
end:


#ApplyDiffOpeOld(ope,x,D1,F): applies the differential operator
#ope(x,D1) to F(x)
ApplyDiffOpeOld:=proc(ope,x,D1,F)
local gu,i,ope1,lu:
ope1:=expand(ope):
gu:=coeff(ope1,D1,0)*F:
lu:=F:

for i from 1 to degree(ope1,D1) do
 lu:=expand(diff(lu,x)):
 gu:=expand(gu+coeff(ope1,D1,i)*lu):
 
od:

normal(gu):
end:



#ApplyDiffOpe(ope,x,D1,F): applies the differential operator
#ope(x,D1) to F(x)
ApplyDiffOpe:=proc(ope,x,D1,F)
local gu,i,ope1,lu:
ope1:=expand(ope):
gu:=hafel(coeff(ope1,D1,0),x,F):
lu:=F:

for i from 1 to degree(ope1,D1) do
 lu:=expand(diff(lu,x)):
 gu:=expand(gu+hafel(coeff(ope1,D1,i),x,lu)):
 
od:

normal(gu):
end:


#GFz2qh(x,h): The generating function
#(1/h!^2)sum(qnh(n,h)/n!^2*x^n,n=0..infinity)
#implied by Phz2 
GFz2qh:=proc(x,h)
local lu,gu,MU,n,N,D1,H:
if h=0 then
 RETURN(1/(1-x)):
fi:
MU:=GFz2qh(x,h-1):
lu:=Phz2(n,N,h-1):
lu:=CODV(lu,n,h,N,H,1/n!^2):
gu:=RecToDiff(lu,N,n,x,D1):
MU:=ApplyDiffOpe(gu,x,D1,MU):
MU:=normal(MU/h^2):
MU:
end:


#GFz3qh(x,h): The generating function
#(1/h!^3)sum(qnh(n,h)/n!^3*x^n,n=0..infinity)
#implied by Phz3 
GFz3qh:=proc(x,h)
local lu,gu,MU,n,N,D1,H:
if h=0 then
 RETURN(1/(1-x)):
fi:
MU:=GFz3qh(x,h-1):
lu:=Phz3(n,N,h-1):
lu:=CODV(lu,n,h,N,H,1/n!^3):
gu:=RecToDiff(lu,N,n,x,D1):
MU:=ApplyDiffOpe(gu,x,D1,MU):
MU:=factor(normal(MU/h^3)):
MU:
end:



#GFLogtqh(x,h): The generating function
#(1/h!)sum(qnh(n,h)/n!*x^n,n=0..infinity)
#implied by PhLogt 
GFLogtqh:=proc(x,h,t)
local lu,gu,MU,n,N,D1,H:
if h=0 then
 RETURN(1/(1-x)):
fi:
MU:=GFLogtqh(x,h-1,t):
lu:=PhLogt(n,N,h-1,t):
lu:=CODV(lu,n,h,N,H,1/n!):
gu:=RecToDiff(lu,N,n,x,D1):
MU:=ApplyDiffOpe(gu,x,D1,MU):
MU:=factor(normal(MU/h)):
MU:
end:




#GFz2ph(x,h,F): The generating function
#(1/h!^2)sum(pnh(n,h)/n!^2*x^n,n=0..infinity)
#implied by Phz2  where F(x)=sum(x^i/i^2,i=1..infinity)
GFz2ph:=proc(x,h,F)
local lu,gu,MU,n,N,D1,H:
if h=0 then
 RETURN(F(x)/(1-x)):
fi:
MU:=GFz2ph(x,h-1,F):
lu:=Phz2(n,N,h-1):
lu:=CODV(lu,n,h,N,H,1/n!^2):
gu:=RecToDiff(lu,N,n,x,D1):
MU:=ApplyDiffOpe(gu,x,D1,MU):
MU:=normal(MU/h^2):
MU:
end:

#FindBGeneral(P,Q,n,L): Uses Apery's method to find a polynomial
#B such that B(n)*(Q(n)+B(n+1))-P(n) is of degree L
#possible. For example: try FindBGeneral(-n^4,2*n^2+2*n+1,n,0);
#This the general version allowing for symbolic or irrational coeffs
#and FindBGeneral(-n^4,2*n^2+2*n+3,n,0);
FindBGeneral:=proc(P,Q,n,L)
local eq,var,B,b,pol,i,var1:
var:={}:
eq:={}:
B:=0:
for i from 0 to degree(Q,n) do
B:=B+b[i]*n^i:
var:=var union {b[i]}:
od:

pol:=expand(B*(Q+subs(n=n+1,B))-P):
eq:={}:
for i from L+1 to degree(pol,n) do
eq:=eq union {coeff(pol,n,i)}:
od:
var1:={solve(eq,var)}:
{seq(subs(var1[i],B),i=1..nops(var1))}:
end:

#FindB1(P,Q,n,L): Uses Apery's method to find a polynomial
#B such that B(n)*(Q(n)+B(n+1))-P(n) is of degree L
#possible. For example: try FindB1(-n^4,2*n^2+2*n+1,n,0);
#and FindB1(-n^4,2*n^2+2*n+3,n,0);
FindB1:=proc(P,Q,n,L)
local eq,var,B,b,pol,i,var1,gu,B1:
var:={}:
eq:={}:
B:=0:
for i from 0 to max(degree(P,n),degree(Q,n)) do
B:=B+b[i]*n^i:
var:=var union {b[i]}:
od:

pol:=expand(B*(Q+subs(n=n+1,B))-P):
eq:={}:
for i from L+1 to degree(pol,n) do
eq:=eq union {coeff(pol,n,i)}:
od:
var1:={solve(eq,var)}:

gu:={}:

for i from 1 to nops(var1) do
B1:=subs(var1[i],B):
if IsRatPol(B1,n) then
 gu:=gu union {B1}:
fi:
od:
gu:
end:

#PQtoOpe(P,Q,N,n): Given functions of n, P(n) and Q(n)
#constructs the operator: N^2-Q(n+1)*N-P(n+1)
PQtoOpe:=proc(P,Q,N,n)
N^2-subs(n=n+1,Q)*N-subs(n=n+1,P)
end:



#IsRatPol(P): Given a polynomial P in n decides whether
#all its coefficients are rational numbers 
IsRatPol:=proc(P)
local i,P1,P11:
P1:=expand(P):
if not type(P1,`+`) then
 if type(evalf(op(1,P1)),numeric) then
  if not type(op(1,P1),rational) then
      RETURN(false):
  fi:
 fi:
RETURN(true):
fi:

for i from 1 to nops(P1) do
  P11:=op(i,P1):
  if type(evalf(op(1,P11)),numeric) then
  if not type(op(1,P11),rational) then
      RETURN(false):
  fi:
 fi:
od:
true:
end:

#IsRatPolOld(P,n): Given a polynomial P in n decides whether
#all its coefficients are rational numbers 
IsRatPolOld:=proc(P,n)
local i:
for i from 0 to degree(P,n) do
 if not type(coeff(P,n,i),rational) then
    RETURN(false):
 fi:
od:
true:
end:

#FindB(P,Q,n,L): Uses Apery's method to find a polynomial
#B such that B(n)*(Q(n)+B(n+1))-P(n) is of degree L
#possible. For example: try FindB(-n^4,2*n^2+2*n+1,n);
#and FindB(-n^4,2*n^2+2*n+3,n);
FindB:=proc(P,Q,n)
local L,gu:

for L from 0 to degree(P,n)-1 do
 gu:=FindB1(P,Q,n,L):
 if gu<>{} then
   RETURN(gu):
 fi:
od:
{}:
end:



#GuessPol(List,h): Given a list of expressions, List,
#guesses a polynomial P,of degree<=nops(List)-3 such that
#it fits List=[P(0),P(1),P(nops(List)-1)], and returns
#0 if none exists
GuessPol:=proc(List,h)
local pol,eq,var,i,a:

if nops(List)-3<0 then
 ERROR(`The first input should have at least length 3`):
fi:
var:={}:
pol:=0:
for i from 0 to nops(List)-3 do
pol:=pol+a[i]*h^i:
var:=var union {a[i]}:
od:
eq:={}:
for i from 1 to nops(List) do
eq:=eq union {expand(subs(h=i-1,pol)-List[i])}:
od:
var:=solve(eq,var):
if var=NULL then
 RETURN(0):
fi:
expand(subs(var,pol)):
end:


#GuessPolNew(List,h,h0,deg): Given a list of expressions, List,
#guesses a polynomial P,of degree deg such that
#it fits List=[P(h0),P(h0+1),P(h0+nops(List)-1)], and returns
#0 if none exists
GuessPolNew:=proc(List,h,h0,deg)
local pol,eq,var,i,a:

var:={}:
pol:=0:
for i from 0 to deg do
pol:=pol+a[i]*h^i:
var:=var union {a[i]}:
od:
eq:={}:
for i from 1 to nops(List) do
eq:=eq union {expand(subs(h=h0+i-1,pol)-List[i])}:
od:
var:=solve(eq,var):
if var=NULL then
 RETURN(0):
fi:
expand(subs(var,pol)):
end:


#Aperyh1(pol,x,n,R,h):Given a polynomial pol in the variable n
#with polynomial coefficient, uses Apery's ORIGINAL method
#and a variable h tries, if possible, to fit the
#accelerating B(n)'s and the coeffs. of the recurrence
#for row h by polynomials of degree h
#by fitting the data from Apery(pol,n,R)
#It returns B(n,h),and [P(n,h),Q(n.h)] if successful
#otherwise it returns 0. For example, try
#Aperyh1(n^2,1,n,8,h); and Aperyh1(n^3,1,n,8,h);
Aperyh1:=proc(pol,x,n,R,h)
local lu,B,P,Q,i,ope,N,ope1:
lu:=Apery(pol,x,n,R):

if lu=0 then
 RETURN(0):
fi:

B:=GuessPol(lu[1],h):
P:=GuessPol(lu[2],h):
Q:=GuessPol(lu[3],h):

if B=0 or P=0 or Q=0 then
RETURN(0):
fi:

ope:=N^2-subs(n=n+1,Q)*N-subs(n=n+1,P):
ope1:=ResultOpe1(ope,N+subs(n=n+1,B),n,N):
if not type(normal(ope1/subs(h=h+1,ope)),integer) then
 print(`The conjectured B,P,Q did not make it`,ope1,subs(h=h+1,ope)):
 RETURN(0):
fi:
B,P,Q:
end:

#qn(n,pol,n0): Given a polynomial in n, computes
#the initial denominators in the partial sum of
#1/pol(1)+1/pol(2)+...+1/pol(n0)
qn:=proc(n,pol,n0) option remember:
if n0=0 then
1:
else
qn(n,pol,n0-1)*subs(n=n0,pol):
fi:
end:

#pn(n,pol,x,n0): Given a polynomial in n, computes
#the initial numerator in the partial sum of
#1/pol(1)+1/pol(2)+...+1/pol(n0)
pn:=proc(n,pol,x,n0) option remember:
if n0=0 then
0:
else
pn(n,pol,x,n0-1)+x^n0/subs(n=n0,pol):
fi:
end:


#qnh(n,h,pol,B,n0,h0) Given a polynomial in the variable n
#and an accelerating polynomial B(n,h), computes the denominator
#in the approximating scheme for the n0-th entry in row h0
#for example, try: qnh(n,h,n^2,-1/2+n-n^2+h*n-h-1/2*h^2,2,2);
qnh:=proc(n,h,pol,B,n0,h0) option remember:
if h0=0 then
  1:
else
qnh(n,h,pol,B,n0,h0-1)*subs({n=n0+1,h=h0-1},B)+
subs(n=n0+1,pol)*qnh(n,h,pol,B,n0+1,h0-1):
fi:
end:

#pnh(n,h,pol,x,B,n0,h0) Given a polynomial in the variable n
#and an accelerating polynomial B(n,h), computes the numerator
#in the approximating scheme for the n0-th entry in row h0
#for example, try: pnh(n,h,n^2,1,-1/2+n-n^2+h*n-h-1/2*h^2,2,2);
pnh:=proc(n,h,pol,x,B,n0,h0) option remember:
if h0=0 then
  pn(n,pol,x,n0):
else
pnh(n,h,pol,x,B,n0,h0-1)*subs({n=n0+1,h=h0-1},B)+
subs(n=n0+1,pol)*pnh(n,h,pol,x,B,n0+1,h0-1):
fi:
end:
 





#AperySeq(pol,x,n,L): 
# Returns the sequence of diagonal approximants
#up to the L-th term,
AperySeq:=proc(pol,x,n,L)
local lu,h,lu1,lu2,n0,B,i:

lu:=Aperyh(pol,x,n,h):

if lu=0 then
RETURN(0):
fi:

B:=lu[1]:
lu1:=[seq(pnh(n,h,pol,x,B,n0,n0),n0=0..L)]:
lu2:=[seq(qnh(n,h,pol,B,n0,n0),n0=0..L)]:
[seq(lu1[i]/lu2[i],i=1..L+1)]:
end:


#delta(sid): Given a sequence of rational numbers, sid
#with at least 31 terms
#finds the smallest delta in the list s.t. |K-an/bn|<=C/bn^(1+delta)
#n<nops(sid)-10
#where K is its last element
#delta([seq(fibonacci(i+1)/fibonacci(i),i=1..100)]):
delta:=proc(sid)
local n,mu,i,K:
Digits:=200:
n:=nops(sid):
if n<31 then
ERROR( `the number of terms should be at least 15`):
fi:
K:=sid[n]:
mu:=[]:
for i from 15 to n-15 do
mu:=[op(mu),evalf(-log(abs(K-sid[i]))/log(denom(sid[i])))-1]:
od:
min(op(mu)):
end:


#CODVpol(ope,n,N,pol): Given a partial recurrence operator
#(with poly. coeffs.), ope in n and the shift operator
#N, and a polynomial pol(in n) 
#outputs the operator annihilating F(n)/(pol(1)...pol(n))
#where F(n) is a function annihilated by ope
#For example, try CODVpol(N-2*(n+1),n,N,n)
CODVpol:=proc(ope,n,N,pol)
local i,ope1,gu,mu,lu,gu1,lu1:
ope1:=expand(ope):
gu:=0:
mu:=1:
for i from 0 to degree(ope1,N) do
lu:=coeff(ope1,N,i):
mu:=mu*subs(n=n+i,pol):
gu:=gu+lu*mu*N^i:
od:
gu:=expand(gu):
lu:=coeff(gu,N,degree(gu,N)):
for i from 0 to degree(gu,N)-1 do
lu1:=coeff(gu,N,i):
if lu1<>0 then
lu:=gcd(lu,lu1):
fi:
od:
gu1:=0:
for i from 0 to degree(gu,N) do
gu1:=gu1+factor(normal(coeff(gu,N,i)/lu))*N^i:
od:
gu1:
end:


#SeqFromRec(ope,n,N,ini,L): Given an (ordinary) recurrence operator
#ope(n,N) in the variable n and the shift operator N (Nf(n):=f(n+1))
#and given initial values [ini0,ini1,..,ini_{ORD-1}],computes
#the first L+1 terms of the sequence a(n) satisfying
#ope(n,N)a(n)=0 and a[i]=ini[i] for i=0,...,ORD-1
#For example, try: SeqFromRec(N^2-N-1,n,N,[0,1],10);
SeqFromRec:=proc(ope,n,N,ini,L)
local gu,lu,j,ORD,n0,i:

ORD:=degree(ope,N):
if ORD<>nops(ini) then
 ERROR(`The length of`,ini, `should be `, degree(ope,N)):
fi:
gu:=ini:
for i from ORD to L do
lu:=0:
n0:=i-ORD:
lu:=0:
for j from 0 to ORD-1 do
lu:=lu+subs(n=n0,coeff(ope,N,j))*gu[n0+j+1]:
od:
lu:=-lu/subs(n=n0,coeff(ope,N,ORD)):
gu:=[op(gu),lu]:
od:
gu:
end:

#SeqQuotFromRec(ope,ini1,ini2,n,N,L): 
#Given an (ordinary) recurrence operator
#ope(n,N) in the variable n and the shift operator N (Nf(n):=f(n+1))
#and given initial values ini1,ini2 finds
#the first L+1 terms of the sequence a(n)/b(n) where both 
#a(n) and b(n) are annihilated by ope(n,N)
#and a[i]=ini1[i] for i=0,...,ORD-1 ; 
#b[i]=ini2[i] for i=0,...,ORD-1
#For example, try: SeqQuotFromRec(N^2-N-1,[0,1],[1,1],n,N,10);
SeqQuotFromRec:=proc(ope,ini1,ini2,n,N,L)
local lu1,lu2,i:
lu1:=SeqFromRec(ope,n,N,ini1,L):
lu2:=SeqFromRec(ope,n,N,ini2,L):
[seq(lu1[i]/lu2[i],i=1..L+1)]:
end:



#Appx(pol,x,n,L): Gives the Lth diag. approximant
#Apery's recurrence for x/pol(1)+x^2/pol(2)+x^3/pol(3)+... 
#followed by the error. For example: try Appx(n^2,1,n,20);
Appx:=proc(pol,x,n,L)
local lu,ku,h,N,lu1p,lu1q:
lu:=AccRec(pol,x,n,N):

if lu=0 then
print(`No recurrence found`):
RETURN(0):
fi:
ku:=SeqQuotFromRec(lu,n,N,L):
ku[L],evalf(ku[L]-ku[L-1]):
end:

#Apery(pol,x,n,R):Given a polynomial pol in the variable n
#with polynomial coefficient, uses Apery's ORIGINAL method
#for TRYING to prove irrationality for the constant
#sum(x^i/pol(i),i=1..infinity)
#as described in his beautiful C.T.H.S. paper. It continues
#to accelerate up to R rows, if possible, and returns
#the list of accelerating B's, the list of P's, the list of Q's
#for the recurrences for the
#rows, and the expression themselves of the
#"easy" solution in row i divided by pol(1)pol(2)...pol(n)
#If it fails, it returns 0
#For example, try Apery(n^2,1,n,5);, Apery(n^3,1,n,5);
Apery:=proc(pol,x,n,R)
local ListB,ListP,ListQ,ListSol,P,Q,B,ope,Sol1,lu,Sol2,Sol2a,i1,N,ope1,i:
ListB:=[]:
ListSol:=[1]:
Sol1:=1:
P:=-x*pol^2:
Q:=x*pol+subs(n=n+1,pol):
ListP:=[]:
ListQ:=[]:
for i from 1 to R do


ListP:=[op(ListP),P]:
ListQ:=[op(ListQ),Q]:
lu:=FindB1(P,Q,n,0):

ope:=expand(N^2-subs(n=n+1,Q)*N-subs(n=n+1,P)):
if nops(lu)<=1 then
# print(`Nothing found`):
  RETURN(0):
fi:

B:=lu[1]:

Sol2:=expand(subs(n=n+1,Sol1)*subs(n=n+1,pol)+Sol1*subs(n=n+1,B)):

for i1 from 2 to nops(lu) do
Sol2a:=expand(subs(n=n+1,Sol1)*subs(n=n+1,pol)+Sol1*subs(n=n+1,lu[i1])):
 if degree(Sol2a,n)>degree(Sol2,n) then
   Sol2:=Sol2a:
   B:=lu[i1]:
 fi:
od:
ListB:=[op(ListB),B]:
if nops(ListB)>1 and
degree(ListB[nops(ListB)],n)>degree(ListB[nops(ListB)-1],n) then
print(`ListB is`,ListB):
print(degree(ListB[nops(ListB)],n),degree(ListB[nops(ListB)-1],n)):
 RETURN(0):
fi:

Sol1:=Sol2:
ListSol:=[op(ListSol),Sol1]:
ope1:=expand(ResultOpe1(ope,N+subs(n=n+1,B),n,N)):

if ope1=0 then
RETURN(0):
fi:

P:=expand(subs(n=n-1,-coeff(ope1,N,0))):
Q:=expand(subs(n=n-1,-coeff(ope1,N,1))):
od:

expand(ListB),expand(ListP),expand(ListQ),expand(ListSol):

end:


#AccRec(pol,x,n,N): Given a polynomial pol
#in the variable n, it tries to uses Apery's method to 
#find a convergence-acceleration recurrence (using the diagonal)
#followed by the initial conditions [[a0,a1],[b0,b1]]
#such that the constant sum(1/pol(n),n=1..infinity) is
#approximated by the quotients of solution of the recurrence
#with the given initial conditions for the numerator and
#denominator. For example, try:
#AccRec(n^2,1,n,N);
AccRec:=proc(pol,x,n,N)
local lu,B,P,Q,ope,lu1,lu2,n0,ope2,ope3,h,degDiff,c,p,r,ku:

lu:=Aperyh(pol,x,n,h):

if lu=0 then
print(`Nothing found`):
RETURN(0):
fi:
B:=lu[1]:
P:=lu[2]:
Q:=lu[3]:
ope:=N^2-subs(n=n+1,Q)*N-subs(n=n+1,P):

lu1:=[seq(pnh(n,h,pol,x,B,n0,n0),n0=0..1)]:
lu2:=[seq(qnh(n,h,pol,B,n0,n0),n0=0..1)]:
ope2:=N+subs(n=n+1,B):
ope3:=MainDiagPQ(ope,ope2,n,h,N,1,1):
ope3:=CODVpol(ope3,n,N,pol):
ope3:=yafe(ope3,N):
lu1:=[lu1[1],lu1[2]]:
lu2:=[lu2[1],lu2[2]]:
degDiff:=(degree(coeff(ope3,N,0),n)-degree(coeff(ope3,N,2),n))/2:

if degDiff>0 and type(degDiff,integer) then
ope3:=CODVpol(ope3,n,N,n^degDiff):
fi:

p:=expand(coeff(ope3,N,2)):
r:=expand(coeff(ope3,N,0)):
if degree(p,n)<>degree(r,n) then
ERROR(`Something is wrong`):
fi:

c:=sqrt(abs(normal(coeff(r,n,degree(r,n))/coeff(p,n,degree(p,n))))):

if c<>1 and type(c,rational) then
ope3:=coeff(ope3,N,2)*N^2*c^2+coeff(ope3,N,1)*c*N+coeff(ope3,N,0):
p:=expand(coeff(ope3,N,2)):
ku:=coeff(p,n,degree(p,n)):
ope3:=ope3/ku:
ope3:=yafe(ope3,N):
lu1:=[lu1[1],lu1[2]/c]:lu2:=[lu2[1],lu2[2]/c]:
fi:
ope3,lu1,lu2:
end:


#delt(a,c): Given a rational number a and a constant
#c, finds the delta such that |a-c|=1/denom(a)^(1+delta)
#For example, try delta(22/7,evalf(Pi)):
delt:=proc(a,c): evalf(-log(abs(c-a))/log(denom(a)))-1: end:


#Aperyh(pol,x,n,h):Given a polynomial pol in the variable n
#with polynomial coefficient, uses Apery's ORIGINAL method
#and a variable h tries, if possible, to fit the
#accelerating B(n)'s and the coeffs. of the recurrence
#for row h by polynomials of degree h
#by fitting the data from Apery(pol,n,R)
#It returns B(n,h),and [P(n,h),Q(n.h)] if successful
#otherwise it returns 0. For example, try
#Aperyh(n^2,1,n,h); and Aperyh(n^3,1,n,h);
Aperyh:=proc(pol,x,n,h)
local R,lu,lu1:
option remember:
lu:=Apery(pol,x,n,6):

if lu=0 then
print(`Nothing Found`):
 RETURN(0):
fi:

for R from 4 to 10 while Aperyh1(pol,x,n,R,h)=0 do
od:
lu:=Aperyh1(pol,x,n,R+1,h):
lu1:=Aperyh1(pol,x,n,R+2,h):

if lu<>lu1 then
print(`Nothing found`):
RETURN(0):
fi:
lu:
end:


#Qnh0(pol,B,n,h,h0): an exact expression, in terms of the symbol n
#for qnh in Apery's scheme for h=h0. For example, try:
#Qnh0(n^2,Aperyh(n^2,1,n,h)[1],n,h,3):
Qnh0:=proc(pol,B,n,h,h0)
local gu:
option remember:
if h0=0 then
 RETURN(1):
fi:
gu:=Qnh0(pol,B,n,h,h0-1):
expand(gu*subs({n=n+1,h=h0-1},B)+subs(n=n+1,gu)*subs(n=n+1,pol)):
end:


#PolNiceBasis(pol,n): Given a polynomial pol in the variable
#n, expresses it in terms of binomial(n+i,d), (i=0..d)
#where d is the degree, e.g., try: PolNiceBasis(n^2+n,n);
PolNiceBasis:=proc(pol,n)
local eq,var,POL,a,i,d,POL1,lu:
d:=degree(pol,n):
POL:=0:
var:={}:

for i from 0 to d do
 POL:=POL+a[i]*binomial(n+i,d):
 var:=var union {a[i]}:
od:

POL1:=expand(expand(POL))-expand(pol):
eq:={}:
for i from 0 to d do
eq:=eq union {coeff(POL1,n,i)}:
od:
var:=solve(eq,var):
POL:=subs(var,POL):
lu:=coeff(POL,binomial(n,d),1):
[lu,POL/lu]:
end:


#Aitken(Lis1): The Aitken transform of the sequence Lis1
Aitken:=proc(Lis1)
local gu,n:
gu:=[]:

for n from 3 to nops(Lis1) do
gu:=[op(gu),(Lis1[n]*Lis1[n-2]-Lis1[n-1]^2)/(Lis1[n]-2*Lis1[n-1]+Lis1[n-2])]:
od:

gu:
end:

#Rnh0(pol,B,n,h,h0): an exact expression, in terms of the symbol n
#for r(n,h) in Apery's scheme for h=h0.
#Where r(n,h) is the polynomial in h such that
#p(n,h)=(1/pol(1)+...+1/pol(n))*q(n,h)+r(n,h)
#For example, try: Rnh0(n^2,Aperyh(n^2,1,n,h)[1],n,h,3):
Rnh0:=proc(pol,B,n,h,h0)
local gu:
option remember:
if h0=0 then
 RETURN(0):
fi:
gu:=Rnh0(pol,B,n,h,h0-1):
expand(gu*subs({n=n+1,h=h0-1},B)+
subs(n=n+1,gu)*subs(n=n+1,pol)
+subs(n=n+1,Qnh0(pol,B,n,h,h0-1))):
end:

#ApplyRec(ope,n,N,Lis): applies the recurrence ope(N,n) to the 
#list Lis, For example, try: ApplyRec(N^2-N-1,n,N,[1,1,2,3,5,8]);
ApplyRec:=proc(ope,n,N,Lis)
local ORD,lu,i,j,mu:
ORD:=degree(ope,N):
mu:=[]:
for i from 1 to nops(Lis)-ORD do
lu:=0:
 for j from 0 to ORD do
   lu:=lu+Lis[i+j]*subs(n=i-1,coeff(ope,N,j)):
 od:
mu:=[op(mu),lu]:
od:
mu:
end:


#Ztrans(pol,n,x): Given a polynomial pol, in n, finds
#the rational function sum(pol(i)*x^i,i=0..infinity)
Ztrans:=proc(pol,n,x)
local gu,i,lu,mu:
if pol=0 then
 RETURN(0):
fi:

lu:=1/(1-x):
gu:=0:
for i from 0 to degree(pol,n) do
mu:=coeff(pol,n,i):
gu:=normal(gu+mu*lu):
lu:=x*diff(lu,x):
od:
gu:
end:

#GF1Q(pol,x0,n,x,L): The list of generating functions for the first
#L+1 rows of the Apery scheme for 1/pol(1)+ ...+1/pol(i)+...
#for the qnh
GF1Q:=proc(pol,x0,n,x,L)
local h,B,h0:
B:=Aperyh(pol,x0,n,h):
if B=0 then
RETURN(0):
fi:
B:=B[1]:
[seq(factor(Ztrans(Qnh0(pol,B,n,h,h0),n,x)),h0=0..L)]:
end:



#GF1R(pol,x0,n,x,L): The list of generating functions for the first
#L+1 rows of the Apery scheme for 1/pol(1)+ ...+1/pol(i)+...
#for the qnh
GF1R:=proc(pol,x0,n,x,L)
local h,B,h0:
B:=Aperyh(pol,x0,n,h):
if B=0 then
RETURN(0):
fi:
B:=B[1]:
[seq(factor(Ztrans(Rnh0(pol,B,n,h,h0),n,x)),h0=0..L)]:
end:




hafel:=proc(pol,x,f)
local i,gu,i1,j,mu,lu,POL:
gu:=0:
POL:=pol:
for i from ldegree(pol,x) to -1 do
POL:=POL-coeff(pol,x,i)*x^i:
i1:=-i:
mu:=f:
lu:=taylor(f,x=0,i1+1):
for j from 0 to i1-1 do
mu:=mu-coeff(lu,x,j)*x^j:
od:
gu:=normal(gu+coeff(pol,x,i)*mu*x^i):
od:

normal(gu+POL*f):
end:

#OpeDh(pol,n,x,D1,h0) The differential operator that yields
#the generating function for the (h0-1)th row to the h0^th row
#for Apery's scheme for 1/pol(1)+...+1/pol(n)+...
OpeDh:=proc(pol,n,x,D1,h0)
local ope,B,h,N:
B:=Aperyh(pol,1,n,h):
if B=0 then
RETURN(0):
fi:
B:=B[1]:
ope:=subs({n=n+1,h=h0-1},B)+subs(n=n+1,pol)*N:
ope:=RecToDiff(ope,N,n,x,D1):
ope:
end:


#GF1Qdirect(pol,n,x,L): Like GF1Q but directly in terms of
#differential operator
GF1Qdirect:=proc(pol,n,x,L)
local i,h,D1,ope,gu,resh,ope0:
ope:=OpeDh(pol,n,x,D1,h):
resh:=[1/(1-x)]:
gu:=1/(1-x):
for i from 1 to L do
ope0:=subs(h=i,ope):
gu:=ApplyDiffOpe(ope0,x,D1,gu):
resh:=[op(resh),gu]:
od:
resh:
end:



#ilcmP(pol,n,n0): The least common multiple of pol(1), ..., pol(n0)
#where pol is a polynomial in n 
ilcmP:=proc(pol,n,n0)
local gu:
option remember:
if not type(n0,integer) or n0<=0 then
ERROR(`n0 integer >=0`):
fi:
if n0=1 then
 RETURN(subs(n=1,pol)):
fi:

gu:=ilcmP(pol,n,n0-1):

ilcm(gu,subs(n=n0,pol)):
end:

#AccRecN(pol,x0,n,N,g): Like AccRec(pol,x0,n,N,g), but with
#a normalizing function g entered
#For example, try:
#AccRecN(n^2,1,n,N,2^n/n!^2); AccRecN(n^3,1,n,N,1/n!^3);
#It also returns the first 11 terms of the numerator and
#denominator sequences
AccRecN:=proc(pol,x0,n,N,g)
local lu,B,P,Q,ope,lu1,lu2,n0,ope2,ope3,h,M,m,gu,gu1,gu2,i:

lu:=Aperyh(pol,x0,n,h):

if lu=0 then
print(`Nothing found`):
 RETURN(0):
fi:

B:=lu[1]: P:=lu[2]: Q:=lu[3]:

ope:=N^2-subs(n=n+1,Q)*N-subs(n=n+1,P):

lu1:=[seq(pnh(n,h,pol,x0,B,n0,n0),n0=0..1)]:
lu2:=[seq(qnh(n,h,pol,B,n0,n0),n0=0..1)]:
ope2:=N+subs(n=n+1,B):
ope3:=MainDiagPQ(ope,ope2,n,h,N,1,1):
ope3:=CODVpol(ope3,n,N,pol):
ope3:=yafe(ope3,N):
lu1:=[lu1[1]*subs(n=0,g),lu1[2]*subs(n=1,g)]:
lu2:=[lu2[1]*subs(n=0,g),lu2[2]*subs(n=1,g)]:
ope3:=CODV(ope3,n,m,N,M,g):
ope3:=expand(ope3):
ope3:=yafe(ope3,N):
gu:=expand(ope3):
gu:=coeff(gu,n,degree(gu,n)):
gu:=fsolve(gu):
if gu<>NULL then
gu:=log(max(gu)):
else
gu:=0:
fi:


gu1:=SeqFromRec(ope3,n,N,lu1,10):
gu2:=SeqFromRec(ope3,n,N,lu2,10):
gu1:=[seq(gu1[i]*ilcmP(pol,n,i),i=1..nops(gu1))]:
ope3,lu1,lu2,gu1,gu2,gu:
end:



#CheckN(pol,x0,n,h,g,L): Given a polynomial pol, in the variable
#n, a no. x0, and another variable h, and a normalization function g
#checks whether q(n,h)*g(n,h) are integers and
#lcm(pol(1), ..., pol(n))*p(n,h)*g(n,h) are integers
#For example, try: CheckN(n^3,n,h,1/h!^3,20);
CheckN:=proc(pol,x0,n,h,g,L)
local n0,h0,B:
B:=Aperyh(pol,x0,n,h);

if B=0 then
 RETURN(0):
fi:

B:=B[1]:

for h0 from 1 to L do
for n0 from h0 to L do

if not type(qnh(n,h,pol,B,n0,h0)*subs({n=n0,h=h0},g),integer) then
 print(`n=`,n0, `h=`,h0, `is a violation for qnh`):
 RETURN(false):
fi:

if not 
type(pnh(n,h,pol,x0,B,n0,h0)*subs({n=n0,h=h0},g)*ilcmP(pol,n,n0),integer)
       then
 print(`n=`,n0, `h=`,h0, `is a violation for pnh`):
RETURN(false):
fi:
od:
od:
true:
end:




#Lucky(n,x,deg,a,b): Given a variable n, and letters a,b, an integer
#deg, finds the general monic polynomial pol of dergee 
#deg that allows Apery-acceleration for sum(1/pol(i),i=1..infinity)
#For example, try: Lucky(n,1,2,a,b); 
Lucky:=proc(n,x,deg,a,b)
local eq,var,B,pol,i,var1,gu,B1,pol1,P,Q,pol1a,pol1aTry,i1:
pol1:=n^deg:
var:={}:
for i from 0 to deg-1 do
 pol1:=pol1+a[i]*n^i:
 var:=var union {a[i]}:
od:

P:=-(x*subs(n=n+1,pol1)+subs(n=n+2,pol1)):
Q:=x*subs(n=n+1,pol1)^2:

eq:={}:
B:=0:
for i from 0 to max(degree(P,n),degree(Q,n)) do
B:=B+b[i]*n^i:
var:=var union {b[i]}:
od:

pol:=expand(B*(P-subs(n=n+1,B))-Q):
eq:={}:
for i from 1 to degree(pol,n) do
eq:=eq union {coeff(pol,n,i)}:
od:
var1:={solve(eq,var)}:
gu:={}:

for i from 1 to nops(var1) do
 B1:=subs(var1[i],B):
 pol1a:=subs(var1[i],pol1):

pol1aTry:=pol1a:
for i1 from 1 to nops(var) do
 pol1aTry:=subs(var[i1]=1/(i1^2+1/3),pol1a):
od:

if diff(expand(B1+subs(n=n+1,pol1a)),n)<>0 and 
Apery(pol1aTry, x,n,2)<>0 then
gu:=gu union {pol1a}:
fi:
od:
gu:
end:




#Lucky1(n,x,deg,a,b): Given a variable n, and letters a,b, an integer
#deg, finds the general monic polynomial pol of dergee 
#deg that allows Apery-acceleration for sum(1/pol(i),i=1..infinity)
#For example, try: Lucky(n,1,2,a,b); 
Lucky1:=proc(n,x,deg,a,b)
local eq,var,B,pol,i,var1,gu,B1,pol1,P,Q,pol1a,pol1aTry,i1:
pol1:=n^deg:
var:={}:
for i from 0 to deg-1 do
 pol1:=pol1+a[i]*n^i:
 var:=var union {a[i]}:
od:

P:=-(x*subs(n=n+1,pol1)+subs(n=n+2,pol1)):
Q:=x*subs(n=n+1,pol1)^2:

eq:={}:
B:=0:
for i from 0 to max(degree(P,n),degree(Q,n)) do
B:=B+b[i]*n^i:
var:=var union {b[i]}:
od:

pol:=expand(B*(P-subs(n=n+1,B))-Q):
eq:={}:
for i from 1 to degree(pol,n) do
eq:=eq union {coeff(pol,n,i)}:
od:
var1:={solve(eq,var)}:
gu:={}:


for i from 1 to nops(var1) do
 B1:=subs(var1[i],B):
 pol1a:=subs(var1[i],pol1):
gu:=gu union {pol1a}:
od:
gu:
end:



#ResultOpe(Ope1,IniTop,IniBot,Ope2,n,N): 
#Given two recurrence operators Ope1(n,N) and 
#Ope2(n,N), and lists IniTop, IniBot, for the
#initial values
#finds the operator Ope3(n,N) of the same order as Ope1 such that
#if f is a solution of Ope1(n,N)f=0 then Ope3(n,N)(Ope2(n,N)f)=0
#followed by the initial values for the top and bottom
#of the sequence Ope2(n,N)f
#For example, try ResultOpe(N^2-N-1,[0,1],[1,0],N+1,n,N) 
ResultOpe:=proc(P,IniTop,IniBot,Q,n,N)
local ope,i,eq,var,a,I1,gu,var1:
I1:=degree(P,N):
ope:=0:
var:={}:
for i from 0 to I1 do
 ope:=ope+a[i]*N^i:
 var:=var union {a[i]}:
od:
gu:=0:
for i from 0 to I1 do
 gu:=expand(gu+a[i]*ApplyOpe(P,expand(subs(n=n+i,Q)*N^i),n,N)):
od:

eq:={}:
for i from 0 to I1 do
eq:=eq union {coeff(gu,N,i)}:
od:
var1:=solve(eq,var):
ope:=subs(var1,ope):
for i from 1 to nops(var) do
ope:=subs(var[i]=1,ope):
od:
if ApplyIni(Q,P,n,N,IniBot)=[0,0] then
 RETURN(0):
fi:
collect(numer(normal(ope)),N),
ApplyIni(Q,P,n,N,IniTop),
ApplyIni(Q,P,n,N,IniBot):
end:






#ApplyIni(ope,Q,n,N,Lis1): Given a recurrence operator 
#ope(n,N) and another one Q(n,N)
#in the variable n and its corresponding shift N
#and a list Lis1 of (size degree(Q,N)) of initial values
#finds the initial values for the sequence ope(n,N) applied to
#the sequence whose inital values are given by Lis1 
#and annihilated by Q(n,N) For example, try:
#ApplyIni(N+2,N-1,n,N,[1]);
ApplyIni:=proc(ope,Q,n,N,Lis1)
local gu:
if degree(Q,N)<>nops(Lis1) then
ERROR(`The order of`, Q, `should equal the numebr of elements in`, Lis1):
fi:
gu:=SeqFromRec(Q,n,N,Lis1,nops(Lis1)+degree(ope,N)):
gu:=ApplyRec(ope,n,N,gu):
[op(1..nops(Lis1),gu)]:
end:



#pBfromC(C,n): Gives a generic good pair (p,B) from C
#For example, try pBfromC(2*n*(n-1),n)
pBfromC:=proc(C,n)
local p:
p:=normal((subs(n=n+1,C)*C+subs(n=n+1,C)+C)/(subs(n=n+1,C)-C)):
p,normal(1-p+C):
end:



#FindAcc(pol,n): Uses Apery's method to find a rational
#function R(n) such that sum(1/pol(i),i=1..n)+R(n)
#is a better appx. to sum(1/pol(i),i=1..infinity) then
#the the partial sum itself
#For example, try FindAcc(n^4+7/16,n,N);
FindAcc:=proc(pol,n)
local B,h:
B:=Apery(pol,1,n,1):
if B=0 then
RETURN(0):
fi:
B:=B[1][1]:
Rnh0(pol,B,n,h,1)/Qnh0(pol,B,n,h,1):
end:

#ParSum(pol,n,n0): sum(1/pol(i),i=1..n0)
ParSum:=proc(pol,n,n0) local i:
convert([seq(1/subs(n=i,pol),i=1..n0)],`+`):
end:

#ImParSum(pol,n,n0): The improved partial sum
#e.g. ImParSum(n^3,n,10)
ImParSum:=proc(pol,n,n0) local C:
C:=FindAcc(pol,n):
if C=0 then
 RETURN(0):
fi:
ParSum(pol,n,n0)+subs(n=n0,C):
end:


#CheckC(C,n,n0): Given an expression C(n) and an integers n0
#computes
#sum((C(i+1)-C(i))/(C(i+1)*C(i)+C(i+1)+C(i)),i=n0+1..infinity)-
#1/(C(n0+1)). It should be small
#For example, try CheckC(2*n*(n-1),n,10);
CheckC:=proc(C,n,n0) 
local i:
evalf(sum( (subs(n=i+1,C)-subs(n=i,C))/
(subs(n=i+1,C)*subs(n=i,C)+subs(n=i+1,C)+subs(n=i,C)),
i=n0+1..infinity)-1/(1+subs(n=n0+1,C))):
end:


#RatImp(pol,n,x,h,h0): The rational functions that
#improves the convergence of the partial sum of
#sum(x^i/pol(i),i=1..n) to the infinite sum
#implied by the h0-th row of Apery's scheme
#for example, try RatImp(n^2,n,1,h,3); 
RatImp:=proc(pol,n,x,h,h0) local gu:
gu:=Aperyh(pol,x,n,h):
if gu=0 then
RETURN(0):
fi:
Rnh0(pol,gu[1],n,h,h0)/Qnh0(pol,gu[1],n,h,h0):
end:

#RatImpD1(pol,n,deg,L): Tries to find the rational function
#R(n), with monic denominator of of degree deg,
#R(n)-R(n-1)+1/pol(n) has numerator of degree L
#For example, try: RatImpD1(n^2,n,3,0);
RatImpD1:=proc(pol,n,d2,L)
local mone,mekh,a,b,eq,var,i,gu,R,d1,Deg1:
Deg1:=degree(numer(pol),n)-degree(denom(pol),n):
d1:=d2-(Deg1-1):
var:={}:
mone:=0:
for i from 0 to d1 do
 mone:=mone+a[i]*n^i:
 var:=var union {a[i]}:
od:

mekh:=n^d2:

for i from 0 to d2-1 do
 mekh:=mekh+b[i]*n^i:
 var:=var union {b[i]}:
od:

R:=mone/mekh:
gu:=expand(numer(normal(R-subs(n=n-1,R)+1/pol))):

eq:={}:

for i from L+1 to degree(gu,n) do
 eq:=eq union {coeff(gu,n,i)}:
od:

var:=solve(eq,var):

if var=NULL then
 RETURN(0):
fi:
if nops({var})>1 then
print(`Too many solutions`):
RETURN(0):
fi:
normal(subs(var,R)):
end:

#RatImpD(pol,n,deg): Tries to find the rational function
#R(n), with monic denominator of of degree deg,
#R(n)-R(n-1)+1/pol(n) has numerator of degree
#as small as possible
#For example, try: RatImpD(n^2,n,3);
RatImpD:=proc(pol,n,d2)
local L,gu:
option remember:
if d2<degree(numer(pol),n)-degree(denom(pol),n)-1 then
ERROR(`d2>=degree(pol,n)-1`):
fi:

gu:=0:
for L from 0 while gu=0 do
gu:=RatImpD1(pol,n,d2,L):
if gu<>0 then
RETURN(gu):
fi:
od:
end:


#Aluf(pol,Rat,n,L): the best appx. (according to delta)
#amongst sum(1/pol(i),i=1..n0)+Rat(n0) to 
#a:=sum(1/pol(i),i=1..infinity) for n0=2..L
#For example, try: Aluf(n^2,1/(n+1/2),n,5);
Aluf:=proc(pol,Rat,n,L) 
local shi, aluf,a,i,n0,shiu,ulai:
Digits:=100:
a:=sum(1/subs(n=i,pol),i=1..infinity):
a:=evalf(a):

aluf:=convert([seq(1/subs(n=i,pol),i=1..2)],`+`)+subs(n=2,Rat):
shi:=delt(aluf,a):
for n0 from 3 to L do
ulai:=convert([seq(1/subs(n=i,pol),i=1..n0)],`+`)+subs(n=n0,Rat):
shiu:=delt(ulai,a):
if shiu>shi then
  shi:=shiu:
  aluf:=ulai:
fi:
od:
aluf,evalf(shi,5):
end:

#Sidra(pol,n,L): The sequence of best appx. followed by the
#delta to sum(1/pol(i),i=1..infinity) furnished by the
#rational-function corrections supplied by 
#RatImpD(pol,n,d2), d2=degree(pol,n)-1..L
#For example, try: Sidra(n^2,n,10);

Sidra:=proc(pol,n,L)
local lu,d2,mu,ku:
lu:=[]:
for d2 from degree(numer(pol),n)-degree(denom(pol),n)-1 to L do
mu:=RatImpD(pol,n,d2):
ku:=[Aluf(pol,mu,n,3*d2)]:
lu:=[op(lu),ku]:
od:
lu:
end:




#RatImpCfD1(pol,Cf,n,deg,L): 
#Given a rational function pol(n), and a closed form sequence
#Cf(n), tries to find the rational function
#R(n), with monic denominator of of degree deg,
#(Cf(n)*R(n)-Cf(n-1)R(n-1)-Cf(n)/pol(n))/Cf(n) has numerator of 
#degree L; For example, try: RatImpCfD1(n^2,1,n,3,0);
RatImpCfD1:=proc(pol,Cf,n,d2,L)
local mone,mekh,a,b,eq,var,i,gu,R,d1,Deg1,ku:
Deg1:=degree(numer(pol),n)-degree(denom(pol),n):
d1:=d2-(Deg1-1):
var:={}:
mone:=0:
for i from 0 to d1 do
 mone:=mone+a[i]*n^i:
 var:=var union {a[i]}:
od:

mekh:=n^d2:

for i from 0 to d2-1 do
 mekh:=mekh+b[i]*n^i:
 var:=var union {b[i]}:
od:

R:=mone/mekh:
ku:=normal(expand(normal(simplify(subs(n=n-1,Cf)/Cf)))):
gu:=expand(numer(normal(R-ku*subs(n=n-1,R)+1/pol))):

eq:={}:

for i from L+1 to degree(gu,n) do
 eq:=eq union {coeff(gu,n,i)}:
od:

var:=solve(eq,var):

if var=NULL then
 RETURN(0):
fi:
if nops({var})>1 then
print(`Too many solutions`):
RETURN(0):
fi:
normal(subs(var,R)):
end:


#RatImpCfD(pol,Cf,n,deg): Given a rational function pol(n)
#and a closed-form function pol(n)
#Tries to find the rational function
#R(n), with monic denominator of of degree deg,
#(Cf(n)*R(n)-Cf(n-1)*R(n-1)-Cf(n)/pol(n))/Cf(n) 
#has numerator of degree
#as small as possible
#For example, try: RatImpCfD(n^2,1,n,3);
RatImpCfD:=proc(pol,Cf,n,d2)
local L,gu:
if d2<degree(numer(pol),n)-degree(denom(pol),n)-1 then
ERROR(`d2>=degree(pol,n)-1`):
fi:

gu:=0:
for L from 0 while gu=0 do
gu:=RatImpCfD1(pol,Cf,n,d2,L):
if gu<>0 then
RETURN(gu):
fi:
od:
end:



#AlufCf(pol,Cf,Rat,n,L): the best appx. (according to delta)
#amongst sum(Cf(i)/pol(i),i=1..n0)+Rat(n0) to 
#a:=sum(1/pol(i),i=1..infinity) for n0=2..L
#For example, try: AlufCf(n^2,1,1/(n+1/2),n,5);
AlufCf:=proc(pol,Cf,Rat,n,L) 
local shi, aluf,a,i,n0,shiu,ulai:
Digits:=100:
a:=sum(subs(n=i,Cf)/subs(n=i,pol),i=1..infinity):
a:=evalf(a):

aluf:=convert([seq(subs(n=i,Cf)/subs(n=i,pol),i=1..2)],`+`)+subs(n=2,Cf*Rat):
shi:=delt(aluf,a):
for n0 from 3 to L do
ulai:=convert([seq(subs(n=i,Cf)/subs(n=i,pol),i=1..n0)],`+`)+subs(n=n0,Cf*Rat):
shiu:=delt(ulai,a):
if shiu>shi then
  shi:=shiu:
  aluf:=ulai:
fi:
od:
aluf,evalf(shi,5):
end:

#SidraCf(pol,Cf,n,L): The sequence of best appx. followed by the
#delta to sum(Cf(i)/pol(i),i=1..infinity) furnished by the
#rational-function corrections supplied by 
#RatImpCfD(pol,Cf,n,d2), d2=degree(pol,n)-1..L
#For example, try: Sidra(n^2,1,n,10);

SidraCf:=proc(pol,Cf,n,L)
local lu,d2,mu,ku:
lu:=[]:
for d2 from degree(numer(pol),n)-degree(denom(pol)) to L do
mu:=RatImpCfD(pol,Cf,n,d2):
ku:=[AlufCf(pol,Cf,mu,n,3*d2)]:
lu:=[op(lu),ku]:
od:
lu:
end:


#FindaB(Q1,Q,pol,n): Given polynomials Q1(n),Q(n),
#pol(n) tries to find a constant a and a polynomial
#B(n) such that a*Q1(n)-Q(n+1)*pol(n+1)-Q(n)*B(n)=0
#If unsuccesful, returns 0
#For example, try: FindaB(2*n^3+3*n^2+3*n+1,n,n^2,n);
FindaB:=proc(Q1,Q,pol,n)
local a,B,gu,i,eq,var,c:
B:=0:
var:={a}:
for i from 0 to degree(pol,n) do
 B:=B+c[i]*n^i:
 var:=var union {c[i]}:
od:
gu:=expand(a*Q1-subs(n=n+1,Q)*subs(n=n+1,pol)-Q*B):
eq:={}:

for i from 0 to degree(gu,n) do
 eq:=eq union {coeff(gu,n,i)}:
od:

var:=solve(eq,var):

if var=NULL then
 RETURN(0):
fi:
a:=subs(var,a):
B:=subs(var,B):
a,B:
end:


#FindOpePol(Q,pol,n,N): tries to find an operator of the
#form ope(N,n)=pol(n+1)*pol(n+2)*N^2+pol(n+1)*P1(n)*N+P2(n)
#annihilating the polynomial Q(n)
#For example, try FindOpePol(1,n,n,N);
FindOpePol:=proc(Q,pol,n,N)
local ope,P1,P2,eq,var,i,c,d,gu:
var:={}:
P1:=0:
for i from 0 to degree(pol,n) do
 P1:=P1+c[i]*n^i:
 var:=var union {c[i]}:
od:

P2:=0:
for i from 0 to 2*degree(pol,n) do
 P2:=P2+d[i]*n^i:
 var:=var union {d[i]}:
od:

gu:=subs(n=n+1,pol)*subs(n=n+2,pol)*subs(n=n+2,Q)+
subs(n=n+1,pol)*P1*subs(n=n+1,Q)+P2*Q;
ope:=subs(n=n+1,pol)*subs(n=n+2,pol)*N^2+subs(n=n+1,pol)*P1*N+P2;

gu:=expand(gu):
eq:={}:

for i from 0 to degree(gu,n) do
eq:=eq union {coeff(gu,n,i)}:
od:

var:=solve(eq,var):

if var=NULL then
RETURN(0):
fi:

yafe(subs(var,ope),N):
end:




#FindOpe(Q,d,n,N): tries to find a second-order operator of the
#form ope(N,n)=P2*N^2+P1(n)*N+P0(n) where P0,P1,P2 are
#polynomials of degree d
#annihilating the polynomial Q(n)
#For example, try FindOpe(n,0,n,N);
FindOpe:=proc(Q,d,n,N)
local ope,P0,P1,P2,eq,var,i,a,b,c,gu:
var:={}:
P0:=0:P1:=0:P2:=0:
for i from 0 to d do
 P0:=P0+a[i]*n^i:
 P1:=P1+b[i]*n^i:
 P2:=P2+c[i]*n^i:
 var:=var union {a[i],b[i],c[i]}:
od:

gu:=P2*subs(n=n+2,Q)+P1*subs(n=n+1,Q)+P0*Q;
ope:=P2*N^2+P1*N+P0;

gu:=expand(gu):
eq:={}:

for i from 0 to degree(gu,n) do
eq:=eq union {coeff(gu,n,i)}:
od:

var:=solve(eq,var):

if var=NULL then
RETURN(0):
fi:

ope:=subs(var,ope):
if ope<>0 then
 RETURN(yafe(ope,N)):
 else
 RETURN(0):
fi:
end:




#AperyOp(pol,x,n,R,N):Given a polynomial pol in the variable n
#with polynomial coefficient, uses Apery's ORIGINAL method
#for TRYING to prove irrationality for the constant
#sum(x^i/pol(i),i=1..infinity)
#as described in his beautiful C.T.H.S. paper. It continues
#to accelerate up to R rows, if possible, and returns
#the list of accelerating B's, followed by the list of operators
#for the recurrences for the
#rows, and the expression themselves of the
#"easy" solution in row i divided by pol(1)pol(2)...pol(n)
#If it fails, it returns 0
#For example, try Apery(n^2,1,n,5);, Apery(n^3,1,n,5);
AperyOp:=proc(pol,x,n,R,N)
local ListB,ListP,ListQ,ListSol,P,Q,B,ope,Sol1,lu,Sol2,Sol2a,i1,ope1,i,
ListOpe:
ListB:=[]:
ListSol:=[1]:
Sol1:=1:
P:=-x*pol^2:
Q:=x*pol+subs(n=n+1,pol):
ope:=expand(N^2-subs(n=n+1,Q)*N-subs(n=n+1,P)):
ListP:=[]:
ListQ:=[]:
ListOpe:=[]:
for i from 1 to R do
ListP:=[op(ListP),P]:
ListQ:=[op(ListQ),Q]:
ListOpe:=[op(ListOpe),ope]:
lu:=FindB1(P,Q,n,0):

if nops(lu)<=1 then
# print(`Nothing found`):
  RETURN(0):
fi:

B:=lu[1]:

Sol2:=expand(subs(n=n+1,Sol1)*subs(n=n+1,pol)+Sol1*subs(n=n+1,B)):

for i1 from 2 to nops(lu) do
Sol2a:=expand(subs(n=n+1,Sol1)*subs(n=n+1,pol)+Sol1*subs(n=n+1,lu[i1])):
 if degree(Sol2a,n)>degree(Sol2,n) then
   Sol2:=Sol2a:
   B:=lu[i1]:
 fi:
od:
ListB:=[op(ListB),B]:
if nops(ListB)>1 and
degree(ListB[nops(ListB)],n)>degree(ListB[nops(ListB)-1],n) then
print(`ListB is`,ListB):
print(degree(ListB[nops(ListB)],n),degree(ListB[nops(ListB)-1],n)):
 RETURN(0):
fi:

Sol1:=Sol2:
ListSol:=[op(ListSol),Sol1]:
ope1:=expand(ResultOpe1(ope,N+subs(n=n+1,B),n,N)):

if ope1=0 then
RETURN(0):
fi:

P:=expand(subs(n=n-1,-coeff(ope1,N,0))):
Q:=expand(subs(n=n-1,-coeff(ope1,N,1))):
ope:=ope1:
od:

expand(ListB),expand(ListOpe):

end:



#Aperyh1Op(pol,x,n,R,h,N):Given a polynomial pol in the variable n
#with polynomial coefficient, uses Apery's ORIGINAL method
#and a variable h tries, if possible, to fit the
#accelerating B(n)'s and the coeffs. of the recurrence
#for row h by polynomials of degree h
#by fitting the data from Apery(pol,n,R)
#It returns B(n,h),and [P(n,h),Q(n.h)] if successful
#otherwise it returns 0. For example, try
#Aperyh1Op(n^2,1,n,8,h,N); and Aperyh1Op(n^3,1,n,8,h,N);
Aperyh1Op:=proc(pol,x,n,R,h,N)
local lu,B,Ope,ope1:
lu:=AperyOp(pol,x,n,R,N):

if lu=0 then
 RETURN(0):
fi:

B:=GuessPol(lu[1],h):
Ope:=GuessPol(lu[2],h):

if B=0 or Ope=0 then
RETURN(0):
fi:

ope1:=ResultOpe1(Ope,N+subs(n=n+1,B),n,N):
if not type(normal(ope1/subs(h=h+1,Ope)),integer) then
 print(`The conjectured B,Ope did not make it`,ope1,subs(h=h+1,Ope)):
 RETURN(0):
fi:
B,yafe(Ope,N):
end:


#AperyhOp(pol,x,n,h,N):Given a polynomial pol in the variable n
#with polynomial coefficient, uses Apery's ORIGINAL method
#and a variable h tries, if possible, to fit the
#accelerating B(n)'s and the coeffs. of the recurrence
#for row h by polynomials of degree h
#by fitting the data from Apery(pol,n,R)
#It returns B(n,h),and [P(n,h),Q(n.h)] if successful
#otherwise it returns 0. For example, try
#AperyhOp(n^2,1,n,h,n); and AperyhOp(n^3,1,n,h,N);
AperyhOp:=proc(pol,x,n,h,N)
local R,lu,lu1:
option remember:
lu:=AperyOp(pol,x,n,6,N):

if lu=0 then
print(`Nothing Found`):
 RETURN(0):
fi:

for R from 4 to 10 while Aperyh1Op(pol,x,n,R,h,N)=0 do
od:
lu:=Aperyh1Op(pol,x,n,R+1,h,N):
lu1:=Aperyh1Op(pol,x,n,R+2,h,N):

if lu<>lu1 then
print(`Nothing found`):
RETURN(0):
fi:
lu:
end:

#BestAcc(ope,IniT,IniB,n,N,ListDegs): Given an operator ope(n,N).
#and a list of degrees ListDegs tries to find
#the best accelerating operator ope1=P0(n)+P1(n)*N+P2(n)*N^2+...
#with ListDegs=[degree(P0(n),n),degree(P1(n),n),..]
#For example, try: BestAcc(IniOpe(n^2,1,n,N),n,N,[0,2]);
BestAcc:=proc(ope,IniT,IniB,n,N,ListDegs)
local L,gu:
for L from 0 to degree(ope,n)+2 do
gu:=BestAcc1(ope,IniT,IniB,n,N,ListDegs,L):
if gu<>0 then
 RETURN(gu):
fi:
od:
gu:
end:



#BestAcc1(ope,IniT,IniB,n,N,ListDegs,L): Given an operator ope(n,N).
#and initial values IniT, IniB (for top and bottom)
#and a list of degrees ListDegs tries to find
#the best accelerating operator (to within L)
#ope1=P0(n)+P1(n)*N+P2(n)*N^2+...
#with ListDegs=[degree(P0(n),n),degree(P1(n),n),..]
#For example, try: BestAcc(IniOpe(n^2,1,n,N).n,N,[0,2],0);
BestAcc1:=proc(ope,IniT,IniB,n,N,ListDegs,L)
local ope1,i,c,var,eq,j,gu,mua,lu,lua,aluf,shi,opemu,Ropemu:
_EnvAllSolutions:=true:
ope1:=0:
var:={}:
for i from 1 to nops(ListDegs) do
for j from 0 to ListDegs[i] do
ope1:=ope1+c[i-1,j]*n^j*N^(i-1):
var:=var union {c[i-1,j]}:
od:
od:
ope1:=subs(c[nops(ListDegs)-1,ListDegs[nops(ListDegs)]]=1,ope1):
var:=var minus {c[nops(ListDegs)-1,ListDegs[nops(ListDegs)]]}:
gu:=ResultOpe(ope,IniT,IniB,ope1,n,N)[1]:
gu:=coeff(gu,N,2):
eq:={}:
for i from L+1 to degree(gu,n) do
eq:=eq union {coeff(gu,n,i)}:
od:

var:=[solve(eq,var)]:

lu:=SeqQuotFromRec(ope,IniT,IniB,n,N,20):
aluf:=0:
shi:=abs(lu[21]-lu[20]):
for i from 1 to nops(var) do
if IsNumericSol(var[i]) then
 opemu:=subs(var[i],ope1):
 Ropemu:=ResultOpe(ope,IniT,IniB,opemu,n,N):
 lua:=SeqQuotFromRec(Ropemu,n,N,20):
 mua:=abs(lua[21]-lua[20]):
  if mua<shi then
   shi:=mua:
   aluf:=opemu:
  fi:
fi:
od:
aluf:
end:


#IsNumericSol(Sol): Given a solution set, Sol,
#returns true iff it is numeric;
#For example, try: IsNumericSol({a=a}); IsNumericSol({a=5});
IsNumericSol:=proc(Sol)
local i:
for i from 1 to nops(Sol) do
 if not type(op(2,Sol[i]),numeric) then
  RETURN(false):
 fi:
od:
true:
end:


#FindOpeNew(Q,d,n,N,pol): tries to find a second-order operator of the
#form ope(N,n)=pol(n+2)*P2*N^2+P1(n)*N+pol(n+1)*P0(n) where P0,P1,P2 are
#polynomials of degree d
#annihilating the polynomial Q(n)
#For example, try FindOpe(n,0,n,N);
FindOpeNew:=proc(Q,d,n,N,pol)
local ope,P0,P1,P2,eq,var,i,a,b,c,gu,i1,j1,ope1:
#option remember:
var:={}:
P0:=0:P1:=0:P2:=0:
for i from 0 to d-degree(pol,n) do
 P0:=P0+a[i]*n^i:
 var:=var union {a[i]}:
od:

for i from 0 to d do
 P1:=P1+b[i]*n^i:
 var:=var union {b[i]}:
od:

for i from 0 to d-degree(pol,n) do
 P2:=P2+c[i]*n^i:
 var:=var union {c[i]}:
od:



gu:=P2*subs(n=n+2,Q*pol)+P1*subs(n=n+1,Q)+P0*Q*subs(n=n+1,pol);
ope:=P2*subs(n=n+2,pol)*N^2+P1*N+P0*subs(n=n+1,pol);

gu:=expand(gu):
eq:={}:

for i from 0 to degree(gu,n) do
eq:=eq union {coeff(gu,n,i)}:
od:

var:=solve(eq,var):

if var=NULL then
RETURN(0):
fi:

ope:=expand(subs(var,ope)):



if ope<>0 then
 ope1:=yafe(ope,N):
 ope:=expand(ope1):

for i1 from 0 to degree(ope,n) do
for j1 from 0 to degree(coeff(ope,n,i1),N) do
 if not type(coeff(coeff(ope,n,i1),N,j1),numeric) then
    RETURN(0):
 fi:
od:
od:

 RETURN(ope1):
 else
 RETURN(0):
fi:
end:

#AperyhOper(pol,x,n,h,N):Like Aperyh(pol,x,n,h,N) but
#returns the accelerating operator followed by 
#the recurrence satistied in the n-direction
#For example, try: AperyhOper(n^2,1,n,h,N);
AperyhOper:=proc(pol,x,n,h,N)
local gu,B,P,Q:
gu:=Aperyh(pol,x,n,h):
B:=gu[1]:P:=gu[2]:Q:=gu[3]:
N+expand(subs(n=n+1,B)),
CODVpol(N^2-expand(subs(n=n+1,Q))*N-expand(subs(n=n+1,P)),n,N,pol):
end:

#AperyNes(pol,x,n,N,h) : Given a polynomial pol, in n
#finds, if possible, an Apery-style miracle pair
#accelarating operator: Acc0=N+B(n+1,h+1) and a row-operator
#ope:=N^2+P(n,h)*N+pol(n+1)^2 such that
#ResultOpe(ope,[1,0],[0,1],Acc,n,N) is ope with
#h replaced by h+1, and B is  poly of (n,h) of total degree
#degree(pol,n), and P is a (not necessarily homog.) 
#polynomial of degree degree(pol,n)
#For example, try:AperyNes(n^2,1,n,N,h) 
AperyNes:=proc(pol,x,n,N,h) 
local Acc,ope,eq,var,P,B,p,b,i,j,k,deg,ope0,ope1,gu,gu1,gu2,gu3:
_EnvAllSolutions:=true:
deg:=degree(pol,n):
var:={}:
B:=0:
for i from 0 to deg do
for j from 0 to deg-i do
var:=var union {b[i,j]}:
B:=B+b[i,j]*n^i*h^j:
var:=var union {b[i,j]}:
od:
od:

ope0:=IniOpe(pol,x,n,N):
P:=coeff(ope0[1],N,1):
for i from 0 to deg do
for j from 1 to deg-i do
var:=var union {p[i,j]}:
P:=P+p[i,j]*n^i*h^j:
od:
od:

ope:=N^2+P*N+x*subs(n=n+1,pol)^2:
Acc:=N+B:
ope1:=ResultOpe(ope,[0,1],[1,0],Acc,n,N)[1]:

eq:={}:
gu:=expand(ope1-coeff(ope1,N,2)*subs(h=h+1,ope)):
for i from 0 to degree(gu,N) do
 gu1:=coeff(gu,N,i):
 for j from 0 to degree(gu1,n) do
  gu2:=coeff(gu1,n,j):
  for k from 0 to degree(gu2,h) do
   gu3:=coeff(gu2,h,k): 
     eq:=eq union {gu3}:
  od:
 od:
od:

var:=[solve(eq,var)]:
gu:=[]:
for i from 1 to nops(var) do
if expand(diff(subs(var[i],ope),h))<>0 
and ResultOpe(ope0,subs(h=0,subs(var[i],Acc)),n,N)<>0 then

gu:=[op(gu),[subs(var[i],Acc),subs(var[i],ope)]]:
fi:
od:

if gu=[] then
 RETURN(0):
fi:
if nops(gu)>1 then
 RETURN(gu):
fi:

op(gu[1]):
end:


#MakePol(sid,x,y,deg): Given a sequence sid, of polynomials
#in x,y tries to find a sequence a[h], of the same
#length of sid, such that sid[h]*a[h] is a polynomial of degree
#deg in h
#For example, try 
#MakePol([(1+n*m+n^2)/2,(1+n*m+n^2)*3,(1+n*m+n^2)*4],n,m,1);
MakePol:=proc(sid,x,y,deg)
local a,eq,sid1,i,lu,i1,mu,mu1,j,k,var:
if nops(sid)-deg<2 then
ERROR(`sequence is too small`):
fi:

lu:=[seq(a[i],i=1..nops(sid))]:
var:=convert(lu,set):
sid1:=[seq(sid[i]*a[i],i=1..nops(sid))]:
for i from 1 to deg+1 do
sid1:=expand([seq(sid1[i1+1]-sid1[i1],i1=1..nops(sid1)-1)]):
od:

eq:={}:
for i from 1 to nops(sid1) do
mu:=sid1[i]:
 for j from 0 to degree(mu,x) do
   mu1:=coeff(mu,x,j):
      for k from 0 to degree(mu,y) do
       eq:=eq union {coeff(mu1,y,k)}:
      od:
 od:
od:
var:=solve(eq,var):
subs(var,lu):
end:


#Halevay1(pol,n,h,N,start1,end1,res1,degn,degh): Given a polynomial
#pol, tries to guess second-order operators of degn
#in n and degh in h for the rows of the output of RatImpD
#starting at start1, ending at end1, and with resolution res1
#For example, try:
#Halevay1(n^2,n,h,N,2,6,1,2,2)
Halevay1:=proc(pol,n,h,N,start1,end1,res1,degn,degh)
local gu,n0,mu,i,lu:
if trunc((end1-start1)/res1)+1-degh<2 then
ERROR(`end1 is too small `):
fi:
gu:=[seq(
denom(RatImpD(pol,n,start1+res1*n0)),n0=0..trunc((end1-start1)/res1))]:

mu:=[seq(FindOpeNew(gu[i],degn,n,N,pol),i=1..nops(gu))]:

if mu[1]=0 then
 RETURN(0):
fi:

lu:=MakePol(mu,n,N,degh):

if lu[1]=0 then
RETURN(0):
fi:

mu:=[seq(mu[i]/lu[i],i=1..nops(mu))]:
lu:=GuessPolNew(mu,h,1,degh):
#mu:
lu:
end:


#Halevay(pol,n,h,N,start1,end1,res1,degn,Maxdegh): Given a polynomial
#pol, tries to guess second-order operators of degn
#in n and at most Maxdegh in h for the rows of the output of RatImpD
#starting at start1, ending at end1, and with resolution res1
#For example, try:
#Halevay(n^2,n,h,N,2,6,1,2,2)
Halevay:=proc(pol,n,h,N,start1,end1,res1,degn,Maxdegh)
local gu,i:

if trunc((end1-start1)/res1)+1-Maxdegh<2 then
ERROR(`end1 is too small `):
fi:

gu:=0:

for i from 0 to Maxdegh while gu=0 do
gu:=Halevay1(pol,n,h,N,start1,end1,res1,degn,i):
if gu<>0 then

 RETURN(yafe(gu,N)):
fi:
od:
0:
end:


#AppxDel(pol,x,n,L): Gives the Lth diag. approximant
#Apery's recurrence for x/pol(1)+x^2/pol(2)+x^3/pol(3)+... 
#followed by the error and the delta. For example: 
#try AppxDel(n^2,1,n,20);
AppxDel:=proc(pol,x,n,L) local gu,mu:

gu:=Appx(pol,x,n,L):

mu:=sum(x^n/pol,n=1..infinity):

[gu[1],gu[2],delt(gu[1],mu)]: 

end: