##############################################################################################
##This is Proj2.txt, a Maple package to generate and investigate integer sequences counting  #
#Lattice Walks to the diagonal                                                               #
#in the 2-Dimensional Manhattan Lattice for an arbitrary set of atomic steps                 #
#It also investigates walks that stay in x>=y `):                                            #
####It is the Maple package created by Team X in Dr. Z.'s Combinatorics Clss at              #
#Rutgers University, Fall 2020.                                                              #
# Save this file as `Proj2.txt`, to use it                                                   #
#Type, in a Maple session                                                                    #
#read `ProjX.txt`):                                                                          #
#and then to get a list of the functions type                                                #
#Help():                                                                                     #
#For Help with any of the functions, type                                                    # 
#Help(FunctionName)                                                                          #
#Team Leader:                                                                                #
#Team members:                                                                               #
##############################################################################################


print(`This is Proj2.txt, a Maple package that is part of Project 2 in Dr. Z.'s Combinatorics Class at Rutgers University`):
print(``):
print(`Team Leader: tbd `):
print(``):
print(`Other Team members: tbd `):
print(`For a list of all the functions type: Help(); `):
print(`For Help with any of the functions, type  Help(FunctionName):`):                                                                        

Help:=proc()
if nargs=0 then
print(`The available procedures are`):
print(` CriExp, FindRec, GrowthConst,  NuGW, NuW, SeqFromRec, SeqGW, SeqW `): 

elif nargs=1 and args[1]=CriExp then
print(`CriExp(ope,n,N): The critical exponent of the operator ope in n and N Try`):
print(`CriExp((n+1)/(n+2)-3*(2*n+3)/(n+2)*N+N^2,n,N):`):


elif nargs=1 and args[1]=Findrec then
print(`Findrec(L,n,N,MaxC): Given a list L 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 in n  `):
print(`of maximum COMPLEXCITY, DEGREE+ORDER<=MaxC`):
print(`If none exists, it returns FAIL. . Try:`):
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]=GrothConst then
print(`GrowthConst(ope,n,N): The growth constant of sequence of integers satistying the recurrence given by the`):
print(`shift operator, ope. Try: `):
print(`GrowthConst(N^2-N-1,n,N);`):



elif nargs=1 and args[1]=NuGW then
print(`NuGW(Pt,A): Input a Lattice Pt pt of the form [m,n],  and a set A of "atomic steps" (pairs with non-negative integers)`):
print(`(except that  [0,0] is not allowed) outputs the number of Good walks from [0,0] to Pt using these atomic steps.`):
print(`A good walk is a walk that always stays in the region x>=y. `):
print(`Try:`):
print(`NuGW([6,4],{[0,1],[1,0]});`):

elif nargs=1 and args[1]=NuW then
print(`NuW(Pt,A): Input a Lattice Pt pt of the form [m,n],  and a set A of "atomic steps" (pairs with non-negative integers)`):
print(`(except that  [0,0] is not allowed) outputs the number of walks from [0,0] to Pt using these atomic steps.`):
print(`Try:`):
print(`NuW([6,4],{[0,1],[1,0]});`):

elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): inputs a linear recurrence operator with constant coefficients`):
print(`in the discrete variable n and the shift operator N, and a list of intetegers Ini whose length is the order of ope`):
print(`outputs the list of length K of the sequence that satisfies the unerlying recurrence`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10);`):


elif nargs=1 and args[1]=SeqGW then
print(`SeqGW(A,K): inputs a set of atomic walks, A, in the 2D lattice , and a positive integer K, outputs the`):
print(`first K terms of the sequence "Number of GOOD (staying in x>=y) walks from [0,0] to [i,i]" using the steps in A`):
print(`Try:`):
print(`SeqGW({[1,0],[0,1]},20);`):

elif nargs=1 and args[1]=SeqW then
print(`SeqW(A,K): inputs a set of atomic walks, A, in the 2D lattice , and a positive integer K, outputs the`):
print(`first K terms of the sequence "Number of walks from [0,0] to [i,i]" using the steps in A`):
print(`Try:`):
print(`SeqW({[1,0],[0,1]},20);`):

else
 print(`There is no Help for`):

print(args):
fi:

end:

 


##Start CODE taken from Dr. Z.'s Maple package Findrec.txt
ezra:=proc()
if args=NULL then

print(` FindRec.txt: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of ONE Variable`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  `):
print(`  findrec, Findrec, FindrecF`):
print():

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]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. .g. try FindrecF(i->i!,n,N);`):


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`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10);`):


fi:
 
end:


###Findrec 
#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);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+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)+2 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
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

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

#GrowthConst(ope,n,N): The growth constant of the operator ope in n and the shift operator N. Try
#GrowthConst(N^2-N-1,n,N);
GrowthConst:=proc(ope,n,N) local ope1,champ,i,gu,rec:
ope1:=lcoeff(numer(ope),n):
gu:=[solve(ope1,N)]:
champ:=gu[1]:
rec:=evalf(abs(gu[1])):

for i from 2 to nops(gu) do
 if evalf(abs(gu[i]))>rec then
  champ:=gu[i]:
 rec:=evalf(abs(gu[i])):
 fi:
od:
champ:
end:


###END CODE taken from Dr. Z.'s Maple package Findrec.txt

#NuW(Pt,A): Input a Lattice Pt pt of the form [m,n],  and a set A of "atomic steps" (pairs with non-negative integers)
#(except that  [0,0] is not allowed) outputs the number of walks from [0,0] to Pt using these atomic steps.
#Try:
#NuW([6,4],{[0,1],[1,0]});
NuW:=proc(Pt,A) local a:
option remember:

if Pt[1]<0 or Pt[2]<0 then
 RETURN(0):
elif Pt=[0,0] then
 RETURN(1):
else
add(NuW(Pt-a,A),a in A):
fi:
end:


#NuGW(Pt,A): Input a Lattice Pt pt of the form [m,n],  and a set A of "atomic steps" (pairs with non-negative integers)
#(except that  [0,0] is not allowed) outputs the number of walks from [0,0] to Pt that always stay in x>=y
#using these atomic steps.
#Try:
#NuGW([6,4],{[0,1],[1,0]});
NuGW:=proc(Pt,A) local a:
option remember:

if (Pt[1]<0 or Pt[2]<0 or Pt[1]<Pt[2]) then
 RETURN(0):
elif Pt=[0,0] then
 RETURN(1):
else
add(NuGW(Pt-a,A),a in A):
fi:
end:



#SeqW(A,K): inputs a set of atomic walks, A, in the 2D lattice , and a positive integer K, outputs the
#first K terms of the sequence "Number of walks from [0,0] to [i,i]" using the steps in A
#Try:
#SeqW({[1,0],[0,1]},20);
SeqW:=proc(A,K) local i:
[seq(NuW([i,i],A),i=1..K)]:
end:


#SeqGW(A,K): inputs a set of atomic walks, A, in the 2D lattice , and a positive integer K, outputs the
#first K terms of the sequence "Number of Good walks from [0,0] to [i,i]" using the steps in A
#Try:
#SeqGW({[1,0],[0,1]},20);
SeqGW:=proc(A,K) local i:
[seq(NuGW([i,i],A),i=1..K)]:
end:


#CriExp(ope,n,N): The critical exponent of the operator ope in n and N Try
#CriExp((n+1)/(n+2)-3*(2*n+3)/(n+2)*N+N^2,n,N):
CriExp:=proc(ope,n,N) local A,ope1,x,c,t,t0,i,eq:

ope1:=numer(ope):

A:=numer(normal(add(subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^t*c^i,i=0..degree(ope1,N)))):

A:=taylor(A,x=0,4):
eq:={coeff(A,x,0),coeff(A,x,1)}:
t0:=solve(eq,{t,c}):
if t0=NULL then
RETURN(FAIL):
fi:
subs(t0,t):
end: