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


print(`First Written: July 2,2004: tested for Maple 8  `):
print(`Version of July 2, 2004:  `):
print():
print(`This is MultiZeilberger, A Maple package`):
print(`accompanying the article `):
print(`The Multi-Zeilberger Algorithm`):
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg .`):
print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`):
print():
print(` type "ezra1();". for a list of all functions.`):
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name)" `):
print():
 
 

ezra1:=proc()
if args=NULL then
print(`The procedures are : hafokh, MulZeil , MulZeil1, MulZeilCb`):
print():
else
 ezra(args):
fi:
end: 
 
ezra:=proc()
if args=NULL then
print(`MultiZeilberger`):
print(`A Maple package for )`):
print(`finding recurrences for hypergeometric multi-sums`):
print(`Using the Mohammed-Zeilberger Extension to Multi-sums of the Zeilberger`):
print(`algorithm`):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures: MulZeil , MulZeil1, MulZeilCb`):
print():
fi:
 
if nargs=1 and args[1]=hafokh then
print(`hafokh(F,kList,n): given a summand F, and a list of summation variables kList`):
print(`and another (discrete) variable n (the argument of the right side)`):
print(`decides whether it is proper-hypegeometric (if it is not, it returns FAIL)`):
print(`if it is, outputs it in canonical notation`):
fi:

if nargs=1 and args[1]=MulZeil1 then
print(`MulZeil1(F,kList,n,N,para): Given a hypegeometric term F`):
print(`in the variable n and the variables of the list kList,`):
print(`and a set of auxiliary parameters para, of other variables`):
print(`(para may always be set to be the empty set {} )`):
print(`finds a recurrence operator ope(n,N) and a list of discrete functions`):
print(`that are the pre-certificates (which are polynomials in kList)`):
print(`For example, try:`):
print(`MulZeil1(n!/i!/j!/(n-i-j)!,[i,j],n,N,{}); and `):
print(`MulZeil1(n!/i!/(m+i)!/j!/(n-i-j)!,[i,j],n,N,{m}); `):
fi:


if nargs=1 and args[1]=MulZeil then
print(`MulZeil(F,kList,n,N,para): Given a hypegeometric term F`):
print(`in the variable n and the variables of the list kList`):
print(`and a set of auxiliary parameters para, of other variables`):
print(`(para may always be set to be the empty set {} )`):
print(`finds a recurrence operator ope(n,N) and a list of discrete functions`):
print(`that are the certificate`):
print(`In other words: ope(n,N)F=sum( (K_i-1)(Cert[i]*F),i=1..nops(kList)`):
print(`where N is the shift operator in n and K_i the shift operator`):
print(`in kList[i]`):
print(`For example, try:`):
print(`MulZeil(n!/i!/j!/(n-i-j)!,[i,j],n,N,{}); and `):
print(`MulZeil(n!/i!/(m+i)!/j!/(n-i-j)!,[i,j],n,N,{m}); `):
fi:



if nargs=1 and args[1]=MulZeilCb then
print(`MulZeilCb(F,kList,kLimits,n,N,para)`): 
print(`First does what MulZeilC does, but then tries to find`):
print(` a smaller-order recurrence, empirically`):
print(`by summing F over the kLimits (i is tacicty assumed that`):
print(`F vanishes outside the summation region`):
print(`If there is no lower-order recurrence, the last output is 1`):
print(`otherwise, the first output is the empirically-found recurrence`):
print(`the last output is the operator that if you left-multiply it`):
print(`by the first output you would get the operator given by MulZeil`):
print(`and the second output is like in MulZeil`):
print(`For example, try:`):
print(`MulZeilCb(n!/i!/j!/(n-i-j)!,[i,j],[[0,n],[0,n-i]],n,N,{}); and `):
print(`Let F:=(i+j)!*(n-i)!*(n-j)!/i!^2/j!^2/(n-i-j)!^2;, then try`):
print(`MulZeilCb(F,[i,j],[[0,n],[0,n-i]],n,N,{}); `):
fi:


end:




rf:=proc(a,n) local i: if n<0 then  0 elif n=0 then  1
else convert([seq(a+i,i=0..n-1)],`*`) fi: end:

RatioH:=proc(Num,Den,n,L,i) local j:
convert([seq( rf(Num[j]+1,i*coeff(Num[j],n,1)),j=1..nops(Num))],`*`)*
convert([seq( rf(
subs(n=n+i,Den[j])+1,(L-i)*coeff(Den[j],n,1)),j=1..nops(Den))],`*`):
end:

c:=proc(Num,Den,POL,e,n,L) local i:convert(
[seq(e[i]*subs(n=n+i,POL)*RatioH(Num,Den,n,L,i),i=0..L)],`+`):
end:

a:=proc(Num,Den,zList,kList,n,L) local j,gu,i:

for i from 1 to nops(kList) do
 gu[i]:=zList[i]:
od:

for j from 1 to nops(Num) do
  for i from  1 to nops(kList) do
    if coeff(Num[j],kList[i],1) > 0 then
     gu[i]:=gu[i]*rf(Num[j]+1,coeff(Num[j],kList[i],1)):
    fi:
  od:
od:

for j from 1 to nops(Den) do
 for i from 1 to nops(kList) do 
  if coeff(Den[j],kList[i],1) < 0 then
   gu[i]:=gu[i]*rf(subs({n=n+L,kList[i]=kList[i]+1},
   Den[j])+1,-coeff(Den[j],kList[i],1)):
  fi:
 od:
od:
[seq(gu[i],i=1..nops(kList))]:
end:



b:=proc(Num,Den,kList,n,L) local i,j,gu:

for i from 1 to nops(kList) do
 gu[i]:=1:
od:

for j from 1 to nops(Num) do
 for i from 1 to nops(kList)  do
  if coeff(Num[j],kList[i],1) < 0 then
   gu[i]:=gu[i]*
       rf(subs(kList[i]=kList[i]+1,Num[j])+1,-coeff(Num[j],kList[i],1)):
  fi:
 od:
od:

for j from 1 to nops(Den) do
 for i from 1 to nops(kList) do 
  if coeff(Den[j],kList[i],1) > 0 then
   gu[i]:=gu[i]*rf(subs({n=n+L},Den[j])+1,coeff(Den[j],kList[i],1)):
  fi:
 od:
od:
[seq(gu[i],i=1..nops(kList))]:
end:

Equ1:=proc(F1,kList,n,L,e,X1,X) 
local Num,Den,POL,zList,aList, bList,i:
Num:=F1[1]: Den:=F1[2]: POL:=F1[3]: zList:=F1[4]:
aList:=a(Num,Den,zList,kList,n,L):
bList:=b(Num,Den,kList,n,L):

add(
aList[i]*X1[i]-
subs(kList[i]=kList[i]-1,bList[i])*X[i],i=1..nops(kList))-
c(Num,Den,POL,e,n,L):
end:

yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:




#IV1(d,n): all the vectors of non-negative integres of length d whose
#sum is n
IV1:=proc(d,n)
local gu,i,gu1,i1:
if d=0 then
 if n=0 then
   RETURN({[]}):
 else
    RETURN({}):
fi:
fi:
 
gu:={}:

for i from 0 to n do
 gu1:=IV1(d-1,n-i):
 gu:=gu union {seq([op(gu1[i1]),i],i1=1..nops(gu1))}:
od:
gu:
end:

#IV(d,n): all the vectors of non-negative integres of length d whose
#sum is <=n
IV:=proc(d,n) local i: {seq(op(IV1(d,i)),i=0..n)}: end:

#GenPol(kList,a,deg): The generic polynomial in kList of
#degree deg,  using the indexed variable a, followed by the set
#of coeffs.
GenPol:=proc(kList,a,deg)
local gu,i,i1:
gu:=IV(nops(kList),deg):
convert([seq(a[i]*convert([seq(kList[i1]^gu[i][i1],i1=1..nops(kList))],`*`),
i=1..nops(gu))],`+`),{seq(a[i],i=1..nops(gu))}:

end:

#Extract1(POL,kList): extracts all the coeffs. of a POL in kList
Extract1:=proc(POL,kList) local k1,kList1,i:
k1:=kList[nops(kList)]:
kList1:=[op(1..nops(kList)-1,kList)]:
if nops(kList)=1 then
 RETURN({seq(coeff(POL,k1,i),i=0..degree(POL,k1))}):

else
 RETURN({seq(
op(Extract1(coeff(POL,k1,i),kList1)), i=0..degree(POL,k1))}):
fi:

end:


Z1Ld:=proc(F,kList,n,L,N,degX,para) local e,x,j,X,gu,var1,
var,ope,eq,zList,ku,i,ope1,kap,
i1,Sols,j1,Sols1,X1,Sols2,kvar,Y,Y1,eqTry,varTry:
global t0:

ope:=convert([seq(e[j]*N^j,j=0..L)],`+`):

zList:=F[4]:
var1:=[seq(e[j],j=0..L)]:

X:=[]:
for i from 1 to nops(kList) do
 ku:=GenPol(kList,x[i],degX[i]):
 X:=[op(X),ku[1]]:
 var1:=[op(var1), op(ku[2])]:
od:



gu:=Equ1(F,kList,n,L,e,Y1,Y) :

for i from 1 to nops(kList) do
 gu:=subs({Y1[i]=subs(kList[i]=kList[i]+1,X[i]),Y[i]=X[i]},gu):
od:


eq:=Extract1(gu,kList):

eqTry:=eq:
eqTry:=subs(n=5/11,eq):
for i from 1 to nops(zList) do
 if type(zList[i],symbol) then
  eqTry:=subs(zList[i]=5/(i^2+3),eqTry):
 fi:
od:

for i from 1 to nops(para) do
   eqTry:=subs(para[i]=7/(i^3+5),eqTry):
od:

varTry:=LS(eqTry,var1):
if varTry={} or
   {seq([op(1..L+1,varTry[i])],i=1..nops(varTry))}={[0$(L+1)]} then
RETURN(FAIL):
fi:

print(`The time is:`, time()-t0):
print(`There is great hope for a recurrence of order`,L,` please be patient`):
var:=LS(eq,var1):
Sols:={
seq(
[subs({seq(var1[j]=var[i1][j],j=1..nops(var1))},ope),
[seq(subs({seq(var1[j]=var[i1][j],j=1..nops(var1))},X[j1]),j1=1..nops(X))]]
,i1=1..nops(var))}:

Sols1:={}:
for i from 1 to nops(Sols) do
 ope1:=normal(Sols[i][1]):
   X1:=normal(Sols[i][2]):
 if ope1<>0 then
  kap:=yafe(ope1,N):
  ope1:=kap[2]:
  Sols1:=Sols1 union
  {[ope1,[seq(normal(X1[i1]/kap[1]),i1=1..nops(X1))]]}:
 fi:
od:


if Sols1={} then
 RETURN(FAIL)
fi:

if nops(Sols1)>1 then
 Sols2:={Sols1[1]}:
 kvar:={Sols1[1][1]}:
 for i from 2 to nops(Sols1) do
  if not member(Sols1[i][1],kvar) then
    kvar:=kvar union {Sols1[i][1]}:
    Sols2:=Sols2 union {Sols1[i]}:
  fi:
 od:
 if nops(Sols2)>1 then
  print(`There are more than one operator`):
 else
   print(`There is just one operator but several cerficates`):
 fi:
print(`The whole thing took:`, time()-t0):
 RETURN(Sols2):
else
print(`The whole thing took:`, time()-t0):
 RETURN(Sols1):
fi:

end:

Z1L:=proc(F,kList,n,L,N,para) 
local degX,gu,i,j,aList,bList,c1,Num,Den,POL,zList,e:
Num:=F[1]: Den:=F[2]: POL:=F[3]: zList:=F[4]:
aList:=a(Num,Den,zList,kList,n,L):
bList:=b(Num,Den,kList,n,L):
c1:=c(Num,Den,POL,e,n,L):

degX:=[]:
for i from 1 to nops(kList) do
degX:=[op(degX),max(0,degree(c1,{op(kList)})-max(degree(aList[i],{op(kList)}),
            degree(bList[i],{op(kList)})))]:
od:

gu:=Z1Ld(F,kList,n,L,N,degX,para):

 if gu<>FAIL then
  RETURN(gu):
 fi:

for j from 1 to 1 do
 degX:=[seq(degX[i]+1,i=1..nops(degX))]:
 gu:=Z1Ld(F,kList,n,L,N,degX,para):

 if gu<>FAIL then
  RETURN(gu):
 fi:
od:
FAIL:
end:


Z1:=proc(F,kList,n,N,para) local L,gu:
gu:=Z1L(F,kList,n,0,N,para):
for L from 1 while gu=FAIL do
 gu:=Z1L(F,kList,n,L,N,para):
od:
gu:
end:



####converting to the input needed by Z1###
FindFactorial:=proc(Prod1)
local Prod2,i,lu,lu1,lu2:
Prod2:=factor(Prod1):
if type(Prod2,function) and op(0,Prod2)=factorial then
 RETURN(Prod1):
fi:

for i from 1 to nops(Prod2) do
 lu:=op(i,Prod2):
 if type(lu,function) and op(0,lu)=factorial  then
  RETURN(lu):
 fi:

 if type(lu,`^`)  then
   lu1:=op(1,lu):
   lu2:=op(2,lu):

     if type(lu1,function) and op(0,lu1)=factorial then
          if not type(lu2,integer) and lu2>0 then
              RETURN(FAIL):
          else
              RETURN(lu1):
          fi:
     fi: 
 fi:
od:
FAIL:
end:


FindPower:=proc(Prod1,kList)
local Prod2,i,lu:
Prod2:=factor(Prod1):
if Prod2=1 then
 RETURN(FAIL):
fi:
if type(Prod2,`^`) and degree(op(2,Prod1),{op(kList)})=1 then
 RETURN(Prod1):
fi:

for i from 1 to nops(Prod2) do
 lu:=op(i,Prod2):
 if type(lu,`^`) and degree(op(2,lu),{op(kList)})=1 then
  RETURN(lu):
 fi:
od:
FAIL:
end:




hafokh:=proc(F,kList,n)
local F1,Num,Den,POL,z,F1t,F1b,lu,lu1,lu2,i:

F1:=factor(normal(F)):

F1t:=factor(numer(F1)):
F1b:=factor(denom(F1)):

Num:=[]:
lu:=FindFactorial(F1t):
while lu<>FAIL do
lu1:=op(lu):
 if degree(lu1,{n,op(kList)})>1 or not type(coeff(lu1,n,1),integer) then
    RETURN(FAIL):
 fi:  

 for i from 1 to nops(kList) do
  if not type(coeff(lu1,kList[i],1),integer) then
   RETURN(FAIL):
  fi:
 od:
Num:=[op(Num),lu1]:
F1t:=normal(F1t/lu):
lu:=FindFactorial(F1t):
od:

for i from 1 to nops(kList) do
 z[i]:=1:
od:

lu:=FindPower(F1t,kList):
while lu<>FAIL do
lu1:=op(1,lu):
lu2:=op(2,lu):

if degree(lu2,{op(kList)})>1 then
  RETURN(FAIL):
else
 for i from 1 to nops(kList) do
  z[i]:=z[i]*lu1^coeff(lu2,kList[i],1):
  F1t:=normal(simplify(normal(F1t/lu1^(kList[i]*coeff(lu2,kList[i]))))):
 od:
fi:

lu:=FindPower(F1t,kList):
od:

 
POL:=F1t:

Den:=[]:

lu:=FindFactorial(F1b):
while lu<>FAIL do
lu1:=op(lu):
 if degree(lu1,{n,op(kList)})>1 or 
   not type(coeff(lu1,n,1),integer) then
    RETURN(FAIL):
 fi:  

 for i from 1 to nops(kList) do
  if not type(coeff(lu1,kList[i],1),integer) then
   RETURN(FAIL):
  fi:

 od:
Den:=[op(Den),lu1]:


F1b:=normal(F1b/lu):
lu:=FindFactorial(F1b):
od:

lu:=FindPower(F1b,kList):
while lu<>FAIL do
lu1:=op(1,lu):
lu2:=op(2,lu):

if degree(lu2,{op(kList)} ) >1 then
  RETURN(FAIL):
else

 for i from 1 to nops(kList) do
  z[i]:=z[i]/lu1^coeff(lu2,kList[i],1):
  F1b:=normal(simplify(normal(F1b/lu1^(kList[i]*coeff(lu2,kList[i],1))))):
 od:

fi:

lu:=FindPower(F1b,kList):
od:

POL:=normal(POL/F1b):

for i from 1 to nops(Num) do 
if coeff(Num[i],n,1)<0 then
 RETURN(FAIL):
fi:
od:

for i from 1 to nops(Den) do 
if coeff(Den[i],n,1)<0 then
 RETURN(FAIL):
fi:
od:


if degree(POL,{n,op(kList)})=FAIL then
 RETURN(FAIL):
fi:

[Num,Den,POL,[seq(z[i],i=1..nops(kList))]]:

end:

##End conversion part


MulZeil1:=proc(F,kList,n,N,para) local F1,lu:
global t0:
t0:=time():
F1:=normal(convert(F,factorial)):
F1:=hafokh(F1,kList,n):
if F1=FAIL then
 print(`Bad input`):
 RETURN(FAIL):
fi:
lu:=Z1(F1,kList,n,N,para):

if nops(lu)>1 then
 print(`There are more than one operator`):
 RETURN(lu):
else
 RETURN(lu[1]):
fi:
end:




#Linear solver
LS:=proc(eq,var)
local i,var1,ind,X,j:
var1:=SolveTools[Linear](eq,{op(var)}):

ind:={}:

for i from 1 to nops(var1) do
 if op(1,var1[i])=op(2,var1[i]) then
   ind:=ind union {op(1,var1[i])}:
 fi:
od:
X:=[seq(expand(subs(var1,var[i])),i=1..nops(var))]:
{seq([seq(coeff(X[i],ind[j],1),i=1..nops(X))],j=1..nops(ind))}:
end:
#End Linear solver


#MulZeil(F,kList,n,N,para):The operators plus the certificates
MulZeil:=proc(F,kList,n,N,para) local F1, lu,kadima,X,ope,i,L,bList:
option remember:
lu:=MulZeil1(F,kList,n,N,para):

if type(lu,set) then
 lu:=lu[1]:
fi:

ope:=lu[1]:
L:=degree(ope,N):

F1:=hafokh(convert(F,factorial),kList,n):

bList:=b(F1[1],F1[2],kList,n,L):

kadima:=RatioH(F1[1],F1[2],n,L,0):


X:=lu[2]:
ope,[seq(X[i]/kadima*subs(kList[i]=kList[i]-1,bList[i]),i=1..nops(kList))]:
end:


#CheckZeilC0(F,kList,n,N,n0,kList0,para): Checks ZeilC at n=n0 and kList=kList0
CheckZeilC0:=proc(F,kList,n,N,n0,kList0) local lu,cert,i,i1,ope,gu:
lu:=MulZeil(F,kList,n,N,para):
ope:=lu[1]:
cert:=lu[2]:

gu:=0:

for i from 0 to degree(ope,N) do
 gu:=gu+subs({n=n0+i,seq(kList[i1]=kList0[i1],i1=1..nops(kList))},F)*
 subs({n=n0,seq(kList[i1]=kList0[i1],i1=1..nops(kList))},coeff(ope,N,i)):
od:
gu:=eval(gu):

for i from 1 to nops(cert) do
gu:=gu+
subs({n=n0,seq(kList[i1]=kList0[i1],i1=1..nops(kList))},F*cert[i]):
gu:=eval(gu):

gu:=gu-
subs(
{n=n0,seq(kList[i1]=kList0[i1],i1=1..i-1),
kList[i]=kList0[i]+1, seq(kList[i1]=kList0[i1],i1=i+1..nops(kList))},
F*cert[i]):
gu:=eval(gu):
od:
gu:
end:

CheckZeilCn:=proc(F,kList,n,N,N0,N1,K0,K1,para) local k0,n0,i1:
{seq(seq(CheckZeilC0(F,kList,n,N,n0,[seq(k0+i1^2,i1=1..nops(kList))],
para),
n0=N0..N1),k0=K0..K1)}:
end:


#CheckZeilCs(F1,kList,n,N,n0,kList0,para):
# Checks ZeilC at n=n0 and kList=kList0
CheckZeilCs:=proc(F1,kList,n,N,para) local lu,cert,i,F,ope,gu:
lu:=MulZeilC(F1,kList,n,N,para):
ope:=lu[1]:
cert:=lu[2]:

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

for i from 0 to degree(ope,N) do
 gu:=gu+normal(simplify(subs(n=n+i,F)/F))*coeff(ope,N,i):
od:

for i from 1 to nops(cert) do
gu:=gu+cert[i]-subs(kList[i]=kList[i]+1,cert[i])*
simplify(subs(kList[i]=kList[i]+1,F)/F):
od:
normal(gu):
end:



#SumSet(n,kList,kLimits,n0): the set of summationn given
#by kLimits
SumSet:=proc(n,kList,kLimits,n0)
local i1,j1,gu,gu1,low1,high1,kLimits1,kList1:
if nops(kList)<>nops(kLimits) or kList=[] then 
ERROR(`Bad input`):
fi:
low1:=subs(n=n0,kLimits[1][1]):
high1:=subs(n=n0,kLimits[1][2]):
if nops(kList)=1 then
 RETURN({seq([i1],i1=low1..high1)}):
fi:

gu:={}:

for i1 from low1 to high1 do
 kLimits1:=subs(n=n0,[op(2..nops(kLimits),kLimits)]):
 kList1:=[op(2..nops(kList),kList)]:
 gu1:=SumSet(kList[1],kList1,kLimits1,i1):
 gu:=gu union {seq([i1,op(gu1[j1])],j1=1..nops(gu1))}:
od:
gu:
end:

#SumF1(F,n,kList,kLimits,n0): The sum of F(n,[kList]) at n=n0
#given by the limits indicated by kLimits, in that order
#For example, try:
#SumF1((i+j)!*(n-i)!*(n-j)!/i!^2/j!^2/(n-i-j)!^2,n,[i,j],[[0,n],[0,n-i]],5);
SumF1:=proc(F,n,kList,kLimits,n0)
local Domain1,i1,j1:
Domain1:=SumSet(n,kList,kLimits,n0):
add(eval(subs({n=n0,seq(kList[i1]=Domain1[j1][i1],i1=1..nops(kList))},F)),
j1=1..nops(Domain1)):
end:


#SumF(F,n,kList,kLimits,N0): The first N0+1 terms
#of sum of F(n,[kList]) at n=n0
#given by the limits indicated by kLimits, in that order
#For example, try:
#SumF((i+j)!*(n-i)!*(n-j)!/i!^2/j!^2/(n-i-j)!^2,n,[i,j],[[0,n],[0,n-i]],5);
SumF:=proc(F,n,kList,kLimits,N0) local n0:
[seq(SumF1(F,n,kList,kLimits,n0),n0=1..N0)]:
end:


###findrec (from SCHUTZENBERGER)
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,a,i,j,n0,kv,var1,eq1,mu:
if (1+DEGREE)*(1+ORDER)+2+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:
ope:=subs(var1,ope):
 
 
if nops(kv)>1 then
print(`Too much slack`):
fi:
 
#for i from 1 to nops(kv) do
#  ope:=subs(op(i,kv)=1,ope):
#od:

if ope<>0 then
RETURN(yafe(ope,N)[2]):
else
 RETURN(0):
fi:
 
end:
 

###End findrec (from SCHUTZENBERGER)

#ApplyOpe(ope,n,N,gu): applies the operator ope(n,N) to the sequence gu
ApplyOpe:=proc(ope,n,N,gu)
local i,i1:
[seq(add(subs(n=i,coeff(ope,N,i1))*gu[i+i1],i1=0..degree(ope,N)),
i=1..nops(gu)-degree(ope,N))]:
end:

#OpeToSeq(ope,n,N,Ini,n0,n1): Given an operator ope(n,N) and 
#Initial conditions [a[n0],a[n0+1], ... a[n0+L-1]] (L=degree(ope,N)
#is the order of ope, outputs the sequence up to n=n0+n1-1
OpeToSeq:=proc(ope,n,N,Ini,n0,n1)
local gu,i,ope1,L,j:

ope1:=expand(normal((subs(n=n-ldegree(ope,N),ope)/N^ldegree(ope,N)))):
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`The size of Ini should be`, L , ` the order of`, ope):
fi:

ope1:=normal(ope/coeff(ope,N,L)):
ope1:=expand(normal(subs(n=n-L,ope1)/N^L)):
gu:=Ini:

for i from n0+L to n0+n1-1 do
 gu:=[op(gu),-add(subs(n=i,coeff(ope1,N,-j))*gu[nops(gu)-j+1],j=1..L)]:
od:

gu:

end:


#MulZeilCb(F,kList,kLimits,n,N,para): firsts applies ZeilC, then
#tries to find a bettter operator (i.e. of smaller 
#order and complexity). 
#if succesful, returns the original [ope,cert] followed by
#the better operator. If it couldn't find a better operator
#it returns, [ope,cert],0
#For example, try:
#MulZeilCb(F,n!/i!/j!/(n-i-j)!,[i,j],[[0,n],[0,n-i]],n,N):
MulZeilCb:=proc(F,kList,kLimits,n,N,para) 
local lu,L,Ini,ORDER,DEGREE,Comp1,lu1,
gu,ope,lu2:
lu:=MulZeil(F,kList,n,N,para):

ope:=lu[1]:
L:=degree(ope,N):
Comp1:=(L+1)*(degree(numer(normal(ope)),n)+1):

Ini:=SumF(F,n,kList,kLimits,L):


gu:=OpeToSeq(ope,n,N,Ini,1,(L+2)*(degree(numer(normal(ope)),n)+2)+3):


for ORDER from 0 to L-1 do 
 for DEGREE from 0 to trunc(Comp1/(ORDER+1)) do
  lu1:=findrec(gu,DEGREE,ORDER,n,N):
  if lu1<>0 then
  lu2:=Khalek(ope,lu1,n,N):
  if lu2=FAIL then
   print(`Something is wrong`):
    RETURN(FAIL):
   else
    RETURN(lu1,lu[2],lu2):
  fi:
  fi:
 od:
od:

lu,1:
end:

#Khalek(ope1,ope2,n,N): given recurrence
#operators ope1 and ope2, both monic in
#finds the operators ope3 such that ope1=ope3*ope2
#or returns FAIL
Khalek:=proc(ope1,ope2,n,N) local ope3,d,lu,mu:
if degree(ope1,N)=degree(ope2,N) then
 if ope1=ope2 then
   RETURN(1):
 else
  RETURN(FAIL):
 fi:
fi:

d:=degree(ope1,N)-degree(ope2,N):
 ope3:=expand(ope1-N^d*subs(n=n+d,ope2)):
 lu:=yafe(ope3,N):
  mu:=Khalek(lu[2],ope2,n,N):
  if mu=FAIL then
   RETURN(FAIL):
  fi:

  N^d+lu[1]*mu:
end: