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


print(`First Written: Nov. 17, 2006: tested for Maple 10 `):
print(`Version of Nov. 17, 2006:  `):
print():
print(`This is OneDimWalks, A Maple package`):
print(`To study walks on the half-line`):
print(`accompanying Manuel Kauers and Doron Zeilberger's  article: `):
print(`"The Quasi-Holonomic Ansatz and Restricted Lattice Walks" `):

print():
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/OneDimWalks .`):
print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
#print(`For a list of the supporting functions type: ezra1();`):
print():
 


ezra1:=proc()
if args=NULL then

print(`The supporting procedures are`):
print(` Aij, ApplyOper1,CheckConjDE, CODVa, CODVm, Comm1, ConjDEX1`):
print(` DiA, DiffToRec1, EvolutionOper, FindHorizRec,`):
print(` findrec, Findrec, FindrecX `):
print(` FindVertRec, FindVertRec1, `):
print(` GR1, HafelDO, HafelOper, Mult1, Pxt, Pxt1,`):
print(` NuWalks, RecToDiff, SeqAnx, SeqFromRec, SeqX, Walks  `):
else
 ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then

print(` OneDimWalks: A Maple package for  `):

print(`The MAIN procedures are`):
print(` ConjDE, ConjRec, ConjRecX, `):
print(` FindGoodmRec, FindGoodmRec1, FunctOper, Seq `):


elif nargs=1 and args[1]=Aij then
print(`Aij(t,Dt,i,j): Expanding (tDt)^i(t^(-j)) as a differential operator`):
print(`in the standard form. For example, try: Aij(t,Dt,1,1);`):


elif nargs=1 and args[1]=ApplyOper1 then
print(`ApplyOper1(ope,n,N,x,X,L): Given an operator ope(n,N,x,X) `):
print(`where n, x are variables, N is the forward shift operator`):
print(`in n and X is evaluation at x=0, and a list of`):
print(`polynomials [a(1,x),a(2,x), ...] applies ope to the`):
print(`list. For example do: `):
print(`ApplyOper1(N+x,n,N,x,X,[x,x^2,x^3]);`):

elif nargs=1 and args[1]=CheckConjDE then
print(`CheckConjDE(Steps,x0,MaxC,t,Dt,K): checks`):
print(`ConjDE(Steps,x0,MaxC,t,Dt) up to K-th term`):
print(`For example, try: `):
print(`CheckConjDE({1,-1},2,6,t,Dt,20);`):

elif nargs=1 and args[1]=CODVa then
print(`CODVa(DE1,A,t,Dt): Given a differential equation DE1=[DO,rhs]`):
print(`phrased in terms of t and Dt, and a rational function, A, of t`):
print(`finds the differential equation satisfied by Z(t):=Y(t)-A(t)`):
print(`if Y(t) is a solution of DE1 (i.e. DO (Y(t))=rhs(t) `):
print(`For example, try: CODVa([Dt,0],t^2,t,Dt);`):

elif nargs=1 and args[1]=CODVm then
print(`CODVm(DE1,A,t,Dt): Given a differential equation DE1=[DO,rhs]`):
print(`phrased in terms of t and Dt, and a rational function, A, of t`):
print(`finds the differential equation satisfied by Z(t):=Y(t)/A(t)`):
print(`if Y(t) is a solution of DE1 (i.e. DO (Y(t))=rhs(t) `):
print(`For example, try: CODVm([Dt,0],t^2,t,Dt);`):

elif nargs=1 and args[1]=Comm1 then
print(`Comm1(ope1,ope2,n,N,x,X): The commutator of ope1 and ope2`):


elif nargs=1 and args[1]=ConjDE then
print(`ConjDE(Steps,x0,MaxC,t,Dt):`):
print(`Given a set of steps (which are integers)`):
print(`and a positive integer, conjectures a homog. linear differential`):
print(`equation satisfied by the generating function  (w.r.t. t) of`):
print(`the sequence:`):
print(`The number of walks with n steps, staying at the half-line x>=0`):
print(`starting at the origin and winding-up and x0.`):
print(`Dt is the differentiation operator w.r.t. t.`):
print(`The output is a pair [ope, rhs], where ope`):
print(`is a linear differential operator with polynomial coefficients`):
print(`and rhs is a polynomial in t, and it means that if f(t)`):
print(`is the above-mentioned generating function, then`):
print(` ope(t,Dt)f(t)=rhs `):
print(` For example, try: ConjDE({1,-1},2,6,t,Dt); `):

elif nargs=1 and args[1]=ConjDEX1 then
print(`ConjDEX1(Steps,MaxC,t,Dt,x): Conjecturs a`):
print(`diferential equation satisfies by F(x,t, Steps)`):
print(`using CODV etc. this program only allows one negative`):
print(`step: -1 .`):
print(`For example, try:`):
print(`ConjDEX1({1,-1},4,t,Dt,x);`):


elif nargs=1 and args[1]=ConjRec then
print(`ConjRec(Steps,x0,n,N,MaxC): Given a set of steps, Steps `):
print(` (which are integers)`):
print(`and a positive integer x0, conjectures a homog. linear recurrence`):
print(`operator, with polynomial coefficients, ope(n,N)`):
print(`of maximum complexity MaxC `):
print(`(where N is the shift operator in n), for the sequence`):
print(`The number of walks with n steps, staying at the half-line x>=0`):
print(`starting at the origin and winding-up and x0.`):
print(`using the steps of the set Steps`):
print(`For example, try: ConjRec({1,-1},2,n,N,6);`):

elif nargs=1 and args[1]=ConjRecX then
print(`ConjRecX(Steps,x,n,N,MaxC): Given a set of steps, Steps `):
print(` (which are integers)`):
print(`and a positive integer x0, conjectures a homog. linear recurrence`):
print(`operator, with polynomial coefficients, ope(n,N)`):
print(`of maximum complexity MaxC `):
print(`(where N is the shift operator in n), for the sequence`):
print(`The  weight-enumerator according to x^(destination)`):
print(` for the walks with n steps, staying at the half-line x>=0`):
print(`using the steps of the set Steps`):
print(` The output is a pair [ope, IniConditions] `):
print(`For example, try: ConjRecX({1,-1},x,n,N,4); `):



elif nargs=1 and args[1]=FindGoodmRec then
print(`FindGoodmRec(Steps,m,n,N, M, degn,degm,degN,degM): `):
print(`Given a set of steps Steps `):
print(`(a set of integers) and`):
print(`discrete variables m and n, the symbol for the shift operator`):
print(` in n, N,`):
print(`and postitive integers degm, degn , degm, degN`):
print(`conjectures a recurrence operator ope(m,n,N),`):
print(`of the form P(n,N)+mQ(n,m,M,N) `):
print(` of degree <=degm in m, degree <=degn in n, degree <=degN in N`):
print(`annilihating F(m,n):= the number of ways`):
print(`of walking n steps on the positive discrete half-line, starting at`):
print(`0, and ending in m. For example, try:`):
print(`FindGoodmRec({1,-1},m,n,N,M,2,2,2,2);`):

elif nargs=1 and args[1]=FindGoodmRec1 then
print(`FindGoodmRec1(Steps,m,n,N, M, degn,degm,degN,degM, MaxC): `):
print(`Given a set of steps Steps `):
print(`(a set of integers) and`):
print(`discrete variables m and n, the symbol for the shift operator`):
print(` in n, N,`):
print(`and postitive integers degm, degn , degm, degN, MaxC`):
print(`conjectures a recurrence operator ope(m,n,N),`):
print(`of the form P(n,N)+mQ(n,m,M,N) `):
print(`where P(n,N) has complexity <=MaxC , and `):
print(` of degree <=degm in m, degree <=degn in n, degree <=degN in N`):
print(`annilihating F(m,n):= the number of ways`):
print(`of walking n steps on the positive discrete half-line, starting at`):
print(`0, and ending in m. For example, try:`):
print(`FindGoodmRec1({1,-1},m,n,N,M,2,2,2,2,4);`):

elif nargs=1 and args[1]=DiA then
print(`DiA(A,t,Dt,i): The differential operator Dt^i(A(t))`):
print(`For example, try: DiA(t^4,t,Dt,3);`):

elif nargs=1 and args[1]=DiffToRec1 then
print(`DiffToRec1(ope,t,Dt,n,N): Given a linear differential operator`):
print(`ope(t,Dt) and symbols n,N, outputs the linear recurrence`):
print(`operator with polynomial coefficients satisfied by`):
print(`the sequence of coeff. of any f.p.s. annihilating ope(t,Dt)`):
print(`where N is the shift in n. For example, try:`):
print(`DiffToRec1(Dt+t, t,Dt,n,N);`):

elif nargs=1 and args[1]=EvolutionOper then
print(`EvolutionOper(x,t,Steps,X): The "evolution operator" corresponding`):
print(`to the set of steps Steps, phrased in terms of`):
print(`a space variable x and a time variable t and where`):
print(`X is the operator such that X^j F(x,t) means`):
print(`coeff of x^(j-1) in F(x,t).  For example: try:`):
print(`EvolutionOper(x,t,{1,-1},X);`):



elif nargs=1 and args[1]=FindHorizRec then
print(`FindHorizRec(Steps,m,n,M, degn,degm,degM): `):
print(` Given a set of steps Steps `):
print(`(a set of integers) and`):
print(`discrete variables m and n, the symbol for the shift operator`):
print(` in m, M,`):
print(`and postitive integers degm, degn, degM`):
print(`conjectures a recurrence operator ope(m,n,M),`):
print(`of degree <=degm in m, degree<=degn in n, degree degM in M,`):
print(`annilihating F(m,n):= the numbr of ways`):
print(`of walking n steps on the positive discrete half-line, starting at`):
print(`0, and ending in m. For example, try:`):
print(`FindHorizRec({1,-1},m,n,N,2,2,2);`);

elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER.`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nargs=1 and args[1]=FindrecX then
print(`FindrecX(f,MaxC,x,n,N,DegX): `):
print(`guesses a recurrence operator of complexity <=MaxC,`):
print(`in n, and DegX in x`):
print(`annihilating the sequence f of polynomials in x .`):
print(`For example, try: FindrecX([seq(x^i,i=1..20)],2,x,n,N,2);`):


elif nargs=1 and args[1]=FindVertRec then
print(`FindVertRec(Steps,m,n,N, degn,degm,degN): Given a set of steps Steps `):
print(`(a set of integers) and`):
print(`discrete variables m and n, the symbol for the shift operator`):
print(` in n, N,`):
print(`and postitive integers degm, degn, degN`):
print(`conjectures a recurrence operator ope(m,n,N),`):
print(`of degree <=degm in m, degree<=degn in n, degree degN in N,`):
print(`annilihating F(m,n):= the numbr of ways`):
print(`of walking n steps on the positive discrete half-line, starting at`):
print(`0, and ending in m. For example, try:`):
print(`FindVertRec({1,-1},m,n,N,2,2,2);`);


elif nargs=1 and args[1]=FindVertRec1 then
print(`FindVertRec1(Steps,m,n,N, maxC, degm): `):
print(`Given a set of steps Steps `):
print(`(a set of integers) and`):
print(`discrete variables m and n, the symbol for the shift operator in n,`):
print(`and postitive integers MaxC and degm, `):
print(`conjectures a recurrence operator ope(m,n,N),`):
print(`of complexity MaxC in (n,N) and degree <=degm in m`):
print(`annilihating F(m,n):= the numbr of ways`):
print(`of walking n steps on the positive discrete half-line, starting at`):
print(`0, and ending in m. For example, try:`):
print(`FindVertRec1({1,-1},m,n,N,4,4);`):

elif nargs=1 and args[1]=FunctOper then
print(`FunctOper(Steps,N,x,X): Given a set of steps Steps, and`):
print(`symbols N,x,X, outputs the function operator`):
print(`P(N,x,X) annihilating a(n,x): the polynomial-generating function`):
print(`(according to final destination) for n-step one-dimensional walks`):
print(`using the set of Steps (a set of intgers)`):
print(`Here N is shift  in the n direction and X[j]f(x) is`):
print(`coeff. of x^j . For example, try:`):
print(` FunctOper({-1,1},N,x,X); `):

elif nargs=1 and args[1]=GR1 then
print(`GR1(S1,x): guesses a polynomial from the data-set S1`):
print(`For example, try:`):
print(` GR1({seq([i,(i+4)/(i+9)],i=1..10)},x);`):

elif nargs=1 and args[1]=HafelDO then
print(`HafelDO(ope,t,Dt,F): applies the differential operator ope(t,Dt)`):
print(`to the function F(t). For example, try:`):
print(`HafelDO(Dt,t,Dt,t);`):

elif nargs=1 and args[1]=HafelOper then
print(`HafelOper(oper1,F1,x,X): Given an operator oper1 in the`):
print(`form [function, Operator] applies it to the polynomial F1`):
print(`in x and t, where X^j F(x) is the coeff. of x^(j-1) in F`):
print(`For example, try: HafelOper([1,t*x+t/x*(1-X)],1,x,X);`):


elif nargs=1 and args[1]=Mult1 then
print(`Mult1(ope1,ope2,n,N,x,X): ope1(n,N,x,X)ope2(n,N,x)`):
print(`where X is the opeartor X f(x):=f(0)`):
print(`for example, try:`):
print(`Mult1(x*N-1-x^2+X,N^2-x*N+1,n,N,x,X);`):

elif nargs=1 and args[1]=NuWalks then
print(`NuWalks(Steps,t,x): Given a set of steps Steps (positive or negative `):
print(`integers or 0) and numbers t and x returns the number of walks`):
print(`(Sequences in Steps) after time t that sum to x and never`):
print(`visited negative sites`):
print(`For example, try: NuWalks({1,-1},6,0);`):

elif nargs=1 and args[1]=Pxt then
print(`Pxt(Steps,x,t,N): The beginning of the expansion (up to power N)`):
print(`of F(Steps,x,t). For example, do:`):
print(`Pxt({1,-1},x,t,10);`):

elif nargs=1 and args[1]=Pxt1 then
print(`Pxt1(Steps,x,t,N): The beginning of the expansion (up to power N)`):
print(`of F(Steps,x,t), using EvolutionOper. `):
print(`It should give the same output as Pxt(Steps,x,t,N) `):
print(` For example, do:`):
print(`Pxt1({1,-1},x,t,10);`):


elif nargs=1 and args[1]=RecToDiff then
print(`RecToDiff(ope,n,N,t,Dt,Start1): turns a linear recurrence operator`):
print(`with polynomial coefficients ope(n,N), with initial conditions`):
print(`into the`):
print(`differential operator ope(tDt,1/t) . For example try:`):
print(`RecToDiff(N-n-1,n,N,t,Dt,[1,1,2,6,24,120]);`):

elif nargs=1 and args[1]=Seq then
print(`Seq(Steps,N,x): the sequence of NuWalks ending at x with`):
print(`set of steps Steps`):


elif nargs=1 and args[1]=SeqAnx then
print(`SeqAnx(Steps,x,N):`):
print(`The sequence of length N, whose i^th entry is the polynomial`):
print(`that is the weight-enumerator of all paths starting at the`):
print(`origin, staying in the non-negative half-line, and using`):
print(`the set of steps, weighted by x^(final destination)`):
print(`For example, try: SeqAnx({1,-1},x,10);`):

elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):

elif nargs=1 and args[1]=SeqX then
print(`SeqX(Steps,x,K): The first K terms in the weight-enumerating`):
print(`sequence for walks in the positive half-line using the steps in`):
print(`in Steps, For example, try: `):
print(`SeqX({1,-1},x,10);`):

elif nargs=1 and args[1]=Walks then
print(`Walks(Steps,t,x): Given a set of steps Steps (positive or negative `):
print(`integers or 0) and numbers t and x returns the set of walks`):
print(`(Sequences in Steps) after time t that sum to x and never`):
print(`visited negative sites`):
print(`For example, try: Walks({1,-1},6,0);`):

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

fi:

end:


#Walks(Steps,t,x): Given a set of steps Steps (positive or negative 
#integers or 0) and numbers t and x returns the set of walks
#(Sequences in Steps) after time t that sum to x and never
#visited negative sites
#For example, try: Walks({1,-1},5,0);
Walks:=proc(Steps,t,x)
local gu,i,gu1,s:
option remember:

if t=0 then
 if x=0 then
  RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:

if x<0 then
 RETURN({}):
fi:

gu:={}:
for s in Steps do
 gu1:=Walks(Steps,t-1,x-s):
  gu:=gu union {seq([op(gu1[i]),s],i=1..nops(gu1))}:
od:
gu:
end:

#NuWalks(Steps,t,x): Given a set of steps Steps (positive or negative 
#integers or 0) and numbers t and x returns the number of walks
#(Sequences in Steps) after time t that sum to x and never
#visited negative sites
#For example, try: NuWalks({1,-1},5,0);
NuWalks:=proc(Steps,t,x)
local s:
option remember:

if t=0 then
 if x=0 then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:

if x<0 then
 RETURN(0):
fi:

add(NuWalks(Steps,t-1,x-s), s in Steps):


end:


#Seq(Steps,N,x): the sequence of NuWalks ending at x with
#set of steps Steps
Seq:=proc(Steps,N,x) local t:
[seq(NuWalks(Steps,t,x),t=0..N)]:
end:

Squeeze:=proc(L) local i,L1:
L1:=[]:
for i from 1 to nops(L) do
 if L[i]<>0 then
  L1:=[op(L1),L[i]]:
 fi:
od:
L1:
end:

#Pxt(Steps,x,t,N): The beginning of the expansion (up to power N in t)
#of F(Steps,x,t). For example, do:
#Pxt({1,-1},x,t,10)
Pxt:=proc(Steps,x,t,N) local i,j:
add(add(NuWalks(Steps,i,j)*t^i*x^j,j=0..N*max(op(Steps))),i=0..N):
end:

################Start Findrec and its offsprings

with(linalg):
#with(SolveTools);
#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+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,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:


#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+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 ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

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:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 1 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:

################End Findrec and its offsprings







#ConjRec(Steps,x0,n,N,MaxC): Given a set of steps (which are integers)
#and a positive integer, conjectures a homog. linear recurrence
#operator, with polynomial coefficients, ope(n,N)
#(where N is the shift operator in n), for the sequence
#The number of walks with n steps, staying at the half-line x>=0
#starting at the origin and winding-up and x0.
#For example, try: ConjRec({1,-1},2,n,N,6):
ConjRec:=proc(Steps,x0,n,N,MaxC) local gu,N0,ope,i,i0:
option remember:
N0:=trunc((4+(MaxC/2))*(2+MaxC/2)+10+x0):
gu:=Seq(Steps,N0,x0):
gu:=[op(2..nops(gu),gu)]:


for i from 1 while gu[i]=0 do od:

i0:=i:

gu:=[op(i0+1..nops(gu),gu)]:
ope:=Findrec(gu,n,N,MaxC):
ope:=Yafe(subs(n=n-i0,ope),N)[2]:

ope:
end:


#EvolutionOper(x,t,Steps,X): The "evolution operator" corresponding
#to the set of steps Steps, phrased in terms of
#a space variable x and a time variable t and where
#X is the operator such that X^j F(x,t) means
#coeff of x^(j-1) in F(x,t).  For example: try:
#EvolutionOper(x,t,{1,-1},X);
EvolutionOper:=proc(x,t,Steps,X) local gu,s,j:
gu:=0:

for s in Steps do
 if s>=0 then
 gu:=gu+t*x^s:
 else
 gu:=gu+t*x^s*(1-add(X^j,j=1..-s)):
 fi:
od:

[1,gu]:
end:


#HafelOper(oper1,F1,x,X): Given an operator oper1 in the
#form [function, Operator] applies it to the polynomial F1
#in x and t, where X^j F(x) is the coeff. of x^(j-1) in F
#For example, try: OneStep([1,t*x+t/x*(1-X)],1,x,X);
HafelOper:=proc(oper1,F1,x,X)
local gu,ope,i:

ope:=expand(oper1[2]):

gu:=expand(oper1[1]+coeff(ope,X,0)*F1):

for i from 1 to degree(ope,X) do
 gu:=expand(gu+coeff(ope,X,i)*coeff(F1,x,i-1)*x^(i-1)):
od:

gu:

end:


#Pxt1(Steps,x,t,N): The same output as Pxt(Steps,x,t,N)
#but by applying EvolutionOper. For example, try:
#Pxt1({1,-1},x,t,10)
Pxt1:=proc(Steps,x,t,N) local ope,gu,i,X:
ope:=EvolutionOper(x,t,Steps,X):
gu:=1:
 for i from 1 to N do
  gu:=HafelOper(ope,gu,x,X):
 od:
gu:
end:







##################Begin GR1


with(SolveTools):
#Bdok1(DS1,P,x): checks that P in the single variable x  
#indeed fits the data set DS1
Bdok1:=proc(DaS,P,x) local i:

for i from 1 to nops(DaS) do
if subs(x=DaS[i][1],denom(P))<>0 and
normal(subs(x=DaS[i][1],P)-DaS[i][2])<>0 then

  RETURN(false):
fi:
od:
true:
end:

#GenRat1(x,dxt,dxb,a): a generic (monic top)
#rational function of degree of top dxt 
#and degree of bottom dxb
#with coeffs.parametrized by a[i] followed by the set
#of coeffs.
GenRat1:=proc(x,dxt,dxb,a) local i:

(add(a[i]*x^i,i=0..dxt))/
add(a[i+dxt+1]*x^i,i=0..dxb)
,{seq(a[i],i=0..dxt+dxb+1)}:
end:


#GR1d1d2(S1,x,d1,d2): guesses a rational function in x of degree d1/d2
#that fits the data in DS1
GR1d1d2:=proc(S1,x,d1,d2)
local eq,var,P,a,var1,i:
P:=GenRat1(x,d1,d2,a):
var:=P[2]:
P:=P[1]:

if nops(S1)<=d1+d2+3 then
# print(`Insufficient data`):
 RETURN(FAIL):
fi:

eq:={seq(numer(normal(subs(x=S1[i][1],P)-S1[i][2])),i=1..d1+d2+3)}:

var1:=Linear(eq,var):
if expand(subs(var1,numer(P)))=0 then
  RETURN(FAIL):
fi:
P:=subs(var1,P):
if Bdok1(S1,P,x) then 
   RETURN(factor(normal(P))):
 else
 RETURN(FAIL):
fi:

end:


#GR1(S1,x): guesses a rational function from the data-set S1
GR1:=proc(S1,x) local d1,d2,gu,L,i1:
if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:

for L from 0 to nops(S1)-4 do
 for d1 from 0 to L do
   d2:=L-d1:
 gu:=GR1d1d2(S1,x,d1,d2):
  if gu<>FAIL then
    RETURN(gu):
  fi:
 od:
od:

FAIL:
end:
##################End GR1




#FindrecX(f,MaxC,x,n,N,DegX): 
#guesses a recurrence operator of complexity <=MaxC,
#in n, and DegX in x
#annihilating the sequence f of polynomials in x .
#For example, try: FindrecX([seq(x^i,i=1..20)],2,x,n,N,2);
FindrecX:=proc(f,MaxC,x,n,N,DegX)
local gu,ope0,S, Deg,Ord,i,S0,j,ope,vu:

gu:=subs(x=3/4,f):

ope0:=Findrec(gu,n,N,MaxC):

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

S:={[3/4,ope0]}:

Deg:=degree(numer(normal(ope0)),n):
Ord:=degree(ope0,N):

for i from 1 to DegX+6 do
ope0:=findrec(subs(x=i+1/2,f),Deg,Ord,n,N):
 if ope0<>FAIL then
   S:=S union {[i+1/2,ope0]}:
 fi:
od:

ope:=0:

for i from 0 to Ord do
S0:=[seq([S[j][1],coeff(S[j][2],N,i)],j=1..nops(S))]:
vu:=GR1(S0,x);
 if vu=FAIL then
    RETURN(FAIL):
  else
    ope:=ope+vu*N^i:
 fi:
od:
ope:
end:




#ConjRecX(Steps,x,n,N,MaxC): Given a set of steps, Steps 
#(which are integers)
#and a positive integer x0, conjectures a homog. linear recurrence
#print(`operator, with polynomial coefficients, ope(n,N)
#of maximum complexity MaxC 
#(where N is the shift operator in n), for the sequence`):
#The  weight-enumerator according to x^(destination)
# for the walks with n steps, staying at the half-line x>=0
#using the steps of the set Steps
#The output is a pair [ope, IniConditions]
#For example, try: ConjRecX({1,-1},x,n,N,4):`):

ConjRecX:=proc(Steps,x,n,N,MaxC) local gu,N0,i,t,ope:
option remember:
N0:=trunc((4+(MaxC/2))*(2+MaxC/2)+14):

gu:=Pxt(Steps,x,t,N0):
gu:=[seq(coeff(gu,t,i),i=1..N0)]:
ope:=FindrecX(gu,MaxC,x,n,N,MaxC):
if ope=FAIL then
 RETURN(FAIL):
fi:

ope:=numer(normal(ope)):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):

[ope, [op(1..degree(ope,N),gu)]]:

end:

#Aij(t,Dt,i,j): Expanding (tDt)^i(t^(-j)) as a differential operator
#in the standard form. For example, try: Aij(t,Dt,1,1);
Aij:=proc(t,Dt,i,j) local gu:
option remember:

if i=0 then 
 RETURN(t^(-j)):
fi:

gu:=Aij(t,Dt,i-1,j):

expand(t*(diff(gu,t)+Dt*gu)):
end:


#RecToDiff(ope,n,N,t,Dt,Start1): turns a linear recurrence operator
#with polynomial coefficients ope(n,N) 
#with initial conditions Start1, into the
#differential equation [ope(tDt,1/t),rhs] . For example try:
#RecToDiff(N-n-1,n,N,t,Dt,[1,2,6,24]);
RecToDiff:=proc(ope1,n,N,t,Dt,Start1) local gu,i,j,mu,ope:
ope:=expand(numer(normal(ope1))):

gu:=0:
for i from 0 to degree(ope,n) do
for j from ldegree(ope,N) to degree(ope,N) do
 gu:=gu+coeff(coeff(ope,n,i),N,j)*Aij(t,Dt,i,j):
od:
od:

gu:=expand(numer(normal(gu))):
mu:=add(Start1[i+1]*t^i,i=0..nops(Start1)-1):
mu:=expand(
coeff(gu,Dt,0)*mu+add(coeff(gu,Dt,i)*diff(mu,t$i),i=1..degree(gu,Dt))):
mu:=add(coeff(mu,t,i)*t^i,i=0..nops(Start1)-1):
[collect(gu,Dt),mu]:
end:



#HafelDO(ope,t,Dt,F): applies the differential operator ope(t,Dt)
#to the function F(t). For example, try:
#HafelDO(Dt,t,Dt,t);
HafelDO:=proc(ope1,t,Dt,F) local i,ope,gu:
ope:=expand(ope1):
gu:=coeff(ope,Dt,0)*F:
gu:=expand(gu+add(coeff(ope,Dt,i)*diff(F,t$i), i=1..degree(ope,Dt))):
end:

#ConjDE(Steps,x0,MaxC,t,Dt):
#Given a set of steps (which are integers)
#and a positive integer, conjectures a homog. linear differential
#equation satisfied by the generating function  (w.r.t. t) of
#the sequence:
#The number of walks with n steps, staying at the half-line x>=0
#starting at the origin and winding-up and x0.
#Dt is the differentiation operator w.r.t. t.
#The output is a pair [ope, rhs], where ope
#is a linear differential operator with polynomial coefficients
#and rhs is a polynomial in t, and it means that if f(t)
#is the above-mentioned generating function, then
#ope(t,Dt)f(t)=rhs
#For example, try: ConjDE{1,-1},2,6,t,Dt):
ConjDE:=proc(Steps,x0,MaxC,t,Dt) local gu,N0,ope,i,i0,gu1,mu,n,N:
option remember:

N0:=trunc((4+(MaxC/2))*(2+MaxC/2)+10+x0):
gu1:=Seq(Steps,N0,x0):
gu:=[op(2..nops(gu1),gu1)]:


for i from 1 while gu[i]=0 do od:

i0:=i:

gu:=[op(i0+1..nops(gu),gu)]:
ope:=Findrec(gu,n,N,MaxC):

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

ope:=Yafe(subs(n=n-i0,ope),N)[2]:

gu:=RecToDiff(ope,n,N,t,Dt,gu1):

mu:=Seq(Steps,4*N0,x0):
mu:=add(mu[i+1]*t^i,i=0..4*N0):
mu:=expand(HafelDO(gu[1],t,Dt,mu)-gu[2]):

if ldegree(mu,t)<N0 then
 print(gu, `failed `):
 RETURN(FAIL):
fi:
gu:
end:





#DiA(A,t,Dt,i): The differential operator Dt^i(A(t))
#For example, try: DiA(t^4,t,Dt,3);
DiA:=proc(A,t,Dt,i) local gu:
if i=0 then
 RETURN(A):
fi:
gu:=DiA(A,t,Dt,i-1):
gu:=expand(diff(gu,t)+Dt*gu):
collect(gu,Dt):
end:



#CODVm(DE1,A,t,Dt): Given a differential equation DE1=[DO,rhs]
#phrased in terms of t and Dt, and a rational function, A, of t
#finds the differential equation satisfied by Z(t):=Y(t)/A(t)
#if Y(t) is a solution of DE1 (i.e. DO (Y(t))=rhs(t) 
#For example, try: CODVm([Dt,0],t^2,t,Dt);
CODVm:=proc(DE1,A,t,Dt)  local gu,i,ope,rhs:
ope:=DE1[1]: rhs:=DE1[2]:

gu:=add(coeff(ope,Dt,i)*DiA(A,t,Dt,i),i=0..degree(ope,Dt)):

gu:=normal(gu):
ope:=expand(numer(gu)):
ope:=collect(ope,Dt):
rhs:=rhs*denom(gu):
[ope,rhs]:

end:


#CODVa(DE1,A,t,Dt): Given a differential equation DE1=[DO,rhs]
#phrased in terms of t and Dt, and a rational function, A, of t
#finds the differential equation satisfied by Z(t):=Y(t)-A(t)
#if Y(t) is a solution of DE1 (i.e. DO (Y(t))=rhs(t) 
#For example, try: CODVa([Dt,0],t^2,t,Dt);
CODVa:=proc(DE1,A,t,Dt)  local top1,bot1,ope,rhs:
ope:=DE1[1]: rhs:=DE1[2]:
rhs:=normal(rhs-HafelDO(ope,t,Dt,A)):
top1:=numer(rhs):
bot1:=denom(rhs):
ope:=expand(bot1*ope):
[collect(ope,Dt),top1]:

end:


 
#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=expand([op(gu),normal(-add(gu[nops(gu)+1-j1]* subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L))]):
od:

gu:

end:





#SeqX(Steps,x,K): The first K terms in the weight-enumerating
#sequence for walks in the positive half-line using the steps in
#in Steps, For example, try: 
#SeqX({1,-1},x,10);
SeqX:=proc(Steps,x,K) local gu,i,t:
gu:=Pxt(Steps,x,t,K):
[seq(coeff(gu,t,i),i=1..K)]:
end:


#CheckConjDE(Steps,x0,MaxC,t,Dt,K): checks
#CheckConjDE(Steps,x0,MaxC,t,Dt) up to K-th term
#For example, try: 
#CheckConjDE({1,-1},2,6,t,Dt,20);
CheckConjDE:=proc(Steps,x0,MaxC,t,Dt,K) local ope,gu,x,i:
ope:=ConjDE(Steps,x0,MaxC,t,Dt):

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



gu:=coeff(expand(Pxt(Steps,x,t,K)),x,x0):

gu:=expand(HafelDO(ope[1],t,Dt,gu) -ope[2]):
evalb(add(coeff(gu,t,i)*t^i,i=0..K)=0):

end:



#ConjDEX1(Steps,MaxC,t,Dt,x): Conjecturs a
#diferential equation satisfies by F(x,t, Steps)
#using CODV etc. this program only allows one negative
#step: -1
#For example, try:
#ConjDEX1({1,-1},4,t,Dt,x);
ConjDEX1:=proc(Steps,MaxC,t,Dt,x) local p,s,gu:

if min(op(Steps))<-1 then
 ERROR(`This program can only take the negative step -1 `):
fi:

p:=ConjDE(Steps,0,MaxC,t,Dt):

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


p:=CODVm(p,-1/t,t,Dt):
p:=CODVa(p,-x,t,Dt):
p:=CODVm(p,expand(x-x*t*add(x^s,s in Steps)),t,Dt):

gu:=Pxt(Steps,x,t,100):


gu:=expand(HafelDO(p[1],t,Dt,gu)-p[2]) :

if ldegree(gu,t)<50 then
 RETURN(ldegree(gu,t), FAIL):
fi:
p:

end:



#DiffToRec1(ope,t,Dt,n,N): Given a linear differential operator
#ope(t,Dt) and symbols n,N, outputs the linear recurrence
#operator with polynomial coefficients satisfied by
#the sequence of coeff. of any f.p.s. annihilating ope(t,Dt)
#where N is the shift in n. For example, try:
#DiffToRec1(Dt+t, t,Dt,n,N);
DiffToRec1:=proc(ope,t,Dt,n,N) local gu,i,j,ope1,d,j1:
ope1:=expand(ope):

gu:=0:

for i from ldegree(ope1,t) to degree(ope1,t) do
 for j from 0 to degree(ope1,Dt) do
   gu:=gu+coeff(coeff(ope1,t,i),Dt,j)*mul(n-i+j1,j1=1..j)*N^(j-i):
 od:
od:


d:=ldegree(gu,N):
gu:=N^(-d)*subs(n=n-d,gu):



Yafe(gu,N)[2]:

end:

 

#FunctOper(Steps,N,x,X): Given a set of steps Steps, and
#symbols N,x,X, outputs the function operator
#P(N,x,X) annihilating a(n,x): the polynomial-generating function
#(according to final destination) for n-step one-dimensional walks
#using the set of Steps (a set of intgers)
#Here N is shift  in the n direction and X[j]f(x) is
#coeff. of x^j . For example, try:
#FunctOper({-1,1},N,x,X);
FunctOper:=proc(Steps,N,x,X) local gu,j,s:

gu:=1:

for s in Steps do
 if s>=0 then
   gu:=gu-x^s/N:
 else
  gu:=gu-x^s/N*(1-add(X[j]*x^j,j=0..-s-1)):
 fi:
od:
numer(gu):
end:


#Mult1(ope1,ope2,n,N,x,X): ope1(n,N,x,X)ope2(n,N,x)
#where X is the opeartor X f(x):=f(0)
#for example, try:
#Mult1(x*N-1-x^2+X,N^2-x*N+1,n,N,x,X);
Mult1:=proc(ope1,ope2,n,N,x,X) local ope1a,ope1b,i,gu:
ope1a:=coeff(ope1,X,0):
ope1b:=coeff(ope1,X,1):
 

gu:=0:

for i from 0 to degree(ope1a,N) do
 gu:=expand(gu+coeff(ope1a,N,i)*subs(n=n+i,ope2)*N^i):
od:

for i from 0 to degree(ope1b,N) do
 gu:=gu+coeff(ope1b,N,i)*subs({x=0,n=n+i},ope2)*N^i*X:
od:
expand(gu):

end:

#Comm1(ope1,ope2,n,N,x,X): The commutator of ope1 and ope2
Comm1:=proc(ope1,ope2,n,N,x,X):
expand(Mult1(ope1,ope2,n,N,x,X)-Mult1(ope2,ope1,n,N,x,X)):
end:



#ApplyOper1(ope,n,N,x,X,L): Given an operator ope(n,N,x,X) 
#where n, x are variables, N is the forward shift operator
#in n and X is evaluation at x=0, and a list of
#polynomials [a(1,x),a(2,x), ...] applies ope to the
#list. For example do: 
#ApplyOper1(N+x,n,N,x,X,[x,x^2,x^3]);
ApplyOper1:=proc(ope,n,N,x,X,L) local i,ord1,gu,gu1,ope0,ope1,j:
ord1:=degree(ope,N):
gu:=[]:
ope0:=coeff(ope,X,0):
ope1:=coeff(ope,X,1):
for i from 1 to nops(L)-ord1 do
gu1:=0:
  for j from 0 to degree(ope0,N) do
   gu1:=expand(gu1+subs(n=i,coeff(ope0,N,j))*L[i+j]):
  od:
  for j from 0 to degree(ope1,N) do
   gu1:=expand(gu1+subs(n=i,coeff(ope1,N,j))*subs(x=0,L[i+j])):
  od:
 gu:=[op(gu),gu1]:
od:
gu:
end:

  



#SeqAnx(Steps,x,N):
#The sequence of length N, whose i^th entry is the polynomial
#that is the weight-enumerator of all paths starting at the
#origin, staying in the non-negative half-line, and using
#the set of steps, weighted by x^(final destination)
#For example, try: SeqAnx({1,-1},x,10);
SeqAnx:=proc(Steps,x,N) local gu,i,t:
gu:=Pxt(Steps,x,t,N+2):
[seq(coeff(gu,t,i),i=1..N)]:
end:




#FindVertRec1(Steps,m,n,N, MaxC, degm): Given a set of steps Steps 
#(a set of integers) and
#discrete variables m and n, the symbol for the shift operator in n,
#N, and postitive integers MaxC and degm, 
#conjectures a recurrence operator ope(m,n,N),
#of complexity MaxC in (n,N) and degree <=degm in m
#annilihating F(m,n):= the numbr of ways
#of walking n steps on the positive discrete half-line, starting at
#0, and ending in m. For example, try:
#FindVertRec1({1,-1},m,n,N,4,4);
FindVertRec1:=proc(Steps,m,n,N, MaxC, degm) local S0,m0,
ope0,Ord1,s,S0a,gu,i,mu:

S0:={}:
for m0 from 0 to degm+4 do
 ope0:=ConjRec(Steps,m0,n,N,MaxC):
  if ope0=FAIL then
    RETURN(FAIL):
  fi:

S0:=S0 union {[m0,ope0]}:
od:



Ord1:=max(seq(degree(s[2],N),s in S0)):

gu:=0:

for i from 0 to Ord1 do

 S0a:={seq([s[1],coeff(s[2],N,i)],s in S0)}:

  mu:=GR1(S0a,m):

  if mu=FAIL then
      RETURN(FAIL):
   fi:
 gu:=gu+factor(mu)*N^i:

od:

gu:

end:




#FindVertRec(Steps,m,n,N, degn,degm,degN): Given a set of steps Steps 
#(a set of integers) and
#discrete variables m and n, the symbol for the shift operator in n, N,
#and postitive integers degn , degm, degN
#conjectures a recurrence operator ope(m,n,N),
# of degree <=degm in m, degree <=degn in n, degree <=degN in N
#annilihating F(m,n):= the number of ways
#of walking n steps on the positive discrete half-line, starting at
#0, and ending in m. For example, try:
#FindVertRec({1,-1},m,n,N,2,2,2);
FindVertRec:=proc(Steps,m,n,N, degn, degm,degN) local eq,eq1,var,a,ope,
g,i,j,k,co,m0,n0,lu:
ope:=0:
var:={}:

for i from 0 to degn do
 for j from 0 to degm do
  for k from 0 to degN do
   var:=var union {a[i,j,k]}:
   ope:=ope+a[i,j,k]*n^i*m^j*N^k :
  od:
 od:
od:


co:=0:

eq:={}:
for g from 3 do 
for  n0 from  trunc(g/2)+1 to g  do
m0:=g-n0:
 co:=co+1:
 eq1:=add(subs({m=m0,n=n0},coeff(ope,N,k))*NuWalks(Steps,n0+k,m0),k=0..degN):
  if eq1<>0 then
  eq:=eq union {eq1}:
 fi:
if co>nops(var)+10 then
 var:=solve(eq,var):
 ope:=subs(var,ope):
   if ope=0 then
    RETURN(FAIL):
    else
      lu:={}:
      for i from 1 to nops(var) do
         if op(1,op(i,var))=op(2,op(i,var)) then
           lu:=lu union {op(1,op(i,var))}:
         fi:
      od:
     if nops(lu)>1 then
      print(`too much slack`):
       RETURN(ope):
      else
     RETURN(Yafe(ope,N)[2]):
     fi:
   fi:
fi:
od:
od:
FAIL:
end:




#FindHorizRec(Steps,m,n,M, degn,degm,degM): Given a set of steps Steps 
#(a set of integers) and
#discrete variables m and n, the symbol for the shift operator in m, M,
#and postitive integers degn , degm, degN
#conjectures a recurrence operator ope(m,n,M),
# of degree <=degm in m, degree <=degn in n, degree <=degM in M
#annilihating F(m,n):= the number of ways
#of walking n steps on the positive discrete half-line, starting at
#0, and ending in m. For example, try:
#FindHorizRec({1,-1},m,n,M,2,2,2);
FindHorizRec:=proc(Steps,m,n,M, degn, degm,degM) local eq,eq1,var,a,ope,
g,i,j,k,co,m0,n0,lu:
ope:=0:
var:={}:

for i from 0 to degn do
 for j from 0 to degm do
  for k from 0 to degM do
   var:=var union {a[i,j,k]}:
   ope:=ope+a[i,j,k]*n^i*m^j*M^k :
  od:
 od:
od:


co:=0:

eq:={}:
for g from 3 do 
for  n0 from  trunc(g/2)+1 to g  do
m0:=g-n0:
 co:=co+1:
 eq1:=add(subs({m=m0,n=n0},coeff(ope,M,k))*NuWalks(Steps,n0,m0+k),k=0..degM):
  if eq1<>0 then
  eq:=eq union {eq1}:
 fi:
if co>nops(var)+10 then
 var:=solve(eq,var):
 ope:=subs(var,ope):
   if ope=0 then
    RETURN(FAIL):
    else
      lu:={}:
      for i from 1 to nops(var) do
         if op(1,op(i,var))=op(2,op(i,var)) then
           lu:=lu union {op(1,op(i,var))}:
         fi:
      od:
     if nops(lu)>1 then
      print(`too much slack`):
       RETURN(ope):
      else
     RETURN(Yafe(ope,M)[2]):
     fi:
   fi:
fi:
od:
od:
FAIL:
end:


#FindCCRec(Steps,M,N,degM,degN): Given a set of steps Steps 
#(a set of integers) and
#discrete variables m and n, the symbol for the shift operator in m, M,
#and in n, N,
#and postitive integers degM , degN
#conjectures a recurrence operator ope(M,N),
# of  degree <=degM in M and degree<=degN in N
#annilihating F(m,n):= the number of ways
#of walking n steps on the positive discrete half-line, starting at
#0, and ending in m. For example, try:
#FindCCRec({1,-1},M,N,2,2);
FindCCRec:=proc(Steps,M,N,degM,degN) local eq,eq1,var,a,ope,
g,i,j,k,co,m0,n0,lu,k1:
ope:=0:
var:={}:

  for i from 0 to degM do
  for j from 0 to degN do
   var:=var union {a[i,j]}:
   ope:=ope+a[i,j]*M^i*N^j :
 od:
od:


co:=0:

eq:={}:
for g from 3 do 
for  n0 from  trunc(g/2)+1 to g  do
m0:=g-n0:
 co:=co+1:
 eq1:=
add(
add(
 coeff(coeff(ope,M,k),N,k1)*NuWalks(Steps,n0+k1,m0+k),
   k1=0..degN)
   ,k=0..degM) :
  if eq1<>0 then
  eq:=eq union {eq1}:
 fi:
if co>nops(var)+10 then
 var:=solve(eq,var):
 ope:=subs(var,ope):
   if ope=0 then
    RETURN(FAIL):
    else
      lu:={}:
      for i from 1 to nops(var) do
         if op(1,op(i,var))=op(2,op(i,var)) then
           lu:=lu union {op(1,op(i,var))}:
         fi:
      od:
     if nops(lu)>1 then
      print(`too much slack`):
       RETURN(ope):
      else
     RETURN(ope):
     fi:
   fi:
fi:
od:
od:
FAIL:
end:








#FindGoodmRec(Steps,m,n,N, M, degn,degm,degN,degM): 
#Given a set of steps Steps 
#(a set of integers) and
#discrete variables m and n, the symbol for the shift operator in n, N,
#and postitive integers degm, degn , degm, degN
#conjectures a recurrence operator ope(m,n,N),
#of the form P(n,N)+mQ(n,m,M,N) 
# of degree <=degm in m, degree <=degn in n, degree <=degN in N
#annilihating F(m,n):= the number of ways
#of walking n steps on the positive discrete half-line, starting at
#0, and ending in m. For example, try:
#FindGoodmRec({1,-1},m,n,N,M,2,2,2,2);
FindGoodmRec:=proc(Steps,m,n,N, M, degn, degm,degN,degM) 
local eq,eq1,var,a,ope1,ope2,ope, 
g,i,j,k,co,m0,n0,lu,k1,b,fu:
ope1:=0:
var:={}:

for i from 0 to degn do
 for j from 1 to degm do
  for k from 0 to degN do
  for k1 from 0 to degM do
   var:=var union {a[i,j,k,k1]}:
   ope1:=ope1+a[i,j,k,k1]*n^i*m^j*N^k*M^k1 :
  od:
 od:
od:
od:

ope2:=0:
for i from 0 to degn do
 for j from 0 to degN do
   var:=var union {b[i,j]}:
   ope2:=ope2+b[i,j]*n^i*N^j :

od:
od:
ope:=expand(ope1+ope2):


co:=0:

eq:={}:
for g from 15 do 
for  n0 from  trunc(g/2) to g-2  do
m0:=g-n0:


 eq1:=
add(
add(
 subs({n=n0,m=m0},coeff(coeff(ope,M,k),N,k1))*NuWalks(Steps,n0+k1,m0+k),
   k=0..degM)
   ,k1=0..degN) :


  if eq1<>0 then
  co:=co+1:
  eq:=eq union {eq1}:
 fi:
if co>nops(var)+15 then

 var:=solve(eq,var):
 ope1:=subs(var,ope1):
 ope2:=subs(var,ope2):
   if ope2=0 then
    RETURN(FAIL):
    else
      lu:={}:
      for i from 1 to nops(var) do
         if op(1,op(i,var))=op(2,op(i,var)) then
           lu:=lu union {op(1,op(i,var))}:
         fi:
      od:

 fu:={}:

 for  i in lu do

  if coeff(ope2,i,1)<>0 then
   fu:= fu union {[ coeff(ope2,i,1), coeff(ope1,i,1)]}:
  fi:

 od:

RETURN(fu):
fi:
fi:
od:
od:
FAIL:
end:




#FindGoodmRec1(Steps,m,n,N, M, degn,degm,degN,degM,MaxC): 
#Given a set of steps Steps 
#(a set of integers) and
#discrete variables m and n, the symbol for the shift operator in n, N,
#and postitive integers degm, degn , degm, degN
#conjectures a recurrence operator ope(m,n,N),
#of the form P(n,N)+mQ(n,m,M,N) 
# of degree <=degm in m, degree <=degn in n, degree <=degN in N
#annilihating F(m,n):= the number of ways
#of walking n steps on the positive discrete half-line, starting at
#0, and ending in m. For example, try:
#FindGoodmRec({1,-1},m,n,N,M,2,2,2,2);
FindGoodmRec1:=proc(Steps,m,n,N, M, degn, degm,degN,degM,MaxC) 
local eq,eq1,var,a,ope1,ope2,ope, 
g,i,j,k,co,m0,n0,lu,k1,fu:
ope1:=0:
var:={}:

for i from 0 to degn do
 for j from 1 to degm do
  for k from 0 to degN do
  for k1 from 0 to degM do
   var:=var union {a[i,j,k,k1]}:
   ope1:=ope1+a[i,j,k,k1]*n^i*m^j*N^k*M^k1 :
  od:
 od:
od:
od:

ope2:=numer(ConjRec(Steps,0,n,N,MaxC)):
ope:=expand(ope1+ope2):


co:=0:

eq:={}:
for g from 15 do 
for  n0 from  trunc(g/2) to g-2  do
m0:=g-n0:


 eq1:=
add(
add(
 subs({n=n0,m=m0},coeff(coeff(ope,M,k),N,k1))*NuWalks(Steps,n0+k1,m0+k),
   k=0..degM)
   ,k1=0..degN) :


  if eq1<>0 then
  co:=co+1:
  eq:=eq union {eq1}:
 fi:
if co>nops(var)+15 then

 var:=solve(eq,var):


   if  var=NULL then
    RETURN(FAIL):
    else
 ope1:=subs(var,ope1):
      lu:={}:
      for i from 1 to nops(var) do
         if op(1,op(i,var))=op(2,op(i,var)) then
           lu:=lu union {op(1,op(i,var))}:
         fi:
      od:

 fu:={}:

 for  i in lu do

   fu:= fu union {coeff(ope1,i,1)}:

 od:

RETURN(fu):
fi:
fi:
od:
od:
FAIL:
end: