read `W:\\mPOL.txt`:

HelpC:=proc(): print(`Celine1(F,n,k,N,K,degn,degN,degK)`):
print(`ApplyOper(F,n,k,P,N,K,n1,k1)`):
print(`CheckCeline(F,n,k,P,N,K,A)`): 
print(`Celine(F,n,k,N,K,LimitOfSum),ApplyOperS(F,n,k,P,N,K) `):
end:


#Celine1(F,n,k,N,K,degn,degN,degK): Given an expression 
#in n and k, variables (symbols) n and k, symbols
#N and K denoting the fundamental shift operators
#in n and k resp. (NF(n,k):=F(n+1,k), KF(n,k):=F(n,k+1)
#outputs a conjectured linear recurrence operator
#of the form P(N,K,n) )of degrees degn,degN,degK
#resp.) in each of
#n,N,K
#(free of k) that annihilates
#F(n,k)

Celine1:=proc(F,n,k,N,K,degn,degN,degK) local P,a,var,k1,n1,eq,c,
var1,eq1,aek,i:


P:=GenPol([n,N,K],[degn,degN,degK],a,0):

var:=P[2]:

P:=P[1]:


eq:={}:

c:=0:

for n1 from 1 do
 for k1 from 0 to n1 do

 eq1:=ApplyOper(F,n,k,P,N,K,n1,k1):
 
  if eq1<>0 then
     c:=c+1:
     eq:=eq union {eq1}:
   fi:


if nops(eq)-nops(var)>8 then
 var1:=solve(eq,var):

   aek:={}:

     for i from 1 to nops(var1) do
        if op(1,var1[i])=op(2,var1[i]) then
          aek:=aek union {op(1,var1[i])}:
        fi:
   od:


#if nops(aek)>1 then
# print(`Too much slack, reduce one of the degrees`):
#fi:



 P:=subs(var1,P):
  if P=0 then
    RETURN(FAIL):
   else

    P:=subs({aek[1]=1,seq(aek[i]=0,i=2..nops(aek))},P):

      if not ApplyOperS(F,n,k,P,N,K)<>0 then
       RETURN(FAIL):
       else
      RETURN(P):
      fi:
   fi:
fi:
 
 od:
od:




end:


#ApplyOper(F,n,k,P,N,K,n1,k1): P(n,N,K)F(n1,k1)
ApplyOper:=proc(F,n,k,P,N,K,n1,k1)
local s,i,j:

s:=0:

for i from ldegree(P,N) to degree(P,N) do
 for j from ldegree(P,K) to degree(P,K) do
   s:=s+
  subs(n=n1,coeff(coeff(P,N,i),K,j))*eval(subs({n=n1+i,k=k1+j},F)):
 od:
od:

s:
   




end:


#CheckCeline(F,n,k,P,N,K,A): Checks that the operator
#P indeed annihilates the discrete (expression) function F(n,k)
#for 0<=k<=n<=A
CheckCeline:=proc(F,n,k,P,N,K,A) local n1,k1:

 evalb({ seq(seq(ApplyOper(F,n,k,P,N,K,n1,k1),k1=0..n1),
n1=1..A) }={0}):


end:


#Celine(F,n,k,N,K,Arvind): inputs an expression F, in n and k,
# and outputs a recurrence operator P(n,N,K) annihilating
#F with degn+degN+degK<=Arvind
Celine:=proc(F,n,k,N,K,Arvind) local A,ope,degn,degK,degN:


for A from 0 to Arvind do
 for degn from 0 to A do
  for degN from 0 to A-degn do
     degK:=A-degn-degN: 

     ope:=Celine1(F,n,k,N,K,degn,degN,degK): 
       if ope<>FAIL then
          RETURN(ope):
       fi:
   od:
 od:
od: 

FAIL:

end:


#ApplyOperS(F,n,k,P,N,K): given a hypergeometric
#discrete function F(n,k) (F(n+i,k+j)/F(n,k) must
#be a rational function of n and k)
#applies the operator P(n,N,K) to F (divided by F)
#P=Sum p_{i,j}(n) N^i *K^j 
ApplyOperS:=proc(F,n,k,P,N,K) local s,i,j:
s:=0:

   for i from 0 to degree(P,N) do
      for j from 0 to degree(P,K) do
          s:=normal(s+coeff(coeff(P,N,i),K,j)*
             expand(subs({n=n+i, k=k+j},F)/F)):
      od:
  od:

s:

end:


#CelineWZ(F,n,k,N,Arvind): inputs a hypergeometric F(n,k)
#and shift operator N (in n) outputs a
#linear recurrence operator with poly coeff.
#P(n,N) that annihilates a(n):=sum(F(n,k),k=-infinity..infinity)
#as well as the certificate ( a rational function of (n,k))
#such that G(n,k):=R(n,k)*F(n,k) satisfies
#the general WZ equation
#ope(n,N)F(n,k)=G(n,k+1)-G(n,k)
CelineWZ:=proc(F,n,k,N,K,Arvind) 

end: