######################################################################
##LOZITSKY.txt: Save this file as  LOZITSKY.txt                      #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read LOZITSKY.txt<Enter>                           #
##Then follow the instructions given there                           #
##                                                                   #
##Written by George Spahn and                                        #
#Doron Zeilberger, Rutgers University ,                              #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Fall 2023
 
print(`Created:  Fall  2023`):
print(` This is  LOZITSKY.txt `):
print(`It a package that was written after the publication of the article`):
print(`  Experimenting with Discrete Dynamical Systems`):
print(`by George Spahn and Doron Zeilberger`):
print(``):
print(`but it still accompanies it`):
print(``):
print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/LOZITSKY.txt  .`):


print(`---------------------------------------`):
 print(`For a list of the Story procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):



print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

#Digits:=50:
with(combinat):
with(linalg):
with(Optimization):

with(numtheory):

ezraS:=proc()

if args=NULL then
 print(` The story procedures are:  Paper `):

 print(``):

else
ezra(args):
fi:

end:


 

ezra1:=proc()
 
if args=NULL then
 print(` The supporting procedures are:  Coes, InvR, OS, YafeR`):

 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  `):
 print(` DB, DB1, PR, PR1, Traj`):




elif nops([args])=1 and op(1,[args])=Coes then
print(`Coes(P,X): Given a polynomial P in the list of variables X finds the set of coefficients. Try:`):
print(`Coes((x+y)^10,[x,y]);`):

elif nops([args])=1 and op(1,[args])=DB then
print(`DB(z,n): A data base of previously known recurrence including complex and algenraic coefficients in the form [k,F] where k is the order and`):
print(`F in the expression for z[n] in terms of z[n-1], ..., z[n-k]. Try:`):
print(`DB(z,n);`):


elif nops([args])=1 and op(1,[args])=DB1 then
print(`DB1(z,n): A data base of previously known recurrence with integer coefficients in the form [k,F] where k is the order and`):
print(`F in the expression for z[n] in terms of z[n-1], ..., z[n-k]. Try:`):
print(`DB1(z,n);`):

elif nops([args])=1 and op(1,[args])=InvR then
print(`InvR(T,k,z,n): The inverse of the recurrence of order k of the form z[n] a function of z[n-1], ..., z[n-k]. Try:`):
print(`InvR((1+z[n-1])/z[n-2],2,z,n);`):

elif nops([args])=1 and op(1,[args])=OS then
print(`OS(L,T,z,n,k): Given an expression T in z[n-k-1], z[n-k],z[n-k+1],..., z[n-1] where L has at least length k+1, finds the next term in the recurrence`):
print(`z[n]=T(z[n-k-1],...,...z[n-1]). Try:`):
print(`OS([x,y],z[n-1]/z[n-2],z,n,2); `):


elif nops([args])=1 and op(1,[args])=Paper then
print(`Paper(K,R,z,n): Inputs a positive integers K, R, outputs  a paper with periodic recurrences of the from`):
print(`z[n]=(1+a1*z[n-1]+...+ak*z[n-k])/z[n-k] with periods from 2 to R (only that exist).`):
print(`It is inspired by Evgeni Lozitsky's paper`):
print(`Fractional Linear Periodic Recurrences of orders two and three`):
print(`https://arxiv.org/abs/2310.13036. Try:`):
print(`Paper(2,12,z,n);`):


elif nops([args])=1 and op(1,[args])=PR then
print(`PR(k,r,z,n): All the peoridic recurrence of period r of the form z[n]=(a[0]+add(a[i]*z[n-i],i=1..k-1))/z[n-k],`):
print(`Try:`):
print(`PR(2,5,z,n);`):


elif nops([args])=1 and op(1,[args])=PR1 then
print(`PR1(k,r,z,n): All the peoridic recurrence of period r of the form z[n]=add(a[i]*z[n-i],i=1..k)/z[n-k],`):
print(`Try:`):
print(`PR1(2,5,z,n);`):

elif nops([args])=1 and op(1,[args])=Traj then
print(`Traj(L,T,z,n,k,K): starting with a list L finds the trajectory of applying OS(L,T,z,n,k). Try:`):
print(`Traj([x,y],z[n-1]/z[n-2],z,n,2,10);`):

elif nops([args])=1 and op(1,[args])=YafeR then
print(`YafeR(T,z,n,k): Given an expression T whose numerator has the form (1+a[1]*z[n-1]+...+a[k]*z[n-k]) retuns it nicely. Try:`):
print(`YafeR(PR1(2,8,z,n)[2],z,n,2);`):
else
print(`There is no ezra for`,args):
fi:
 
end:


#ez:=proc(): print(`DB(z,n), DB1(z,n), OS(L,T,z,n,k), Traj(L,T,z,n,k,K), Prec(k,r,z,n), Prec1(k,r,z,n), Vecs(k,N) , PrecS(k,r,z,n,m,N), Prec1S(k,r,z,n,m,N) `):end:

#DB(z,n): A data base of previously known recurrence including complex and algenraic coefficients in the form [k,F] where k is the order and
#F in the expression for z[n] in terms of z[n-1], ..., z[n-k]. Try:
#DB(z,n)
DB:=proc(z,n):
[
[1,1/z[n-1]],

[2,z[n-1]/z[n-2]],

[2,(1+z[n-1])/z[n-2]],

[3,(1+z[n-1]+z[n-2])/z[n-3]],

[3,(-1-z[n-1]+z[n-2])/z[n-3]],

[2,(I*z[n-2]+z[n-1]+(1-I)/2 )/z[n-2]],

[2,(-I*z[n-2]+z[n-1]+(1+I)/2 )/z[n-2]],

#[2,(-2*I*z[n-2]+2*z[n-1]+(2-sqrt(3)+I))/(2*z[n-2])],

[2,(-2*I*z[n-2]+2*z[n-1]+(2+sqrt(3)+I))/(2*z[n-2])],

[3,(-1/2-z[n-3]+z[n-2]-z[n-1])/z[n-3]],

[3,((1+sqrt(3)*I)*z[n-3]+2*z[n-2]+(1-sqrt(3)*I)*z[n-1]+(-1-sqrt(3)*I))/(2*z[n-3])],

[3,((1-sqrt(3)*I)*z[n-3]+2*z[n-2]+(1+sqrt(3)*I)*z[n-1]+(-1+sqrt(3)*I))/(2*z[n-3])]
]:
end:


#DB1(z,n): A data base of previously known recurrence with only integer coefficients in the form [k,F] where k is the order and
#F in the expression for z[n] in terms of z[n-1], ..., z[n-k]. Try:
#DB1(z,n);
DB1:=proc(z,n):
[
[1,1/z[n-1]],

[2,z[n-1]/z[n-2]],

[2,(1+z[n-1])/z[n-2]],

[3,(1+z[n-1]+z[n-2])/z[n-3]],

[3,(-1-z[n-1]+z[n-2])/z[n-3]],

[3,(-1/2-z[n-3]+z[n-2]-z[n-1])/z[n-3]],

[2,z[n-1]*(-1+2*z[n-1])/2/(-z[n-2]+z[n-1]+z[n-1]*z[n-2])]
]:
end:



#OS(L,T,z,n,k): Given an expression T in z[n-k-1], z[n-k],z[n-k+1],..., z[n-1] where L has at least length k+1, finds the next term in the recurrence
#z[n]=T(z[n-k-1],...,...z[n-1]). Try:
#OS([x,y],z[n-1]/z[n-2],z,n,2).
OS:=proc(L,T,z,n,k) local kha,i:
if nops(L)<k then
 RETURN(FAIL):
fi:
 if member(0,{op(L)}) then
  RETURN(FAIL):
 fi:
kha:=normal(subs({seq(z[n-i]=L[-i],i=1..k)},T)):
[op(L),kha]:

end:

#Traj(L,T,z,n,k,K): starting with a list L finds the trajectory of applying OS(L,T,z,n,k). Try:
#Traj([x,y],z[n-1]/z[n-2],z,n,1,10).
Traj:=proc(L,T,z,n,k,K) local L1,i:
L1:=L:
for i from 1 to K do
 L1:=OS(L1,T,z,n,k):
 if L1=FAIL then
   RETURN(FAIL):
 fi:
od:
L1:
end:



#Coes(P,X): Given a polynomial P in the list of variables X finds the set of coefficients. Try:
#Coes((x+y)^10,[x,y]);
Coes:=proc(P,X) local k,x,i,X1:
k:=nops(X):
x:=X[k]:

 if k=1 then
  RETURN({seq(coeff(P,x,i),i=ldegree(P,x)..degree(P,x))}) minus {0}:
 fi:
X1:=[op(1..k-1,X)]:
{seq(op(Coes(coeff(P,x,i),X1)),i=ldegree(P,x)..degree(P,x))} minus {0}:
  

end:

#PRold(k,r,z,n): all the k-th order periodic recurrences of order k of the form z[n]=(a[0]+...a[k-1]*z[n-k+1])/z[n-k] of period r. Try:
#PRold(2,5,z,n);
PRold:=proc(k,r,z,n) local f,a,i,F,F1,var,eq,x,var1:
var:={seq(a[i],i=0..k-1)}:
#f:=(a[0]+add(a[i]*z[n-i],i=1..k-1))/z[n-k]:
f:=(1+add(a[i]*z[n-i],i=1..k-1))/z[n-k]:
F:=[seq(x[i],i=1..k)]:
F1:=Traj(F,f,z,n,k,r):
eq:={seq(op(Coes(numer(F1[i]-F1[i+r]),F)),i=1..k)}:
#var1:=solve(eq,var):
eq:=Groebner[Basis](eq,plex(seq(a[i],i=0..k-1))):
var1:=solve(eq,var):
if var1=NULL then
 RETURN(FAIL):
fi:
if f=0 then
 RETURN(FAIL):
fi:

f:=subs(var1,f):

if f=0 then
 RETURN(FAIL):
fi:
f:
end:


#PR1old(k,r,z,n): all the k-th order periodic recurrences of order k of the form z[n]=(a[0]+...a[k-1]*z[n-k])/z[n-k] of period r. Try:
#PR1old(2,5,z,n);
PR1old:=proc(k,r,z,n) local f,a,i,F,F1,var,eq,x,var1,var11,gu,f1:
var:={seq(a[i],i=1..k)}:
f:=(1+add(a[i]*z[n-i],i=1..k))/z[n-k]:
F:=[seq(x[i],i=1..k)]:
F1:=Traj(F,f,z,n,k,r):
eq:={seq(op(Coes(numer(F1[i]-F1[i+r]),F)),i=1..k)}:

eq:=Groebner[Basis](eq,plex(seq(a[i],i=1..k))):

var1:={solve({op(eq)},var)}:




if var1=NULL then
 RETURN(FAIL):
fi:

gu:={}:

for var11 in var1 do
 f1:=subs(var11,f):
 gu:=gu union {f1}:
od:

gu:
end:

#InvR(T,k,z,n): The inverse of the recurrence of order k of the form z[n] a function of z[n-1], ..., z[n-k]. Try:
#InvR((1+z[n-1])/z[n-2],2,z,n);
InvR:=proc(T,k,z,n) local gu,i:
gu:=solve(z[n]=T,z[n-k]):
subs({seq(z[n-i]=z[n-(k-i)],i=0..k-1)},gu):
end:


#PR(k,r,z,n): all the k-th order periodic recurrences of order k of the form z[n]=(1+...a[k-1]*z[n-k+1])/z[n-k] of period r. Try:
#PR(2,5,z,n);
PR:=proc(k,r,z,n) local f,a,i,F,var,eq,x,var1,var11,gu1,gu2,f1,Finv, r1,r2,fInv,gu,ka,ka1:
option remember:
var:={seq(a[i],i=1..k-1)}:
f:=(1+add(a[i]*z[n-i],i=1..k-1))/z[n-k]:
fInv:=InvR(f,k,z,n):

if r mod 2=0 then
 r1:=r/2:
 r2:=r/2:
else
 r1:=(r-1)/2:
 r2:=(r+1)/2:
fi:

F:=[seq(x[i],i=1..k)]:
Finv:=[seq(x[k-i+1],i=1..k)]:
gu1:=Traj(F,f,z,n,k,r1):
gu2:=Traj(Finv,fInv,z,n,k,r2):


gu1,gu2:

gu1:=[op(nops(gu1)-k+1..nops(gu1),gu1)]:
gu2:=[op(nops(gu2)-k+1..nops(gu2),gu2)]:
gu2:=[seq(gu2[k-i+1],i=1..k)]:
gu1,gu2:

eq:={seq(op(Coes(numer(gu1[i]-gu2[i]),F)),i=1..k)}:

eq:=Groebner[Basis](eq,plex(seq(a[i],i=1..k))):

var1:={solve({op(eq)},var)}:




if var1=NULL then
 RETURN(FAIL):
fi:

gu:={}:

for var11 in var1 do
 f1:=subs(var11,f):
 gu:=gu union {f1}:
od:

ka:=numtheory[divisors](r) minus {1, r}:
gu minus {seq(op(PR(k,ka1,z,n)),ka1 in ka)}:

end:



#PR1(k,r,z,n): all the k-th order periodic recurrences of order k of the form z[n]=(a[0]+...a[k-1]*z[n-k])/z[n-k] of period r. Try:
#PR1(2,5,z,n);
PR1:=proc(k,r,z,n) local f,a,i,F,var,eq,x,var1,var11,gu1,gu2,f1,Finv, r1,r2,fInv,gu,ka,ka1,lu:
option remember:
var:={seq(a[i],i=1..k)}:
f:=(1+add(a[i]*z[n-i],i=1..k))/z[n-k]:
fInv:=InvR(f,k,z,n):

if r mod 2=0 then
 r1:=r/2:
 r2:=r/2:
else
 r1:=(r-1)/2:
 r2:=(r+1)/2:
fi:

F:=[seq(x[i],i=1..k)]:
Finv:=[seq(x[k-i+1],i=1..k)]:
gu1:=Traj(F,f,z,n,k,r1):
gu2:=Traj(Finv,fInv,z,n,k,r2):


gu1,gu2:

gu1:=[op(nops(gu1)-k+1..nops(gu1),gu1)]:
gu2:=[op(nops(gu2)-k+1..nops(gu2),gu2)]:
gu2:=[seq(gu2[k-i+1],i=1..k)]:
gu1,gu2:

eq:={seq(op(Coes(numer(gu1[i]-gu2[i]),F)),i=1..k)}:

eq:=Groebner[Basis](eq,plex(seq(a[i],i=1..k))):

var1:={solve({op(eq)},var)}:




if var1=NULL then
 RETURN(FAIL):
fi:

gu:={}:

for var11 in var1 do
 f1:=subs(var11,f):
 gu:=gu union {f1}:
od:

ka:=numtheory[divisors](r) minus {1, r}:
gu:=gu minus {seq(op(PR1(k,ka1,z,n)),ka1 in ka)}:

lu:={}:

for gu1 in gu do
if {seq(coeff(numer(gu1),z[n-i],1),i=1..k-1)}<>{0} then
 lu:=lu union {gu1}:
fi:
od:
lu:

end:



#YafeR(T,z,n,k): Given an expression T whose numerator has the form (1+a[1]*z[n-1]+...+a[k]*z[n-k]) retuns it nicely. Try:
#YafeR(PR1(2,8,z,n)[2],z,n,2);
YafeR:=proc(T,z,n,k) local i,gu1,gu2:
gu1:=numer(T):
gu2:=denom(T):
gu1:=gu1-1:

(1+add(identify(evalf(coeff(gu1,z[n-i],1)))*z[n-i],i=1..k))/gu2:

end:

#Paper(K,R,z,n): Inputs a positive integers K, R, outputs  a paper with periodic recurrences of the from
#z[n]=(1+a1*z[n-1]+...+ak*z[n-k])/z[n-k] with periods from 2 to R (only that exist).
#It is inspired by Lozitsky's paper
#Fractional Linear Periodic Recurrences of orders two and three
#https://arxiv.org/abs/2310.13036. Try:
#Paper(2,12,z,n);

Paper:=proc(K,R,z,n) local t0,k,r,gu,gu1,a,i:
t0:=time():

print(`Searching for Lynnes-stype Periodic Ratioanl Recurrences inspired ny Evgeni Lozitsky`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(``):
print(`In the delightful paper: "Fractional Linear Periodic Recurrences of orders two and three", by Evgeni Lozitsky`):
print(`he came up with several rational peridoic recurrencs of orders 2 and 3 with periods up to 12. Here we emulate him`):
print(`and reproduce, in Maple, all his recurrences and try to find new ones, by extending the search for all recurrences of the form`):
print(`z[n]=(1+a1*z[n-1]+...+ak*z[n-k])/z[n-k] for k from 2 to`, K, `with periods from 2 to`,  R ):
print(``):

for k from 2 to K do
print(`--------------------------------`):
print(`Searching for recurrences of order`, k, ` of the form `):
print(``):
print(z[n]=(1+add(a[i]*z[n-i],i=1..k))/z[n-k]):
print(``):
print(`with periods up to R`):
print(``):
for r from 2 to R do
gu:=PR1(k,r,z,n):


if gu<>{} then
print(`--------------------------------`):
print(`There exist periodic Lynnes-style recurrences of order`, k, `and period`, r, `here they are `):
print(``):

for gu1 in gu do
print(z[n]=gu1):
print(``):
od:

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


for gu1 in gu do
lprint(z[n]=gu1):
print(``):
od:

fi:
od:
od:

print(`-------------------------------------`):
print(`This ends this article that took`, time()-t0, `seconds to generate `):

end: