Help:=proc():
	print(`The Main procedures are:`):
	print(`walksn(S,n)`, `numwalksn(S,n)`,`numwalksK(S,K)`, `findrecwalks(S,K,n,N,maxC), findrecwalksE(S,n,N)`): 
	print(`walksnpos(S,n)`, `numwalksnpos(S,n)`, `numwalksKpos(S,K)`, 	`findrecposwalks(S,K,n,N,maxC)`):
        print(`AZd(F,x,n,N) `):
	#print(`walksnm(S,n,m)`, `walksn(S,n)`, `walksK(S,K)`):
	#print(`walksnmpos(S,n,m)`, `walksnpos(S,n)`, `walksKpos(S,K)`):
	#print(`numwalksnm(S,n,m)`, `numwalksn(S,n)`, `numwalksK(S,K)`):
	#print(`numwalksnmpos(S,n,m)`, `numwalksnpos(S,n)`, `numwalksKpos(S,K)`):
end:


#walksnm(S,n,m): the set of walks of length n with steps in S={[1,b1],[1,b2],...} from [0,0] to [n,m]
walksnm:=proc(S,n,m) local W, S2,s,W1,w: option remember:
	W:={}:
	S2:={seq(s[2],s in S)}:

	if n=0 then
		RETURN({[]}):

	elif n=1 then 	
			if member(m,S2)=true then
				W:={[m]}:
			else
				RETURN({}):
			fi:
	else
	
	for s in S2 do
		W1:=op({walksnm(S,n-1,m-s)} minus {{}}): 
		if W1<>{} then
			W:=W union {seq([op(w),s],w in W1)}:
		fi:
	od:
	fi:
	W:
end:


#walksn(S,n): the set of walks of length n with steps in S={[1,b1],[1,b2],...} from [0,0] to [n,0]
walksn:=proc(S,n) local W,w,i:
	W:=walksnm(S,n,0):
	{seq([seq([1,w[i]],i=1..nops(w))],w in W)}:
end:

#walksK(S,K): the sequence enumerating walks of length n with steps in S={[1,b1],[1,b2],...} from [0,0] to [n,0] from n=1 to n=k
walksK:=proc(S,K) local W,n:
	seq(nops(walksn(S,n)),n=1..K):
end:


#walksnmpos(S,n,m): the set of walks of length n with steps in S={[1,b1],[1,b2],...} from [0,0] to [n,m]that have non-negative height
walksnmpos:=proc(S,n,m) local W, S2,s,W1,w: option remember:
	W:={}:
	S2:={seq(s[2],s in S)}:

	if n=0 then
		RETURN({[]}):

	elif n=1 then 	
			if (member(m,S2)=true and m>=0) then
				W:={[m]}: 
			else
				RETURN({}):
			fi:
	else
	
	for s in S2 do 
		if m-s>=0 then
		W1:=op({walksnmpos(S,n-1,m-s)} minus {{}}): 
		if W1<>{} then
			W:=W union {seq([op(w),s],w in W1)}:
		fi: fi:
	od:
	fi:
	W:
end:



#walksnpos(S,n): the set of walks of length n with steps in S={[1,b1],[1,b2],...} from [0,0] to [n,0]that never have negative height
walksnpos:=proc(S,n) local W,i,w:
	W:=walksnmpos(S,n,0):
	{seq([seq([1,w[i]],i=1..nops(w))],w in W)}:
end:

#walksKpos(S,K): the sequence enumerating walks of length n with steps in S={[1,b1],[1,b2],...} from [0,0] to [n,0] from n=1 to n=k that never have negative height
walksKpos:=proc(S,K) local W,n:
	seq(nops(walksnpos(S,n)),n=1..K):
end:





numwalksnm:=proc(S,n,m) local W, S2,s,W1: option remember:
	W:=0:
	S2:={seq(s[2],s in S)}:

	if n=0 then
		RETURN(1):

	elif n=1 then 	
			if member(m,S2)=true then
				RETURN(1):
			else
				RETURN(0):
			fi:
	else
	
	for s in S2 do
		W1:=numwalksnm(S,n-1,m-s): 
		if W1<>{} then
			W:=W +W1:
		fi:
	od:
	fi:
	W:
	
end:


#numwalksn(S,n): the number of walks of length n with steps in S={[1,b1],[1,b2],...} from [0,0] to [n,0]
numwalksn:=proc(S,n) local W:
	numwalksnm(S,n,0):
end:


#numwalksK(S,K): the sequence enumerating walks of length n with steps in S={[1,b1],[1,b2],...} from [0,0] to [n,0] from n=1 to n=k
numwalksK:=proc(S,K) local W,n:
	[seq(numwalksnm(S,n,0),n=1..K)]:
end:


numwalksnmpos:=proc(S,n,m) local W, S2,s,W1: option remember:
	W:=0:
	S2:={seq(s[2],s in S)}:

	if n=0 then
		RETURN(1):

	elif n=1 then 	
			if member(m,S2)=true and m>=0 then
				RETURN(1):
			else
				RETURN(0):
			fi:
	else
	
	for s in S2 do
		if m-s>=0 then
			W1:=numwalksnmpos(S,n-1,m-s): 
			if W1<>{} then
				W:=W +W1:
			fi:
		fi:
	od:
	fi:
	W:
	
end:


#numwalksnpos(S,n): the number of walks of length n with steps in S={[1,b1],[1,b2],...} from [0,0] to [n,0]
numwalksnpos:=proc(S,n) local W,w,i:
	numwalksnmpos(S,n,0):
end:


#numwalksKpos(S,K): the sequence enumerating walks of length n with steps in S={[1,b1],[1,b2],...} from [0,0] to [n,0] from n=1 to n=k
numwalksKpos:=proc(S,K) local W,n:
	[seq(numwalksnpos(S,n),n=1..K)]:
end:

findrecposwalks:=proc(S,K,n,N,maxC):
	Findrec(numwalksKpos(S,K),n,N,maxC):
end:
	

findrecwalks:=proc(S,K,n,N,maxC):
	Findrec(numwalksK(S,K),n,N,maxC):
end:



#findrecwalksE(S,n,N,): a recurrence for the number of walks from [0,0] to [n,0] w/o guessing, using the Almkvist-Zeilberger algorithm. Try:
#findrecwalksE({[1,0],[1,1],[1,-1]},n,N):
findrecwalksE:=proc(S,n,N) local t,ope,i,f,lu:
if {seq(S[i][1],i=1..nops(S))}<>{1} then
 print(`All x-coordinates of the steps of`, S , `must be 1 `):
RETURN(FAIL):
fi:

f:=(add(t^S[i][2],i=1..nops(S)))^n/t:
ope:=AZd(f,t,n,N):
if ope=0 then
 RETURN(FAIL):
else
ope:=ope[1]:
fi:

lu:=lcoeff(ope,N):
ope:=add(normal(coeff(ope,N,i)/lu)*N^i,i=0..degree(ope,N)):
ope:
end:









####################	From Findrec.txt	####################

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





 
#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:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:
 
#end Findrec


with(linalg):

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



##FROM EKHAD.txt

 
goremd:=proc(N,ORDER)
local i, gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*N^i:
od:
RETURN(gu):
end:

pashetd:=proc(p,N)
local gu,p1,mekh,i:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu,mekh):
end:

duisd:= proc(A,ORDER,y,n,N)
local gorem, conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,shad,KAMA,i1:
KAMA:=10:
gorem:=goremd(N,ORDER):
 
 
conj:=gorem:
yakhas:=0:
 
for i1 from 0 to degree(conj,N) do
 
yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve(eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):
 
 
shad:=pashetd(gorem,N):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:
 
RETURN(shad[1],S):
end:
  
AZd:= proc(A,y,n,N)
local ORDER,gu,KAMA:
KAMA:=20:
 
for ORDER from 0 to KAMA do
gu:=duisd(A,ORDER,y,n,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
0:
end:

#END from EKHAD.txt