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


print(`First Written: Feb. 2022: tested for Maple 2020 `):
print():
print(`This is PureRec.txt, A Maple package`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(`Linear-Time and Constant-Space Algorithms to compute Multi-Sequences that arise in Enumerative Combinatorics  `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/PureRec.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(`------------------------------------`):

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

 print(`------------------------------------`):




ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(` CheckScheme, CoefD, CoefDc, INIc, INItay,  SchemeD`):

else
ezra(args):

fi:


end:


ezraS:=proc()
if args=NULL then

print(`The STORY procedures are`):
print(` SchemeV`):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` PureRec.txt: A Maple package for computing with pure recurrence relations  using the Apagodu-Zeilberger multivariate extension of the Almkvist-Zeilberger algorithm `):

print(`The MAIN procedures are`):
print(` EvalScheme, Scheme, SchemeC `):


elif nargs=1 and args[1]=CheckScheme then
print(`CheckScheme(R,n,x,K1,K2): Checks Scheme(R,n,x,m,M) for K2 different values in [0,K1]^n. Try:`):
print(`CheckScheme(1/(1-x[1]-x[2]-x[3]-x[1]*x[2]*x[3])^(1/2),3,x,10,10);`):

elif nargs=1 and args[1]=CoefD then
print(`CoefD(F,x,m): Inputs a function of x[1], ..., x[n] (where m:=nops(m)) that has a nice Taylor expansion around the origin, and a list m of n (say) non-neg. integers, outputs`):
print(`the Taylor coefficient (around the origin, in other words, the MacLaurin coefficient) of x[1]^m[1]*...*x[n]^m[n]. Try:`):
print(`CoefD(1/(1-x[1]-x[2]-x[3]-x[4]+x[1]*x[2]*x[3]*x[4]),x,[2,3,4,5]);`):




elif nargs=1 and args[1]=EvalScheme then
print(`EvalScheme(S,m,M,m0): inputs a scheme S (let n:=nops(S)) in n dimension, phrased in terms of m and M (i.e. the discrete variables are m[1],...,m[n] and the corresponding shift operators are`):
print(`M[1],..., M[n], and a list of non-negative integers m0. Evaluats it at m0. Try:`):

print(`EvalScheme([ [2/M[1],3/M[2],4/M[3]],[[[1]]],m,M,[2,2,2]);`):

print(`EvalScheme([[1/M[1]+1/M[1]^2],[1,1]],m,M,[3]);`):

elif nargs=1 and args[1]=INIc then
print(`INIc(P,x,L): Given a Laurent polynomial in the list of variables x, (of size k, say) and a list L of size k+1, outputs`):
print(`a list-of-lists-....L[m1][m2]..[m[k+1]] is the coefficient of x[1]^m[1]*...*x[k]^m[k] in P^m[k+1]`):
print(`The initial conditions up. Try: `):
print(`INIc(x[1]+1+x[1],[x[1]],[2,3]);`);

elif nargs=1 and args[1]=INItay then
print(`INItay(R,x,A): Given a function R in the list of variables x[1],..., x[n] (for some n:=nops(x)) finds the list of lists...of lists, let's call it L, such that for each 0<=a[1]<A[1], ..., 0<=a[n<A[n]`):
print(`The coefficient of x[1]^a[1]*...*x[n]^a[n] in the Taylor expansion is L[a[1]+1]...[a[n]+1]. Try:`):
print(`INItay(1/(1-x-y-z+x^2*y^2+z^2),[x,y,z],[3,3,3]);`):


elif nargs=1 and args[1]=Scheme then
print(`Scheme(R,k,x,m,M): Let k be a positive integer (the second input parameter) and let R be a function of x[1],..., x[k]`):
print(`that has a "nice" Taylor expansion around the origin in n-space, uses the Apagodu-Zeilberger multivariable extension of`):
print(`the Almkvist-Zeilberger algorithm to output a PURE scheme phrased in terms of the discrete variables m[1],...,m[k] `):
print(`and the respective shift operators M[1],...,M[k]. Try:`):
print(`Scheme(1/(1-x[1]-x[2]-x[3]-x[1]*x[2]*x[3])^(1/2),3,x,m,M);`):


elif nargs=1 and args[1]=SchemeC then
print(`SchemeC(P,k,x,m,M,n,N): Let k be a positive integer (the second input parameter) and let P be a Laurent polynomial in x[1],..., x[k]`):
print(`uses the Almkvist-Zeilberger algorithm to output a PURE scheme phrased in terms of the discrete variables m[1],...,m[k], m[k+1]`):
print(`and the respective shift operators M[1],...,M[k], M[k+1], for the discrete function F(m[1],..., m[k+1]):=`):
print(`The coeff. of x[1]^m[1]*...*x[k]^m[k] in P^m[k+1]. Try:`):
print(`SchemeC(x[1]+1/x[1]+1+x[2]+1/x[2],2,x,m,M);`):



elif nargs=1 and args[1]=SchemeD then
print(`SchemeD(R,n,x,m,M): Let n be a positive integer (the second input parameter) and let R be a function of x[1],..., x[n]`):
print(`that has a "nice" Taylor expansion around the origin in n-space, uses the Apagodu-Zeilberger multivariable extension of`):
print(`the Almkvist-Zeilberger algorithm to output a PURE recurrence in terms of the discrete variables m`):
print(`and the  shift operator M for the coefficient of (x[1]*...*x[n])^m`):
print(`SchemeD(1/(1-x[1]-x[2]-x[3]-4*x[1]*x[2]*x[3]),3,x,m,M);`):

elif nargs=1 and args[1]=SchemeV then
print(`SchemeV(R,n,x): A verbose version of Scheme(R,n,x,m,M) with illustration. try:`):
print(`SchemeV(1/(1-x[1]-x[2]-x[1]*x[2])^(1/3),2,x);`):

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

fi:


end:



#START FROM MultiAlmkvistZeilberger

 
ezraMAZ:=proc()
if nargs=0 then
print(`MultiAlmkvistZeilberger`):
print(`A Maple package for  finding recurrences for multi-integrals of `):
print(` Hypergeometric/Hyperexponential Integrands.`):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  MAZ, MAZpaper `):
 
elif nargs=1 and args[1]=MAZ then
print(`MAZ(POL,H,rat,x,n,N,pars): Given a polynomial, POL, an hyper-exponential function, H, and`):
print(`a rational function, rat, of the continuous variables in the list x, and a discrete variable`):
print(`n, and the shift-operator in n, N, and a set of auxiliary parameters, par`):
print(`ouputs a recurrence operator, let's call it  ope(n,N), and a multi-certificate, let's call it R`):
print(`such that the function F(n;x):=POL*H*rat^n satisfies`):
print(` ope(n,N)F=sum(D_{x[i]} (R[i]H*rat^n),i=1..nops(R)`):
print(`In particular, try:`):
print(`MAZ(1,((1-1/x)*(1-y))^b*((1-1/x/y)*(1-1/y))^c/x/y,(1-x)*(1-x*y),[x,y],a,A,{b,c});`):
print(`MAZ(1,exp(-x^2/2-y^2/2),(x-y)^2,[x,y],n,N,{});`):



elif nargs=1 and args[1]=MAZpaper then
print(`MAZpaper(POL,H,rat,x,n,pars,A,R): Given a polynomial, POL, an hyper-exponential function, H, and`):
print(`a rational function, rat, of the continuous variables in the list x, and a discrete variable`):
print(`n, and the shift-operator in n, N, and a set of auxiliary parameters, par`):
print(`and letters A and R for expressing the statement of the theorem`):
print(`and proof`):
print(`ouputs a paper stating and proving `): 
print(`A recurrence equation for the (contour or vanishing-at-endpoints) integral`):
print(`of F:=POL*H*rat^n `):
print(`with respect to the variables in the list [x] `):
print(`Try:`):
print(`MAZpaper(1,((1-1/x)*(1-y))^b*((1-1/x/y)*(1-1/y))^c/x/y,(1-x)*(1-x*y),[x,y],a,{b,c},A,R);`):
print(`MAZpaper(1,exp(-x^2/2-y^2/2),(x-y)^2,[x,y],n,{},A,R);`):


else print(`There is no help for`, args):

fi:



end:





_EnvAllSolutions:=true:

####From MultiZeilberger
#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:

yafe:=proc(ope,N) local i:
add(factor(coeff(ope,N,i))*N^i,i=ldegree(ope,N)..degree(ope,N)):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:

###########End from MultiZeilberger


FindDeg:=proc(POL,H,rat,x,xSet,n,L)
local s,t,h,e,i,Hbar,q,r,gu:
s:=numer(rat): t:=denom(rat):

h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`):

Hbar:=H*s^n/t^(n+L):

gu:=normal(simplify(diff(Hbar,x)/Hbar)):
q:=numer(gu): r:=denom(gu):
max(degree(h,xSet)-degree( diff(r,x)+q,xSet)+degree(r,xSet) ,degree(h,xSet)-degree(r,xSet)+1+degree(r,xSet)):
end:

#MAZ1(POL,H,rat,x,n,N,L,pars): Given a polynomial POL of the discrete variable n 
#and the continuous variables in x
#and a pure-hyperexponential function H of the variables in x, and a rational function
#Rat of x, and a shift operator N, and a non-neg. integer L
#returns FAIL or, if there is one,  an operator ope(N,n) of order L and
#a list of rational functions R, such that F:=POL*H*rat^n satisfies
#ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) .
#For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{})
MAZ1:=proc(POL,H,rat,x,n,N,L,pars)
local deg,deg1,m,gu,i,i1:
deg:=[]:
for i from 1 to nops(x) do
deg:=[op(deg),FindDeg(POL,H,rat,x[i],{seq(x[i1],i1=1..nops(x))},n,L)]:
od:

m:=min(op(deg)):

for i from 0 to m do
 deg1:=[seq(deg[i1]-m+i,i1=1..nops(deg))]:

 gu:=MAZ1deg(POL,H,rat,x,n,N,L,pars,deg1):

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


MAZ:=proc(POL,H,rat,x,n,N,pars) local gu,L:

for L from 0 do
gu:=MAZ1(POL,H,rat,x,n,N,L,pars):

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


CheckMAZ:=proc(POL,H,rat,x,n,N,pars)
local F,gu,luL,luR,R,ope,i:
F:=H*rat^n:
gu:=MAZ(POL,H,rat,x,n,N,pars):
if gu=FAIL then
RETURN(FAIL):
fi:
ope:=gu[1]:
R:=gu[2]:
luL:=normal(add(subs(n=n+i,POL)*coeff(ope,N,i)*normal(simplify(subs(n=n+i,F)/F)),i=0..degree(ope,N))):
luR:=normal(add( normal(diff(R[i]*F,x[i])/F ),
              i=1..nops(R))):
evalb(normal(luL/luR)=1):
end:



#MAZ1deg(POL,H,rat,x,n,N,L,pars,deg): Given a polynomial POL of the discrete variable n 
#and the continuous variables in x
#and a pure-hyperexponential function H of the variables in x, and a rational function
#Rat of x, and a shift operator N, and a non-neg. integer L
#returns FAIL or, if there is one,  an operator ope(N,n) of order L and
#a list of rational functions R, such that F:=POL*H*rat^n satisfies
#ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) .
#For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{})
MAZ1deg:=proc(POL,H,rat,x,n,N,L,pars,deg)
local s,t,h,e,i,Hbar,gu,X,a,var,q,r,eq,ope,var1,ku,i1,eqN,var1N,opeN,bilti,meka:

s:=numer(rat): t:=denom(rat):
ope:=add(e[i]*N^i,i=0..L):
h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`):

Hbar:=H*s^n/t^(n+L):
gu:=[]:

for i from 1 to nops(x) do
 gu:=[op(gu),normal(simplify(diff(Hbar,x[i])/Hbar))]:
od:

q:=[]: r:=[]:
for i from 1 to nops(gu) do
q:=[op(q),numer(gu[i])]: 
r:=[op(r),denom(gu[i])]:
od:

var:={}:
X:=[]:
for i from 1 to nops(x) do
 ku:=GenPol(x,a[i],deg[i]):
 X:=[op(X),ku[1]]:
 var:=var union ku[2]:
od:

 var:=var union {seq(e[i],i=0..L)}:

gu:=h:

for i from 1 to nops(x) do
gu:=expand(normal(gu-(diff(r[i],x[i])+q[i])*X[i]/r[i]-r[i]*diff(X[i]/r[i],x[i]))):
od:
eq:=Extract1(numer(gu),x):

eqN:=subs(n=9/17,eq):
eqN:=subs({seq(pars[i]=1/(i^2+3),i=1..nops(pars))},eqN):

var1N:=solve(eqN,var):

opeN:=subs(var1N,ope):

if opeN=0 then
  RETURN(FAIL):
fi:



var1:=solve(eq,var):

ope:=subs(var1,ope):
X:=subs(var1,X):

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


bilti:={}:

for i from 1 to nops(var1)  do

 if op(1,op(i,var1))=op(2,op(i,var1)) then
   bilti:=bilti union {op(1,op(i,var1))}:
 fi:
  
od:




for i from 1 to nops(bilti) do
 if coeff(ope,bilti[i],1)<>0 then
    ope:=coeff(ope,bilti[i],1):
    X:=[seq(coeff(X[i1],bilti[i],1),i1=1..nops(X))]:
   
    meka:=coeff(ope,N,degree(ope,N)):
     ope:=expand(normal(ope/meka)):
     ope:=yafe(ope,N):
     X:=[seq(normal(X[i1]/meka/t^L),i1=1..nops(X))]: 
    RETURN(ope,X):
 fi:
od:

end:



#MAZ(POL,H,rat,x,n,pars,F,A,R)
MAZpaper:=proc(POL,H,rat,x,n,pars,A,R) local gu,F,F1,ope,i,N:

print(``):

F1:=H*rat^n:

gu:=MAZ(POL,H,rat,x,n,N,pars):
ope:=gu[1]:
if gu=FAIL then
 RETURN(FAIL):
fi:
print(`A Linear Recurrence Equation for an Integral in`, nops(x), ` variables `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let `):
print(F(n,op(x))=F1):
print(`and let`, A(n), `be its definite integral  of`):
print(POL*F(n,op(x))):
print(`with respect to the continuous`):
print(` variables `, op(x)):
print(``):
print(A(n) , `satisfies the following homogeneous linear recurrence with`):
print(`polynomial coefficients: `):
print(``):
print(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0):
print(``):
print(`Proof: We cleverly constuct the rational functions`):

for i from 1 to nops(gu[2]) do
print(R[i](n,op(x))=gu[2][i]):
od:
print(`With the motive that`):
print(add(coeff(ope,N,i)*F(n+i,op(x)),i=0..degree(ope,N))=
add(D[x[i]](R[i](n,op(x))*F(n,op(x))),i=1..nops(x))):
print(``):
print(`(Check!)`):
print(`and the theorem follows upon integrating with respect to`, op(x)):
print(` QED! `):



end:

###END FROM MultiAlmkvistZeilberger




#CoefD(F,x,m): Inputs a function of x[1], ..., x[n] (where m:=nops(m)) that has a nice Taylor expansion around the origin, and a list m of n (say) non-neg. integers, outputs
#the Taylor coefficient (around the origin, in other words, the MacLaurin coefficient) of x[1]^m[1]*...*x[n]^m[n]. Try:
#CoefD(1/(1-x[1]-x[2]-x[3]-x[4]+x[1]*x[2]*x[3]*x[4]),x,[2,3,4,5]);

CoefD:=proc(F,x,m) local n,i,F1:

n:=nops(m):


F1:=F:
for i from 1 to n do
F1:=taylor(F1,x[i]=0,m[i]+2):
F1:=expand(normal(coeff(F1,x[i],m[i]))):
od:
F1:

end:




#INItay(R,x,A): Given a function R in the list of variables x[1],..., x[n] (for some n:=nops(x)) finds the list of lists...of lists, let's call it L, such that for each 0<=a[1]<A[1], ..., 0<=a[n<A[n]
#The coefficient of x[1]^a[1]*...*x[n]^a[n] in the Taylor expansion is L[a[1]+1]...[a[n]+1]. Try:
#INItay(1/(1-x-y-z+x^2*y^2+z^2),[x,y,z],[3,3,3]);
INItay:=proc(R,x,A) local n,i,R1:
n:=nops(x):
if not (type(x,list) and nops(A)=n) then
 RETURN(FAIL):
fi:

R1:=taylor(R,x[1]=0,A[1]):
if nops(x)=1 then
 RETURN([seq(expand(coeff(R1,x[1],i)),i=0..A[1]-1)]):
else
[seq(INItay(coeff(R1,x[1],i),[op(2..n,x)],[op(2..n,A)]),i=0..A[1]-1)]:
fi:

end:

#Scheme(R,n,x,m,M): Let n be a positive integer (the second input parameter) and let R be a function of x[1],..., x[n]
#that has a "nice" Taylor expansion around the origin in n-space, uses the Apagodu-Zeilberger multivariable extension of
#the Almkvist-Zeilberger algorithm to output a PURE scheme phrased in terms of the discrete variables m[1],...,m[n] 
#and the respective shift operators M[1],...,M[n]. Try:
#Scheme(1/(1-x[1]-x[2]-x[3]-x[1]*x[2]*x[3])^(1/2),3,x,m,M);
Scheme:=proc(R,n,x,m,M) local  i,j,lu,GU,INI:

GU:=[]:

for i from 1 to n do
lu:=MAZ(1,R/(mul(x[j]^(m[j]+1),j=1..i-1)*mul(x[j]^(m[j]+1),j=i+1..n)*x[i]),1/x[i],[seq(x[j],j=1..n)],m[i],M[i],{seq(m[j],j=1..i-1), seq(m[j],j=i+1..n)})[1]:
lu:=expand(1-subs(m[i]=m[i]-degree(lu,M[i]),lu)/M[i]^degree(lu,M[i])):
lu:=add(factor(coeff(lu,M[i],-j))*M[i]^(-j),j=1..-ldegree(lu,M[i])):
GU:=[op(GU),lu]:
od:

INI:=INItay(R,[seq(x[i],i=1..n)],[seq(-ldegree(GU[i],M[i]),i=1..n)]):
[GU,INI]:

end:



EvalScheme1:=proc(S,m,M,m0) local ope,INI,L,a,kha,j:
ope:=S[1][1]:
INI:=S[2]:
L:=INI:

if m0<=nops(INI)-1 then
 RETURN(INI[m0+1]):
fi:

for a from nops(INI) to m0 do

if member(0,
{seq(subs(m=a,denom(coeff(ope,M,-j))),j=1..-ldegree(ope,M))}
) then
print(a, `is a singularity for`, ope):
 RETURN(FAIL):
else

kha:=normal(add(subs(m=a,coeff(ope,M,-j))*L[-j],j=1..-ldegree(ope,M))):
L:=[op(2..nops(L),L),kha]:
fi:
od:

L[-1]:

end:


#EvalScheme(S,m,M,m0): inputs a scheme S (let n:=nops(S)) in n dimension, phrased in terms of m and M (i.e. the discrete variables are m[1],...,m[n] and the corresponding shift operators are
#M[1],..., M[n], and a list of non-negative integers m0. Evaluats it at m0. Try
#EvalScheme([[2/M[1],3/M[2],4/M[3]],[[[1]]],m,M,[2,2,2]);
#EvalScheme([[m[1]/M[1]],[1]],m,M,[3]);

EvalScheme:=proc(S,m,M,m0) local n,ope1,GU,S1,A1,INI1,a,i1,m01,i,ka:

n:=nops(S[1]):
if nops(m0)<>n then
 RETURN(FAIL):
fi:

if n=1 then
RETURN(EvalScheme1(S,m[1],M[1],m0[1])):
fi:


ope1:=S[1][1]:
A1:=-ldegree(ope1,M[1]):




if m0[1]<A1 then
INI1:=S[2][m0[1]+1]:

S1:=[op(2..n,S[1])]:
S1:=subs(m[1]=m0[1],S1):
S1:=subs({seq(m[i1]=m[i1-1],i1=2..n),seq(M[i1]=M[i1-1],i1=2..n)},S1):
m01:=[op(2..n,m0)]:
ka:=EvalScheme([S1,INI1],m,M,m01):
if ka=FAIL then
print(m01 , `is a singularity `):
RETURN(FAIL):
else
RETURN(EvalScheme([S1,INI1],m,M,m01)):
fi:
fi:

GU:=[]:

for a from 0 to A1-1 do

INI1:=S[2][a+1]:


S1:=[op(2..n,S[1])]:
S1:=subs(m[1]=a,S1):
S1:=subs({seq(m[i1]=m[i1-1],i1=2..n),seq(M[i1]=M[i1-1],i1=2..n)},S1):
m01:=[op(2..n,m0)]:

ka:=EvalScheme([S1,INI1],m,M,m01):

if ka=FAIL then
print(m01 , `is a singularity `):
RETURN(FAIL):
else
GU:=[op(GU),ka]:
fi:
od:

ope1:=subs({seq(m[i]=m0[i],i=2..n)},ope1):

EvalScheme1([[ope1],GU],m[1],M[1],m0[1]):
end:



#CheckScheme(R,n,x,K1,K2): Checks SchemeS(R,n,x,m,M) for K2 different values in [0,K1]^n. Try:
#CheckScheme(1/(1-x[1]-x[2]-x[3]-x[1]*x[2]*x[3])^(1/2),3,x,10,10);
CheckScheme:=proc(R,n,x,K1,K2) local  S,m,M,ra,pt,i,i1:
S:=Scheme(R,n,x,m,M):
ra:=rand(0..K1):

for  i from 1 to K2 do

pt:=[seq(ra(),i1=1..n)]:

if EvalScheme(S,m,M,pt)<>CoefD(R,x,pt)  then
 print(`pt failed `):
 print( [EvalScheme(S,m,M,pt),CoefD(R,x,pt)]):
 RETURN(false):
fi:
od:

true:

end:




#SchemeV(R,n,x): A verbose version of Scheme(R,n,x,m,M) with illustration. try:
#SchemeV(1/(1-x[1]-x[2]-x[1]*x[2])^(1/3),2,x);
SchemeV:=proc(R,n,x) local  gu,m,M,F,i,t0,t1, ope,i1,i2,i3,INI,pt:

print(`A Pure Recurrence Scheme that enables Linear-Time and Constant-Space Calculation of the  Taylor coefficients of the function with`, n, `variables `, seq(x[i],i=1..n)):
print(``):
print(R):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let `, F(seq(m[i],i=1..n)), ` be the coefficient of`, mul(x[i]^m[i],i=1..n), ` in the Taylor expansion (around the origing) of the multivariable function (that may also be viewed as a formal power series) of`, seq(x[i],i=1..n) ):
print(``):
print(R):
print(``):

t0:=time():
gu:=Scheme(R,n,x,m,M):
t1:=time()-t0:

print(`The following PURE recurrence relations hold in each of the`, n, `discrete variables`, seq(m[i],i=1..n)):
print(``):

for i from 1 to n do
ope:=gu[1][i]:
print(``):
print(F(seq(m[i1],i1=1..n))= add(coeff(ope,M[i],-i1)*F(seq(m[i2],i2=1..i-1),m[i]-i1,seq(m[i2],i2=i+1..n)),i1=1..-ldegree(ope,M[i]) ) ):
print(``):
od:

print(``):
print(`and in Maple notation `):
print(``):

for i from 1 to n do
ope:=gu[1][i]:
print(``):
lprint(F(seq(m[i1],i1=1..n))= add(coeff(ope,M[i],-i1)*F(seq(m[i2],i2=1..i-1),m[i]-i1,seq(m[i2],i2=i+1..n)),i1=1..-ldegree(ope,M[i]) ) ):
print(``):
od:

print(``):
print(`Subject to the initial conditions`):



INI:=gu[2]:

if n=1 then
 print(seq(F(i1)=INI[i1+1],i1=0..nops(INI)-1)):
elif n=2 then
 print(seq(seq(F(i1,i2)=INI[i1+1][i2+1],i2=0..nops(INI[i1+1])-1),i1=0..nops(INI)-1)):
elif n=3 then
 print(seq(seq(seq(F(i1,i2,i3)=INI[i1+1][i2+1][i3+1],i3=0..nops(INI[i1+1][i2+1])-1),i2=0..nops(INI[i1+1])-1),i1=0..nops(INI)-1)):
else
print(INI):
fi:

print(``):
print(`--------------------------------------------`):
print(``):

print(`By the way, it took`, t1, `seconds to geneate the scheme. `):
print(``):
print(`--------------------------------------------`):
print(``):

t0:=time():
print(`Using this scheme we can compute any value very fast. For example`):
pt:=[seq(1000*i1,i1=1..n)]:
print(``):
print(F(op(pt))=EvalScheme(gu,m,M,pt)):
print(``):
print(`This took`, time()-t0, `seconds. `):
print(``):

end:






#SchemeD(R,n,x,m,M): Let n be a positive integer (the second input parameter) and let R be a function of x[1],..., x[n]
#that has a "nice" Taylor expansion around the origin in n-space, uses the Apagodu-Zeilberger multivariable extension of
#the Almkvist-Zeilberger algorithm to output a PURE scheme phrased in terms of the discrete variables m
#and the  shift operator M for the coefficient of (x[1]*...*x[n])^m
#SchemeD(1/(1-x[1]-x[2]-x[3]-x[1]*x[2]*x[3])^(1/2),3,x,m,M);
SchemeD:=proc(R,n,x,m,M) local  i,j,lu,INI:



lu:=MAZ(1,R/mul(x[j],j=1..n),1/mul(x[j],j=1..n),[seq(x[j],j=1..n)],m,M,{})[1]:
lu:=expand(1-subs(m=m-degree(lu,M),lu)/M^degree(lu,M)):
lu:=add(factor(coeff(lu,M,-j))*M^(-j),j=1..-ldegree(lu,M)):
INI:=[seq(CoefD(R,x,[i$n]),i=0..-ldegree(lu,M)-1)]:
lu:=subs({M=M[1],m=m[1]},lu):
[[lu],INI]:

end: