######################################################################
## AbelCeline.txt Save this file as  AbelCeline.txt  to use it,      #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read  `AbelCeline.txt` <Enter>                    #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################


print(`First Written: April 2, 2024: tested for Maple 2020 `):
print(`Version  : April 2, 2024 `):
print():
print(`This is AbelCeline.txt, one of the Maple packages`):
print(`accompanying Gil Kalai and Doron Zeilberger's  article: `):
print(` Bijective and Automated Approaches to Abel Sums `):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/mamarim/mamarimhtml/abelKZ.html `):
print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/AbelCeline.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`--------------------------`):
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`--------------------------`):
print(`For a list of the supporting functions type: ezra1();`):
print(`For specific help type "ezra(procedure_name);" `):
print(`--------------------------`):

print(`For a list of the Checking functions type: ezraC();`):
print(`For specific help type "ezra(procedure_name);" `):
print(`--------------------------`):

print(`For a list of the Story functions type: ezraS();`):
print(`For specific help type "ezra(procedure_name);" `):
print(`--------------------------`):

print():





ezraC:=proc()
if args=NULL then

print(`The Checking procedures are: CheckDiffPair, CheckOpe, CheckOpeSeq `):
print(``):

else
ezra(args):

fi:


end:

ezraS:=proc()
if args=NULL then

print(`The Story procedures are: Paper, PaperD, PaperWP, PrintOpe, PrintOpeD `):
print(``):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are:  ApplyOpeD, ApplyOpeSeq, ApplyOpeSeqD, FindOpe1,  FindOpeD1g, PolA, PolAg `):
print(``):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` AbelCeline.txt: A Maple package for  automatically generating functional equations for Abel sums`):
print(`It is is based on extensions, and differntial analogs, of John Majewicz's PhD from 1996 (under the direction of Doron Zeilberger)`):

print(`The MAIN procedures are: DiffPair, FindAbelRec, FindOpe, FindOpeDg, NiceOpe, SeqA, SeqAg, SeqAf `):


print(``):


elif nargs=1 and args[1]=ApplyOpeSeq then
print(`ApplyOpeSeq(ope,N,R,S,n,r,s,L): applies the operator ope(N,R,S,n,r,s) to the sequence L. Try:`):
print(`F:=n!/(k!*(n-k)!): ope:=FindOpe(F,n,k,r,s,R,S,N,K,2): ope:=subs(K=1,ope): L:=SeqA(F,n,k,r,s,2,3,10): ApplyOpeSeq(ope,N,R,S,n,r,s,L);`):

elif nargs=1 and args[1]=ApplyOpeSeqD then
print(`ApplyOpeSeqD(ope,N,Dr,Ds,n,r,s,L): applies the operator ope(N,R,S,n,r,s) to the sequence L. Try: `):


elif nargs=1 and args[1]=CheckOpe then
print(`CheckOpe(ope,n,k,N,K,r,s,R,S,F,p,q): checks the operator ope on F`):
print(`Try: `):
print(` F:=n!/(k!*(n-k)!): ope:=FindOpe(F,n,k,r,s,R,S,N,K,2): CheckOpe(ope,n,k,N,K,r,s,R,S,F,2,3);`):


elif nargs=1 and args[1]=ApplyOpeD then
print(`ApplyOpeD(ope,N,Dr,Ds,n,F): applies the operator ope(N,Dr,Ds) to F(n,r,s) divided by F(n,r,s). Try:`):
print(`ApplyOpeD(N+Dr+Ds,N,Dr,Ds,n,r,s,n*r*s);`):


elif nargs=1 and args[1]=CheckDiffPair then
print(`CheckDiffPair(F,n,k,r,s,N,Dr,Ds,MaxOrd,KER,KAMA): checks DiffPair(F,n,k,r,s,N,Dr,Ds,MaxOrd,KER) up to KAMA terms Try:`):
print(`CheckDiffPair(binomial(n,k),n,k,r,s,N,Dr,Ds,2,(r+k)^(k-1)*(s-k)^(n-k),20);`):

elif nargs=1 and args[1]=CheckOpeSeq then
print(`CheckOpeSeq(F,n,k,ORD,p,q,KAMA): checks that the operator induced by F does what it is supposed to do. Try:`):
print(`CheckOpeSeq(n!/k!/(n-k)!,n,k,2,0,0,20);`):


elif nargs=1 and args[1]=DiffPair then
print(`DiffPair(F,n,k,r,s,N,Dr,Ds,MaxOrd,KER) :`):
print(` inputs a hypergeometric term F in n and k, also a kernel KER, and outputs a k-free operators  of (n,r,N,Dr) and (n,s,N,Ds) respectively`):
print(`where sum of recurrence order and the differntial order is <=MaxOrd`):
print(`annihilating F*KER`):
print(`If none exists, then it returns FAIL. Try:`):
print(`DiffPair(binomial(n,k),n,k,r,s,N,Dr,Ds,2, (r+k)^(k-1+p)*(s-k)^(n-k+q));`):

elif nargs=1 and args[1]=FindAbelRec then
print(`FindAbelRec(M,n,r,s,N,R,S,x,MaxOrd): finds all hypergemetric terms of the from x^k/(k!^i*(n-k)^j with 1<=i,j<=M that yields functional recurrences of order <=MaxOrd. `):
print(`it also returns the operators annihilating the Abel-sums induced by F. Try:`):
print(`FindAbelRec(2,n,r,s,N,R,S,x,4);`):

elif nargs=1 and args[1]=FindOpe then
print(`FindOpe(F,n,k,r,s,R,S,N,K,MaxOrd): inputs a hypergeometric term F in n and k, and outputs a k-free operator, of order <=MaxOrd,  annihilating`):
print(`F*(r+k)^(k-1+p)*(s-k)^(n-k+q), for any p or q. If it fails, it retunrs FAIL. Try:`):
print(`FindOpe(n!/k!/(n-k)!,n,k,r,s,R,S,N,K,2); `):

elif nargs=1 and args[1]=FindOpe1 then
print(`FindOpe1(F,r,s,n,k,R,S,N,K,ORDn,ORDk): inputs a hypergeometric term F in n and k, and outputs a k-free operator annihilating`):
print(`F*(r+k)^(k-1+p)*(s-k)^(n-k+q), for any p or q. If it exits. Otherwise it returns FAIL. Try:`):
print(`FindOpe1(binomial(n,k),n,k,r,s,R,S,N,K,2,1)`):


elif nargs=1 and args[1]=FindOpeD1g then
print(`FindOpeD1g(F,n,k,r,N,Dr,ORDn,ORDr,KER) :`):
print(` inputs a hypergeometric term F in n and k, also a kernel KER, and outputs a k-free operator  of (n,r,N,Dr)`):
print(`of order ORDn in n (i.e. degree in N), order ORDr in r (i.e. degree in Dr))`):
print(`annihilating F*KER`. Try):
print(`FindOpeD1g(binomial(n,k),n,k,r,N,Dr,1,1, (r+k)^(k-1+p)*(s-k)^(n-k+q));`):


elif nargs=1 and args[1]=FindOpeDg then
print(`FindOpeDg(F,n,k,r,N,Dr,MaxOrd,KER) :`):
print(` inputs a hypergeometric term F in n and k, also a kernel KER, and outputs a k-free operator  of (n,r,N,Dr)`):
print(`sum of theorder ORDn in n (i.e. degree in N) and the order ORDr in r (i.e. degree in Dr))<=MaxOrd, or returns FAIL`):
print(`annihilating F*KER. Try:`):
print(`FindOpeDg(binomial(n,k),n,k,r,N,Dr,2, (r+k)^(k-1+p)*(s-k)^(n-k+q));`):

elif nargs=1 and args[1]=NiceOpe then
print(`NiceOpe(F,n,k,r,s,R,S,N,MaxOrd): an operator of form Oper(N^(-1),n,R,S) such that the sequence of expressions in r,s`):
print(`a(n,r,s):=Sum(F*(r+k)^(k-1+p)*(s-k)^(n-k+q),k=0..n) satisfies a(n,r,s)=Oper(N^(-1),n,R,S)a(n,r,s) enabling a fast`):
print(`computation of a(n,r,s). Returns FAIL if it does not exist. Try:`):
print(`NiceOpe(n!/k!/(n-k)!,n,k,r,s,R,S,N,2);`):

elif nargs=1 and args[1]=Paper then
print(`Paper(F,n,k,r,s,MaxOrd): inputs a hypergeometric term F in discrete variables n and k and symbols. outputs a paper with just the statement of the functional recurrence equation.`):
print(`satisfied by`):
print(`A[n](r,s):=Sum(F*(r+k)^(k-1+p)*(s-k)^(n-k+q), k=0..n); for any integers p and q. Try: `):
print(`Paper(binomial(n,k)*2^k,n,k,r,s,2);`):


elif nargs=1 and args[1]=PaperD then
print(`PaperD(F,n,k,r,s,MaxOrd,KER): inputs a hypergeometric term F in discrete variables n and k and symbols, also a kernel, an expression in (r,s,n,k). outputs a paper with just the statement of the differential recurrence equations, `):
print(` order of recurrence plus order of differntiations <=MaxOrder one with respect to r and s`):
print(`satisfied by`):
print(`A[n](r,s):=Sum(F*KER, k=0..n); for any integers p and q, and number x. Try: `):
print(`PaperD(binomial(n,k),n,k,r,s,2,(r+k)^(k-1+p)*(s-k)^(n-k+q)*x^k);`):


elif nargs=1 and args[1]=PaperWP then
print(`PaperWP(F,n,k,r,s,MaxOrd): inputs a hypergeometric term F in discrete variables n and k and symbols. outputs a paper with  the statement AND proof, of the functional recurrence equation.`):
print(`satisfied by`):
print(`A[n](r,s):=Sum(F*(r+k)^(k-1+p)*(s-k)^(n-k+q), k=0..n); for any integers p and q. Try: `):
print(`PaperWP(binomial(n,k)*2^k,n,k,r,s,2);`):


elif nargs=1 and args[1]=PolA then
print(`PolA(F,n,k,r,s,n1,p,q): Given a bi-variate hypergemetric term F in n,k, outputs`):
print(`Sum(F(n1,k)*(r+k)^(k-1+p)*(s-k)^(n1-k+q),k=0..n1); Try:`):
print(`PolA(binomial(n,k),n,k,r,s,0,0,5);`):


elif nargs=1 and args[1]=PolAg then
print(`PolAg(F,n,k,KER,n1): Given a bi-variate hypergemetric term F in n,k, and a kernel KER=KER(n,k), depending on n,k,outputs`):
print(`Sum(F(n1,k)*KER(n,k),k=0..n1); Try:`):
print(`PolAg(binomial(n,k),n,k,(r+k)^(k-1)*(s-k)^(n-k),5);`):


elif nargs=1 and args[1]=PrintOpe then
print(`PrintOpe(ope,n,r,s,N,R,S,A): inputs a difference operator in n,r,s, with shift operators N,R,S, applied to A[n](r,s): outputs it in plain Engllsh. Try:`):
print(`PrintOpe(r*N*R*S+s*R+r*S,n,r,s,N,R,S,A);`):

elif nargs=1 and args[1]=PrintOpeD then
print(`PrintOpeD(ope,n,r,s,N,Dr,Ds,A): inputs a recurrence (in n) -and (differential) operator in,r and s, with shift operator N and differnetiation Dr, Ds, respectively applied to A[n](r,s): outputs it in plain Engllsh. Try:`):
print(`PrintOpeD(r*N*Dr+s*Ds,n,r,s,N,Dr,Ds,A);`):


elif nargs=1 and args[1]=SeqA then
print(`SeqA(F,n,k,r,s,p,q,N1): Given a bi-variate hypergemetric term F in n,k, outputs`):
print(`Sum(F(n1,k)*(r+k)^(k-1+p)*(s-k)^(n1-k+q),k=1..n1); Try:`):
print(`SeqA(binomial(n,k),n,k,r,s,0,0,10);`):


elif nargs=1 and args[1]=SeqAf then
print(`SeqAf(F,n,k,r,s,p,q,MaxOrd,N1): Same output as SeqA(F,n,k,r,s,p,q,N1)  (q.v.) but done fast,`):
print(`using the explicit operator with maximum order MaxOrd.  If there is no operator it uses SeqA. Try:`):
print(`SeqAf(n!/(k!*(n-k)!),n,k,r,s,0,0,2,10);`):


elif nargs=1 and args[1]=SeqAg then
print(`SeqAg(F,n,k,KER,N1): Given a bi-variate hypergemetric term F in n,k, and a kernel KER=KER(n,k) outputs `):
print(`[seq(Sum(F(n1,k)*KER(n1,k),k=0..n1),n1=0..N1); Try`):
print(`SeqAg(binomial(n,k),n,k,(r+k)^(k-1)*(s-k)^(n-k),10);`):

else
 print(`There is no such thing as`, args):

fi:


end:

ez:=proc(): 
print(`CheckOpeSeq(F,n,k,ORD,p,q,KAMA)`):
end:


#FindOpe1(F,r,s,n,k,R,S,N,K,ORDn,ORDk): inputs a hypergeometric term F in n and k, and outputs a k-free operator annihilating
#F*(r+k)^(k-1)*(s-k)^(n-k). Try:
#FindOpe1(binomial(n,k),n,k,r,s,R,S,N,K,2,1)
FindOpe1:=proc(F,n,k,r,s,R,S,N,K,ORDn,ORDk) local F1,i,j,a,eq,var,ope,EXPR,eq1,tov,var1:

F1:=convert(F,factorial):

ope:=add(add(a[i,j]*(S*K/R)^i*N^(j),i=0..ORDk),j=0..ORDn):
var:={seq(seq(a[i,j],i=0..ORDk),j=0..ORDn)}:

EXPR:=expand(numer(normal(add(add(a[i,j]*simplify(subs({n=n+j,k=k+i},F1)/F1)*(r+k)^i*(s-k)^(j-i),j=0..ORDn),i=0..ORDk)))):

eq:={seq(coeff(EXPR,k,i),i=0..degree(EXPR,k))}:

eq1:=subs({r=1/3,s=2/5},eq):

if subs(solve(eq1,var),ope)=0 then
RETURN(FAIL):
fi:


var:=solve(eq,var):


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

tov:={}:

for var1 in var do
 if op(1,var1)=op(2,var1) then
   tov:=tov union {op(1,var1)}:
 fi:
od:

if nops(tov)<>1 then
 RETURN(FAIL):
fi:


ope:=subs(var,ope):
ope:=subs(tov[1]=1,ope):

ope:=add(add(factor(coeff(coeff(ope,N,i),K,j))*N^i*K^j,i=0..degree(ope,N)),j=0..degree(ope,K)):
ope:
end:

#FindOpe(F,n,k,r,s,R,S,N,K,MaxOrd) : an opeator with order ORDk<=ORDn<=MaxOrd or FAIL. Try:
#FindOpe(binomial(n,k),n,k,r,s,R,S,N,K,2);
FindOpe:=proc(F,n,k,r,s,R,S,N,K,MaxOrd)  local gu,ORDn,ORDk:
option remember:
for ORDn from 1 to MaxOrd do
 for ORDk from 0 to ORDn do
  gu:=FindOpe1(F,n,k,r,s,R,S,N,K,ORDn,ORDk):
   if gu<>FAIL then
     RETURN(gu):
   fi:
 od:
od:
FAIL:
end:


#PolA(F,n,k,r,s,p,q,n1): Given a bi-variate hypergemetric term F in n,k, outputs
#Sum(F(n1,k)*(r+k)^(k-1+p)*(s-k)^(n1-k+q),k=0..n1); Try:
#PolA(binomial(n,k),n,r,s,5,0,0)
PolA:=proc(F,n,k,r,s,p,q,n1) local k1:
factor(normal(add(eval(subs({n=n1,k=k1},F))*(r+k1)^(k1-1+p)*(s-k1)^(n1-k1+q),k1=0..n1))):
end:

#SeqA(F,n,k,r,s,p,q,N1): Given a bi-variate hypergemetric term F in n,k, outputs 
#[seq(Sum(F(n1,k)*(r+k)^(k-1+p)*(s-k)^(n1-k+q),k=0..n1),n1=0..N1); Try
#SeqA(binomial(n,k),n,k,r,s,0,0,10);
SeqA:=proc(F,n,k,r,s,p,q,N1) local n1:
[seq(PolA(F,n,k,r,s,p,q,n1),n1=1..N1)]:
end:



#CheckOpe(ope,n,k,N,K,r,s,R,S,F,p,q): checks the operator ope on F
CheckOpe:=proc(ope,n,k,N,K,r,s,R,S,F,p,q) local F1,i1,i2,i3,i4,ope1,ope2,ope3,ope4,lu:
F1:=F*(r+k)^(k-1+p)*(s-k)^(n-k+q):

lu:=0:
for i1 from 0 to degree(ope,N) do
ope1:=coeff(ope,N,i1):
for i2 from 0 to degree(ope1,K) do
ope2:=coeff(ope1,K,i2):
for i3 from ldegree(ope2,R) to degree(ope2,R) do
ope3:=coeff(ope2,R,i3):
for i4 from ldegree(ope3,S) to degree(ope3,S) do
ope4:=coeff(ope3,S,i4):
lu:=normal(lu+ope4*normal(simplify(subs({n=n+i1,k=k+i2,r=r+i3,s=s+i4},F1)/F1))):
od:
od:
od:
od:
evalb(lu=0):
end:




#ApplyOpe(ope,s,r,R,S,f): given an expression f in r,s and an operator ope(r,s,R,S) applies it fo f
ApplyOpe:=proc(ope,r,s,R,S,f) local i1,i2:
normal(add(add(coeff(coeff(ope,R,i1),S,i2)*subs({r=r+i1,s=s+i2},f),i2=ldegree(ope,S)..degree(ope,S)),i1=ldegree(ope,R)..degree(ope,R))):
end:


#ApplyOpeSeq(ope,N,R,S,n,r,s,L): applies the operator ope(N,R,S,n,r,s) to the sequence L. 

ApplyOpeSeq:=proc(ope,N,R,S,n,r,s,L) local lu,lu1,n1,i1:
lu:=[]:
for n1 from 1 to nops(L)-degree(ope,N) do
 lu1:=normal(add(ApplyOpe(subs(n=n1,coeff(ope,N,i1)),r,s,R,S,L[n1+i1]),i1=0..degree(ope,N))):
 lu:=[op(lu),lu1]:
 od:
lu:
end:

#CheckOpeSeq(F,n,k,ORD,p,q,KAMA): checks that the operator induced by F does what it is supposed to do. Try:
#CheckOpeSeq(n!/k!/(n-k)!,n,k,2,0,0,20);
CheckOpeSeq:=proc(F,n,k,ORD,p,q,KAMA) local r,s,R,S,N,K,ope,L,L1,ope1:

ope:=FindOpe(F,n,k,r,s,R,S,N,K,ORD); 
if ope=FAIL then
 print(`No operator of order`, ORD):
 RETURN(FAIL):
fi:

ope1:=subs(K=1,ope);  

L:=SeqA(F,n,k,r,s,p,q,KAMA): 
L1:=ApplyOpeSeq(ope1,N,R,S,n,r,s,L):
evalb(L1=[0$nops(L1)]):
end:



#NiceOpe(F,n,k,r,s,R,S,N,MaxOrd): an operator of form Oper(N^(-1),n,R,S) such that the sequence of expressions in r,s
#a(n,r,s):=Sum(F*(r+k)^(k-1+p)*(s-k)^(n-k+q),k=0..n) satisfies a(n,r,s)=Oper(N^(-1),n,R,S)a(n,r,s) enabling a fast
#computation of a(n,r,s). Returns FAIL if it does not exist. Try:
#NiceOpe(n!/k!/(n-k)!,n,k,r,s,R,S,N,2);
NiceOpe:=proc(F,n,k,r,s,R,S,N,MaxOrd) local K,ope,ORD,ordR,ordS,coe0:
ope:=FindOpe(F,n,k,r,s,R,S,N,K,MaxOrd):

if ope=FAIL then
  RETURN(FAIL):
fi:

ope:=subs(K=1,ope):

ORD:=degree(ope,N):

if type(coeff(ope,N,ORD),`+`) then
 RETURN(FAIL):
fi:

ope:=expand(subs(n=n-ORD,ope)/N^ORD):

coe0:=coeff(ope,N,0):
ordR:=degree(coe0,R):
ordS:=degree(coe0,S):
coe0:=normal(coe0/(R^ordR*S^ordS)):

if normal(diff(coe0,R))<>0 or normal(diff(coe0,S))<>0  then
 RETURN(FAIL):
fi:

ope:=expand(ope/coe0):

ope:=subs({r=r-ordR,s=s-ordS},ope)/(R^ordR*S^ordS):
ope:=expand(ope):
ope:=1-ope:

collect(ope,N):

end:





#SeqAf(F,n,k,r,s,p,q,MaxOrd,N1): Same output as SeqA(F,n,k,r,s,p,q,N1)  (q.v.) but done fast,
#using the explicit operator with maximum order MaxOrd.  If there is no operator it uses SeqA. Try:
#SeqAf(binomial(n,k),n,k,r,s,0,0,2,10);
SeqAf:=proc(F,n,k,r,s,p,q,MaxOrd,N1) local n1,ope,gu,N,R,S,k1,i1,newg:

ope:=NiceOpe(F,n,k,r,s,R,S,N,MaxOrd):

if ope=FAIL then
 print(`No explicit operator, doing it the slow way`):
 RETURN(SeqA(F,n,k,r,s,p,q,N1)):
fi:

k1:=-ldegree(ope,N):


gu:=SeqA(F,n,k,r,s,p,q,k1):

for n1 from k1+1 to N1 do
 newg:=normal(add(ApplyOpe(subs(n=n1,coeff(ope,N,-i1)),r,s,R,S,gu[n1-i1]) ,i1=1..k1)):
 gu:=[op(gu),newg]:

od:

gu:
end:


#Paper(F,n,k,r,s,R,S,N,K,MaxOrd): inputs a hypergeometric term F in discrete variables n and k and symbols. Try:
#Paper(binomial(n,k)*2^k,n,k,r,s,R,S,N,K,2);
Paper:=proc(F,n,k,r,s,MaxOrd) local R,S,N,K,ope,ope1,p,q,ope1N,A,i1,j1,j2,t0:
t0:=time():

ope:=FindOpe(F,n,k,r,s,R,S,N,K,MaxOrd):

if ope=FAIL then
 print(`Try to make MaxOrd bigger`):
 RETURN(FAIL):
fi:

ope1:=subs(K=1,ope):

ope1N:=NiceOpe(F,n,k,r,s,R,S,N,MaxOrd):


if ope1N<>FAIL then
print(``):
 print(`Theorem: For any integers p and q define the Abel-sum type sequence  by`):
print(``):
 print(A[n](r,s)=Sum(F*(r+k)^(k-1+p)*(s-k)^(n-k+q),k=0..n)):
print(``):
print(`Then we have the following functional recurrence expressing`, A[n](r,s) , `in terms of `, seq(A[n-i1](r,s),i1=1..-ldegree(ope1N,N))):
print(``):
print(`:`):
print(``):
print(A[n](r,s)=add(add(add(normal(coeff(coeff(coeff(ope1N,N,-i1),R,j1),S,j2))* A[n-i1](r+j1,s+j2),j1=ldegree(ope1N,R)..degree(ope1N,R)),j2=ldegree(ope1N,S)..degree(ope1N,S)),i1=1..-ldegree(ope1N,N))):
print(``):
print(`and in Maple notation`):
print(``):
lprint(A[n](r,s)=add(add(add(normal(coeff(coeff(coeff(ope1N,N,-i1),R,j1),S,j2))* A[n-i1](r+j1,s+j2),j1=ldegree(ope1N,R)..degree(ope1N,R)),j2=ldegree(ope1N,S)..degree(ope1N,S)),i1=1..-ldegree(ope1N,N))):
print(``):
print(`-------------------------------------------------`):
print(``):
print(`This took`, time()-t0, `seconds. `):
print(``):

fi:


if ope1N=FAIL then
print(``):
 print(`Theorem: For any integers p and q define the Abel-sum type sequence  by`):
print(``):
 print(A[n](r,s)=Sum(F*(r+k)^(k-1+p)*(s-k)^(n-k+q),k=0..n)):
print(``):
print(`Then we have the following functional recurrence relating `, seq(A[n+i1](r,s),i1=0..degree(ope1,N))):

print(add(add(add(normal(coeff(coeff(coeff(ope1,N,i1),R,j1),S,j2))* A[n+i1](r+j1,s+j2),j1=ldegree(ope1 ,R)..degree(ope1 ,R)),j2=ldegree(ope1 ,S)..degree(ope1 ,S)),i1=0..degree(ope1 ,N))=0):
print(``):
print(`and in Maple notation`):
print(``):
lprint(add(add(add(normal(coeff(coeff(coeff(ope1,N,i1),R,j1),S,j2))* A[n+i1](r+j1,s+j2),j1=ldegree(ope1 ,R)..degree(ope1 ,R)),j2=ldegree(ope1 ,S)..degree(ope1 ,S)),i1=0..degree(ope1 ,N))=0):
fi:


end:



#PaperWP(F,n,k,r,s,R,S,N,K,MaxOrd): inputs a hypergeometric term F in discrete variables n and k and symbols. Try:
#PaperWP(binomial(n,k)*2^k,n,k,r,s,R,S,N,K,2);
PaperWP:=proc(F,n,k,r,s,MaxOrd) local F1,R,S,N,K,ope,ope1,p,q,ope1N,A,i1,j1,j2,i2,t0:
t0:=time():

ope:=FindOpe(F,n,k,r,s,R,S,N,K,MaxOrd):

if ope=FAIL then
 print(`Try to make MaxOrd bigger`):
 RETURN(FAIL):
fi:

ope1:=subs(K=1,ope):

ope1N:=NiceOpe(F,n,k,r,s,R,S,N,MaxOrd):


if ope1N<>FAIL then
 print(`Theorem: For any integers p and q define the Abel-sum type sequence  by`):
print(``):
 print(A[n](r,s)=Sum(F*(r+k)^(k-1+p)*(s-k)^(n-k+q),k=0..n)):
print(``):
print(`Then we have the following functional recurrence expressing A[n](r,s) in terms of `, seq(A[n-i1](r,s),i1=1..-ldegree(ope1N,N))):
print(``):
print(`:`):
print(``):

print(A[n](r,s)=add(add(add(normal(coeff(coeff(coeff(ope1N,N,-i1),R,j1),S,j2))* A[n-i1](r+j1,s+j2),j1=ldegree(ope1N,R)..degree(ope1N,R)),j2=ldegree(ope1N,S)..degree(ope1N,S)),i1=1..-ldegree(ope1N,N))):
print(``):
print(`and in Maple notation`):
print(``):
lprint(A[n](r,s)=add(add(add(normal(coeff(coeff(coeff(ope1N,N,-i1),R,j1),S,j2))* A[n-i1](r+j1,s+j2),j1=ldegree(ope1N,R)..degree(ope1N,R)),j2=ldegree(ope1N,S)..degree(ope1N,S)),i1=1..-ldegree(ope1N,N))):
fi:


if ope1N=FAIL then
 print(`Theorem: For any integers p and q define the Abel-sum type sequence  by`):
print(``):
 print(A[n](r,s)=Sum(F*(r+k)^(k-1+p)*(s-k)^(n-k+q),k=0..n)):
print(``):
print(`Then we have the following functional recurrence relating `, seq(A[n+i1](r,s),i1=0..degree(ope1,N))):

print(add(add(add(normal(coeff(coeff(coeff(ope1,N,i1),R,j1),S,j2))* A[n+i1](r+j1,s+j2),j1=ldegree(ope1 ,R)..degree(ope1 ,R)),j2=ldegree(ope1 ,S)..degree(ope1 ,S)),i1=0..degree(ope1 ,N))=0):
print(``):
print(`and in Maple notation`):
print(``):
lprint(add(add(add(normal(coeff(coeff(coeff(ope1,N,i1),R,j1),S,j2))* A[n+i1](r+j1,s+j2),j1=ldegree(ope1 ,R)..degree(ope1 ,R)),j2=ldegree(ope1 ,S)..degree(ope1 ,S)),i1=0..degree(ope1 ,N))=0):
fi:

print(`Proof: We claim that: `):
print(``):
print(`Let a(n,k,r,s) be the summand on the sum defining A[n](r,s)`):
print(``):
print(`in other words `):
print(``):
print(a(n,k,r,s)=F*(r+k)^(k-1+p)*(s-k)^(n-k+q)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(a(n,k,r,s)=F*(r+k)^(k-1+p)*(s-k)^(n-k+q)):
print(``):
print(`Then the following identity is true`):
print(``):

print(add(add(add(add(normal(coeff(coeff(coeff(coeff(ope,N,i1),R,j1),S,j2),K,i2))* a(n+i1,k+i2,r+j1,s+j2),j1=ldegree(ope ,R)..degree(ope ,R)),j2=ldegree(ope ,S)..degree(ope ,S)),i2=0..degree(ope,K)),i1=0..degree(ope ,N))=0):
print(``):
print(`and in Maple notation:`):
print(``):
lprint(add(add(add(add(normal(coeff(coeff(coeff(coeff(ope,N,i1),R,j1),S,j2),K,i2))* a(n+i1,k+i2,r+j1,s+j2),j1=ldegree(ope ,R)..degree(ope ,R)),j2=ldegree(ope ,S)..degree(ope ,S)),i2=0..degree(ope,K)),i1=0..degree(ope ,N))=0):
print(``):
print(`The proof of this identity is routine (divide by a(n,k,r,s), simplify each term,and now each term is  a rational function. Now  add them all up and verify that they add up to zero.)`):
print(``):

if ope1N=FAIL then
print(`Now sum it from k=0 to k=n, which is the same as from k=-infinity to k=infinity (since it vanishes for k<0 and k>n).  Since the above recurrence is free from k, we get the claimed recurrence. QED.`):
fi:


if ope1N<>FAIL then

print(`Now sum it from k=0 to k=n, which is the same as from k=-infinity to k=infinity (since it vanishes for k<0 and k>n`):
print(``):
print(add(add(add(normal(coeff(coeff(coeff(ope1,N,i1),R,j1),S,j2))* A[n+i1](r+j1,s+j2),j1=ldegree(ope1 ,R)..degree(ope1 ,R)),j2=ldegree(ope1 ,S)..degree(ope1 ,S)),i1=0..degree(ope1 ,N))=0):
print(``):
print(`replacing n by`, n+ldegree(ope1N,N), ` changing variables, and moving a[n](r,s) to the left side,  yields the statement of the thereom. QED.`):
fi:
print(``):
print(`-------------------------------------------------`):
print(``):
print(`This took`, time()-t0, `seconds. `):
print(``):


end:



#FindAbelRec(M,n,r,s,N,R,S,x,MaxOrd): finds all hypergemetric terms of the from x^k/(k!^i*(n-k)^j with 1<=i,j<=M that yields functional recurrences of order <=MaxOrd. 
#it also returns the operators annihilating the Abel-sums induced by F. Try:
#FindAbelRec(2,n,r,s,N,R,S,x,4);
FindAbelRec:=proc(M,n,r,s,N,R,S,x,MaxOrd) local i,j,F,K,ope,gu,k:

gu:={}:

for i from 1 to M do
 for j from 1 to M do
 F:=1/k!^i/(n-k)!^j*x^k:
  ope:=FindOpe(F,n,k,r,s,R,S,N,K,MaxOrd):
 if ope<>FAIL then
  gu:=gu union {[F,subs(K=1,ope)]}:
 fi:
od:
od:
gu:
end:




####start differential version

diff1:=proc(f,x,i) 
if i=0 then
 f:
else
 diff(f,x$i):
fi:
end:


#ApplyOpeD(ope,N,Dr,Ds,n,r,s,F): applies the operator ope(N,Dr,Ds) to F(n,r,s) divided by F(n,r,s). Try:
#ApplyOpeD(N+Dr+Ds,N,Dr,Ds,n,r,s,n*r*s);
ApplyOpeD:=proc(ope,N,Dr,Ds,n,r,s,F) local gu,i1,i2,i3,mu,mu1,f2:
gu:=0:

for i1 from 0 to degree(ope,N) do

for i2 from 0 to degree(ope,Dr) do

for i3 from 0 to degree(ope,Ds) do
mu:=coeff(coeff(coeff(ope,N,i1),Dr,i2),Ds,i3):
f2:=subs(n=n+i1,F):
f2:=diff1(f2,r,i2):
f2:=diff1(f2,s,i3):

mu1:=normal(f2/F):

 gu:=gu+mu*mu1:

 gu:=normal(gu):
od:
od:
od:
simplify(gu):

end:


#ApplyOpeDnaked(ope,N,Dr,Ds,n,r,s,F): applies the operator ope(N,Dr,Ds) to F(n,r,s) divided by F(n,r,s). Try:
#ApplyOpeDnaked(N+Dr+Ds,N,Dr,Ds,n,r,s,n*r*s);
ApplyOpeDnaked:=proc(ope,N,Dr,Ds,n,r,s,F) local gu,i1,i2,i3,mu,mu1,f2:
gu:=0:

for i1 from 0 to degree(ope,N) do

for i2 from 0 to degree(ope,Dr) do

for i3 from 0 to degree(ope,Ds) do
mu:=coeff(coeff(coeff(ope,N,i1),Dr,i2),Ds,i3):
f2:=subs(n=n+i1,F):
f2:=diff1(f2,r,i2):
f2:=diff1(f2,s,i3):

mu1:=normal(f2):

 gu:=gu+mu*mu1:

 gu:=normal(gu):
od:
od:
od:
simplify(gu):

end:



#ApplyOpeSeqD(ope,N,Dr,Ds,n,r,s,L): applies the operator ope(N,R,S,n,r,s) to the sequence L.  

ApplyOpeSeqD:=proc(ope,N,Dr,Ds,n,r,s,L) local lu,lu1,n1,i1:
lu:=[]:
for n1 from 1 to nops(L)-degree(ope,N) do

 lu1:=normal(add(ApplyOpeDnaked(subs(n=n1,coeff(ope,N,i1)),N,Dr,Ds,n,r,s,L[n1+i1]),i1=0..degree(ope,N))):
 lu:=[op(lu),lu1]:
 od:
lu:
end:



#PrintOpe(ope,n,r,s,N,R,S,A): inputs a difference operator in n,r,s, with shift operators N,R,S, applied to A[n](r,s): outputs it in plain Engllsh. Try:
#PrintOpe(r*N*R*S+s*R+r*S,n,r,s,N,R,S,A);
PrintOpe:=proc(ope,n,r,s,N,R,S,A) local i1,j1,j2:
add(add(add(normal(coeff(coeff(coeff(ope,N,i1),R,j1),S,j2))* A[n+i1](r+j1,s+j2),j1=ldegree(ope ,R)..degree(ope ,R)),j2=ldegree(ope ,S)..degree(ope ,S)),i1=ldegree(ope,N)..degree(ope ,N)):
end:



#PrintOpeD(ope,n,r,s,N,Dr,Ds,A): inputs a recurrence (in n) -and (differential) operator in,r and s, with shift operator N and differnetiation Dr, Ds, respectively applied to A[n](r,s): outputs it in plain Engllsh. Try:
#PrintOpeD(r*N*Dr+s*Ds,n,r,s,N,Dr,Ds,A);
PrintOpeD:=proc(ope,n,r,s,N,Dr,Ds,A) local i1,j1,j2,gu:
gu:=0:
for i1 from ldegree(ope,N) to degree(ope,N) do
 gu:=gu+coeff(coeff(coeff(ope,N,i1),Dr,0),Ds,0)*A[n+i1](r,s):
 for j1 from 1 to degree(ope,Dr) do
   gu:=gu+coeff(coeff(coeff(ope,N,i1),Dr,j1),Ds,0)*diff(A[n+i1](r,s),r$j1):
 od:
 for j2 from 1 to degree(ope,Ds) do
   gu:=gu+coeff(coeff(coeff(ope,N,i1),Dr,0),Ds,j2)*diff(A[n+i1](r,s),s$j2):
 od:
for j1 from 1 to degree(ope,Dr) do
 for j2 from 1 to degree(ope,Dr) do
 gu:=gu+coeff(coeff(coeff(ope,N,i1),Dr,j1),Ds,j2)*diff(diff(A[n+i1](r,s),s$j2),r$j1):
od:
od:
od:
gu:
end:




#PaperD(F,n,k,MaxOrd,KER): inputs a hypergeometric term F in discrete variables n and k and symbols. outputs a paper with just the statement of the differential recurrence equations,
#first order with respect to n and up to order MaxOrder one with respect to r and s
#satisfied by
#A[n](r,s):=Sum(F*KDER, k=0..n); for any integers p and q. Try: 
#PaperD(binomial(n,k),n,k,r,s,2, (r+k)^(k-1+p)*(s-k)^(n-k+q));`):


PaperD:=proc(F,n,k,r,s,MaxOrd,KER) local N,Dr,Ds,A,ope,t0,x:
t0:=time():

ope:=DiffPair(F,n,k,r,s,N,Dr,Ds,MaxOrd,KER):


if ope=FAIL then
 print(`Try to make MaxOrd bigger`):
 RETURN(FAIL):
fi:


 print(`Theorem:  define the Abel-sum type sequence  by`):

print(``):
 print(A[n](r,s)=Sum(F*KER,k=0..n)):
print(``):
print(`and in Maple notation`):
print(``):
 lprint(A[n](r,s)=Sum(F*KER,k=0..n)):
print(``):

print(`Then we have the following two differential-recurrence  equations, relating A[n](r,s) and A[n+1](r,s) , the first one with respect to r, the second with respect to s`):

print(``):
print(PrintOpeD(ope[1],n,r,s,N,Dr,Ds,A)=0):
print(``):
print(PrintOpeD(ope[2],n,r,s,N,Dr,Ds,A)=0):
print(``):


print(`and in Maple notation`):

print(``):
lprint(PrintOpeD(ope[1],n,r,s,N,Dr,Ds,A)=0):
print(``):
lprint(PrintOpeD(ope[2],n,r,s,N,Dr,Ds,A)=0):
print(``):



print(`-------------------------------------------------`):
print(``):
print(`This took`, time()-t0, `seconds. `):
print(``):




end:



#FindOpeD1g(F,n,k,r,N,Dr,ORDn,ORDr,KER) :
# inputs a hypergeometric term F in n and k, also a kernel KER, and outputs a k-free operator  of (n,r,N,Dr)
#of order ORDn in n (i.e. degree in N), order ORDr in r (i.e. degree in Dr))
#annihilating
#F*KER
#FindOpeD1g(binomial(n,k),n,k,r,s,N,Dr,1,1, (r+k)^(k-1+p)*(s-k)^(n-k+q));
FindOpeD1g:=proc(F,n,k,r,N,Dr,ORDn,ORDr,KER) local ope,i1,i2,a,var,gu,F1,i,eq,tov,var1:

ope:=add(add(a[i1,i2]*N^i1*Dr^i2, i1=0..ORDn),i2=0..ORDr):

var:={seq(seq(a[i1,i2], i1=0..ORDn),i2=0..ORDr)}:

F1:=convert(F,factorial)*KER:
gu:=0:

for i1 from 0 to ORDn do
 for i2 from 0 to ORDr do
    gu:=gu+ a[i1,i2]*normal(simplify(diff1(subs(n=n+i1,F1),r,i2)/F1)):
od:
od:


gu:=numer(gu):


eq:={seq(coeff(gu,k,i),i=0..degree(gu,k))}:

var:=solve(eq,var):

if subs(var,ope)=0 then
 RETURN(FAIL):
fi:


tov:={}:

for var1 in var do
 if op(1,var1)=op(2,var1) then
   tov:=tov union {op(1,var1)}:
 fi:
od:

if nops(tov)>1 then
 print(`Too much slack`):
fi:

ope:=subs(var,ope):

ope:=subs(tov[1]=1,ope):
ope:
end:






#FindOpeDg(F,n,k,r,N,Dr,MaxOrd,KER) :
# inputs a hypergeometric term F in n and k, also a kernel KER, and outputs a k-free operator  of (n,r,N,Dr)
#sum of theorder ORDn in n (i.e. degree in N) and the order ORDr in r (i.e. degree in Dr))<=MaxOrd, or returns FAIL
#annihilating
#F*KER
#FindOpeDg(binomial(n,k),n,k,r,N,Dr,2, (r+k)^(k-1+p)*(s-k)^(n-k+q));
FindOpeDg:=proc(F,n,k,r,N,Dr,MaxOrd,KER) local  MaxOrd1,ORDn, ORDr,gu:

for MaxOrd1 from 2 to MaxOrd do
 for ORDn from 1 to MaxOrd1 do
  ORDr:=MaxOrd1-ORDn:
 gu:=FindOpeD1g(F,n,k,r,N,Dr,ORDn,ORDr, KER):
  if gu<>FAIL then
    RETURN(gu):
  fi:
 od:
od:

FAIL:
end:




#DiffPair(F,n,k,r,s,N,Dr,Ds,MaxOrd,KER) :
# inputs a hypergeometric term F in n and k, also a kernel KER, and outputs a k-free operators  of (n,r,N,Dr) and (n,s,N,Ds) respectively
#where sum of recurrence order and the differntial order is <=MaxOrd
#annihilating
#F*KER
#If none exists, then it returns FAIL. Try:
#DiffPair(binomial(n,k),n,k,r,s,N,Dr,Ds,,2, (r+k)^(k-1+p)*(s-k)^(n-k+q));
DiffPair:=proc(F,n,k,r,s,N,Dr,Ds,MaxOrd,KER)  local ope1,ope2:
option remember:

ope1:=FindOpeDg(F,n,k,r,N,Dr,MaxOrd,KER):

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

ope2:=FindOpeDg(F,n,k,s,N,Ds,MaxOrd,KER):


if ope2=FAIL then
 RETURN(FAIL):
fi:

[ope1,ope2]:

end:



#PolAg(F,n,k,n1,KER): Given a bi-variate hypergemetric term F in n,k, and a kernel KER=KER(n,k,r,s) depending on n,k,r,s, outputs
#Sum(F(n1,k)*KER(n1,k,r,s),k=0..n1); Try:
#PolAg(binomial(n,k),n,k,(r+k)^(k-1)*(s-k)^(n-k),5):
PolAg:=proc(F,n,k,KER,n1) local k1:
factor(normal(add(eval(subs({n=n1,k=k1},F))*subs({n=n1,k=k1},KER),k1=0..n1))):
end:

#SeqAg(F,n,k,KER,N1): Given a bi-variate hypergemetric term F in n,k, and a kernel KER=KER(n,k) outputs 
#[seq(Sum(F(n1,k)*KER(n1,k),k=0..n1),n1=0..N1); Try
#SeqAg(binomial(n,k),n,k,(r+k)^(k-1)*(s-k)^(n-k),10);
SeqAg:=proc(F,n,k,KER,N1) local n1:
[seq(PolAg(F,n,k,KER,n1),n1=1..N1)]:
end:



#CheckDiffPair(F,n,k,r,s,N,Dr,Ds,MaxOrd,KER,KAMA): checks DiffPair(F,n,k,r,s,N,Dr,Ds,MaxOrd,KER) up to KAMA terms Try
#CheckDiffPair(binomial(n,k),n,k,r,s,N,Dr,Ds,2,(r+k)^(k-1)*(s-k)^(n-k),20);
CheckDiffPair:=proc(F,n,k,r,s,N,Dr,Ds,MaxOrd,KER,KAMA)  local ope,lu:
ope:=DiffPair(F,n,k,r,s,N,Dr,Ds,MaxOrd,KER):
if ope=FAIL then
 RETURN(FAIL):
fi:
lu:=SeqAg(F,n,k,KER,KAMA);
[ApplyOpeSeqD(ope[1],N,Dr,Ds,n,r,s,lu) , ApplyOpeSeqD(ope[2],N,Dr,Ds,n,r,s,lu) ]:
end: