######################################################################
##MetaFib: Save this file as    MetaFib                              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read MetaFib<Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Nathan Fox and Doron Zeilberger, Rutgers University ,   #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: April, 2015
 
print(`Created:  April 2015`):
print(` This is MetaFib `):
print(`It is one of the packages that accompany the article `):
print(`Automatic Search (and Automatic Proof) for C-finite Solutions to Hofstadter's Q-Meta-Fibonacci Recurrence `):
print(`by Nathan Fox and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

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(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  DD, IsLinear, FindPeriod, GuessPeriod, GuessRec `):
 print(` IsQL, LtoRmax, LtoRmin  `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: GuessGF,  Qg, SeqQg, SeqTg, SeqToQL, Tg, Tovim2,  Tovim2QL, Tovim2V, Tovim3, Tovim3QL,Tovim3V `):
 print(` `):



elif nops([args])=1 and op(1,[args])=DD then
print(`DD(L,m): inputs a sequence L and a positive integer m, outputs the sequence of differences L[i+m]-L[i] from i=1 to nops(L)-m`):
print(`try:`):
print(`DD([seq(i,i=1..30)],2);`):


elif nops([args])=1 and op(1,[args])=FindPeriod then
print(`FindPeriod(L,P): given a sequence L, finds the smallest period <=P or returns FAIL`):
print(`Try:`):
print(`FindPeriod([seq(1+3*(-1)^i,i=1..40)],10);`):

elif nops([args])=1 and op(1,[args])=GuessGF then
print(`GuessGF(L,t): inputs a sequence L, starting at L[1], and tries to guess`):
print(`a rational generating function in t. Try:`):
print(`GuessGF([1,1,1,1,1,1,1,1,1],t);`):


elif nops([args])=1 and op(1,[args])=GuessPeriod then
print(`GuessPeriod(Lis): Given a list, guesses an ultimate period. Try:`):
print(`GuessPeriod([2,3,seq((-1)^i+5,i=1..10)]);`):


elif nops([args])=1 and op(1,[args])=GuessRec then
print(`GuessRec(L): inputs a sequence L and tries to guess`):
print(`a recurrence operator with constant cofficients `):
print(`satisfying it. It returns the initial  values and the operator`):
print(` as a list. For example try:`):
print(`GuessRec([1,1,1,1,1,1]);`):
print(``):

elif nops([args])=1 and op(1,[args])=IsLinear then
print(`IsLinear(L,n): is the sequence L linear?`):

elif nops([args])=1 and op(1,[args])=IsQL then
print(`IsQL(f,t,M): inputs a rational generating function f, in the variable t, and a positive integer M,`):
print(`and outputs the smallest i<=M such (1-t^i)^2*f is a polynomial, together with the`):
print(`that polynomial or FAIL. Try: `):
print(`IsQL(1/(t^2-t+1),t,10); `):

elif nops([args])=1 and op(1,[args])=LtoRmax then
print(`LtoRmax(L): inputs a seqeunce L, and outputs the left-to-right maxima [place,value], try:`):
print(` LtoRmax([3,6,2,4,9,2]); `):

elif nops([args])=1 and op(1,[args])=LtoRmin then
print(`LtoRmin(L): inputs a seqeunce L, and outputs the left-to-right minima [place,value], try:`):
print(` LtoRmin([3,6,2,4,9,2]); `):

elif nops([args])=1 and op(1,[args])=Qg then
print(`Qg(L,n): Hofstadter recurrence with initial conditions L, try:`):
print(` Qg([1,1],10); `):

elif nops([args])=1 and op(1,[args])=SeqQg then
print(`SeqQg(L,N): The first N terms of the Hofstadter recurrence for Q with initial conditions L, try:`):
print(` SeqQg([1,1],100); `):

elif nops([args])=1 and op(1,[args])=SeqTg then
print(`SeqTg(L,N): The first N terms of the solutions to Tanny's  recurrence for T with initial conditions L, starting at n=0. Try:`):
print(` SeqTg([1,1,1],100); `):


elif nops([args])=1 and op(1,[args])=SeqToQL then
print(`SeqToQL(L,n,P): given a sequence L, converts it to a linear quasi-polynomial of quasi-period r<=P, say.`):
print(`The output is a list of initial values,INI, of length a multiple r, and a list of linear`):
print(`expressions (of length r), let's call it M, such that L[n] is equal to  M[i] if n mod r is i (i=1..r)`):
print(` and n>=nops(INI) `):
print(` Try: `):
print(` SeqToQL([seq(4+3*(-1)^i,i=1..200)],n,10); `):

elif nops([args])=1 and op(1,[args])=Tg then
print(`Tg(L,n): Tanny's recurrence with initial conditions L, try:`):
print(` Tg([1,1,1],10); `):

elif nops([args])=1 and op(1,[args])=Tovim2 then
print(`Tovim2(K,t,N): all the pairs with parts<=K that lead to rational functions generating functions, using N terms`):
print(`try:`);
print(`Tovim2(7,t,60);`):

elif nops([args])=1 and op(1,[args])=Tovim2QL then
print(`Tovim2QL(K,n,P): all the pairs with parts<=K that lead to Quasi-linear functions with at most period P`):
print(` expressed in terms of n. `):
print(`Try: `):
print(`Tovim2QL(4,n,30); `):

elif nops([args])=1 and op(1,[args])=Tovim2V then
print(`Tovim2V(K,t,N): a verbose version of Tovim2(K,t,N) (q.v.). `):
print(`try:`);
print(`Tovim2V(4,t,60);`):

elif nops([args])=1 and op(1,[args])=Tovim3 then
print(`Tovim3(K,t,N): all the triplets with parts<=K that lead to rational functions generating functions, using N terms`):
print(`try:`);
print(`Tovim3(3,t,60);`):

elif nops([args])=1 and op(1,[args])=Tovim3QL then
print(`Tovim3QL(K,n,P): all the triples with parts<=K that lead to Quasi-linear functions with at most period P`):
print(` expressed in terms of n. `):
print(`Try: `):
print(`Tovim3QL(4,n,12); `):

elif nops([args])=1 and op(1,[args])=Tovim3V then

print(`Tovim3V(7,t,60);`):
print(`Tovim2V(K,t,N): a verbose version of Tovim3(K,t,N) (q.v.). `):
print(`try:`);
print(`Tovim3V(3,t,60);`):

else
print(`There is no ezra for`,args):
fi:
 
end:


#Golumb
ez:=proc(): print(` Tovim2(K,t,N), Tovim3(K,t,N) `): end:

#GuessRec1(L,d): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients of order d
#satisfying it. It returns the initial d values and the operator
# as a list. For example try:
#GuessRec1([1,1,1,1,1,1],1);
GuessRec1:=proc(L,d) local eq,var,a,i,n:

if nops(L)<=2*d+2 then
 print(`The list must be of size >=`, 2*d+3 ):
 RETURN(FAIL):
fi:

var:={seq(a[i],i=1..d)}:

eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):
fi:

end:


#GuessRec(L): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients 
#satisfying it. It returns the initial  values and the operator
# as a list. For example try:
#GuessRec([1,1,1,1,1,1]);
GuessRec:=proc(L) local gu,d:

for d from 1 to trunc(nops(L)/2)-2 do
 gu:=GuessRec1(L,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:

end:


#Qg(L,n): Hofstadter recurrence with initial conditions L, try:
##Qg([1,1],10);
Qg:=proc(L,n) option remember:

if n<1 then
 0:
elif n<=nops(L) then
 L[n]:

elif
(Qg(L,n-1)<=0 or Qg(L,n-2)<=0 or Qg(L,n-1)=FAIL or Qg(L,n-2)=FAIL) then
 FAIL:

else
Qg(L,n-Qg(L,n-1))+Qg(L,n-Qg(L,n-2)):
fi:

end:

SeqQg:=proc(L,N) local n,lu,lu1:
lu:=[]:
for n from 1 to N do
lu1:=Qg(L,n):

if lu1=FAIL then
  RETURN(lu):
else
 lu:=[op(lu),lu1]:
fi:

od:

lu:

end:

#FindN(L,a,d): inputs initial conditions for the generalized Hofstadter sequence, L, and positive integers
#outputs the set of all triples [i,j,Rec], such that the subsequence Qg(i*n+j) is C-finite together with its recurrence
#Try:
#FindN([3,2,1],3,5);
FindN:=proc(L,a,d) local gu,i,j,gu1,lu,vu,n:
gu:=SeqQg(L, (2*d+10)*(a+1)):

lu:={}:

for i from 1 to a do
 for j from 0 to i-1 do
  gu1:=[seq(gu[n*i+j],n=1..2*d+10)]:
  vu:=GuessRec(gu1):
   if vu<>FAIL then
     lu:=lu union {[i,j,vu]}:
   fi:
 od:
od:

lu:

end:


#GuessGF(L,t): inputs a sequence L, starting at L[1], and tries to guess
#a rational generating function in t. Try:
#GuessGF([1,1,1,1,1,1,1,1,1],t);
GuessGF:=proc(L,t) local gu,f,Den,i:
gu:=GuessRec(L):
if gu=FAIL then
  RETURN(FAIL):
fi:

gu:=gu[2]:

Den:=1-add(gu[i]*t^i,i=1..nops(gu)):
f:=expand(add(L[i]*t^i,i=1..nops(L))*Den):
f:=add(coeff(f,t,i)*t^i,i=0..nops(gu)+1):
f:=f/Den:
if [seq(coeff(taylor(f,t=0,nops(L)+1),t,i),i=1..nops(L))]<>L then
 FAIL:
else
 f:
fi:

end:


#Tovim2(K,t,N): all the pairs with parts<=K that lead to rational functions generating functions, using N terms
#try:
#Tovim2(7,t,60);
Tovim2:=proc(K,t,N) local K1,K2,L,gu,mu:

mu:={}:

for K1 from 1 to K do
for K2 from 1 to K do
 L:=[K1,K2]:
  gu:=GuessGF(SeqQg(L,N),t):

  if gu<>FAIL then
    mu:=mu union {[L,gu]}:
  fi:
od:
od:
mu:
end:

#Tovim3(K,t,N): all the triples with parts<=K that lead to rational functions generating functions, using N terms
#try:
#Tovim3(7,2,t,60);
Tovim3:=proc(K,t,N) local K1,K2,K3,L,gu,mu:

mu:={}:

for K1 from 1 to K do
for K2 from 1 to K do
for K3 from 1 to K do
 L:=[K1,K2,K3]:
  gu:=GuessGF(SeqQg(L,N),t):

  if gu<>FAIL then
    mu:=mu union {[L,gu]}:
  fi:
od:
od:
od:

mu:
end:
 



#Tovim2V(K,t,N): Verbose version of Tovim2(K,t,N)
#try:
#Tovim2V(7,t,60);
Tovim2V:=proc(K,t,N) local gu,t0,i:
t0:=time():

gu:=Tovim2(K,t,N):

print(`All the Solutions to Hofsdater's Q-recurrence with initial 2 initial conditions between 1 and`, K ):
print(`That happen to have Rational Generating Functions`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Recall that the famous (or, more accurately, infamous), Hofsdater's Q-sequence, Sequence A005185 in the OEIS`):
print(`is defined by the intial conditions Q(1)=1, Q(2)=1, and for n>2 by `):
print(`Q(n)=Q(n-Q(n-1))+Q(n-Q(n-2)) `):
print(`It seems to be completely "chaotic". Nevertheless, if instead of the initial conditions [1,1], one  takes the`):
print(`following one, one gets sequences whose generating functions are rational.`):

for i from 1 to nops(gu) do

print(`INI= `, gu[i][1], `Generating Function=`, gu[i][2]):
print(``):
print(`and in Maple notation: `):
lprint( gu[i][2] ):
od:

print(`To sum up, here are the succesful ones, in Maple notation. `):
print(``):
lprint(gu):
print(``):
print(`-----------------------------------------------------------------`):
print(``):
print(`This ends this article, that took`, time()-t0, `seconds. `):
end:

#Tovim3V(K,t,N): Verbose version of Tovim3(K,t,N)
#try:
#Tovim3V(7,t,60);
Tovim3V:=proc(K,t,N) local gu,t0,i:
t0:=time():

gu:=Tovim3(K,t,N):

print(`All the Solutions to Hofsdater's Q-recurrence with initial 2 initial conditions between 1 and`, K ):
print(`That happen to have Rational Generating Functions`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Recall that the famous (or, more accurately, infamous), Hofsdater's Q-sequence, Sequence A005185 in the OEIS`):
print(`is defined by the intial conditions Q(1)=1, Q(2)=1, and for n>2 by `):
print(`Q(n)=Q(n-Q(n-1))+Q(n-Q(n-2)) `):
print(`It seems to be completely "chaotic". Nevertheless, if instead of the initial conditions [1,1], one  takes the`):
print(`following one, one gets sequences whose generating functions are rational.`):

for i from 1 to nops(gu) do

print(`INI= `, gu[i][1], `Generating Function=`, gu[i][2]):
print(``):
print(`and in Maple notation: `):
lprint( gu[i][2] ):
od:

print(`To sum up, here are the succesful ones, in Maple notation. `):
print(``):
lprint(gu):
print(``):
print(`-----------------------------------------------------------------`):
print(``):
print(`This ends this article, that took`, time()-t0, `seconds. `):
end:

#DD(L,m): inputs a sequence L and a positive integer m, outputs the sequence of differences L[i+m]-L[i] from i=1 to nops(L)-m
#try:
#DD([seq(i,i=1..30)],2);
DD:=proc(L,m) local i:

[seq(L[i+m]-L[i],i=1..nops(L)-m)]:

end:

#GuessPeriod1(Lis): Given a list, guesses an ultimate period. Try:
#GuessPeriod1([2,3,seq((-1)^i+5,i=1..10)]);
GuessPeriod1:=proc(Lis)
local s0,p,n,i,j,gu:
n:=nops(Lis)-4:
for s0 from 1 to trunc(n/3) do
for p from 1 to trunc((n-s0)/3) do
 if {seq(nops({seq(Lis[s0+p*i+j],i=0..trunc((n-j-s0)/p))}),j=1..p)}={1} 
   then
 gu:=[[op(1..s0,Lis)],[op(s0+1..s0+p,Lis)]]:
  if nops(gu[1])=1 and nops(gu[2])>=1 and gu[1][1]=gu[2][1] then
    gu:=[[],gu[2]]:
  fi:

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



#GuessPeriod(Lis): Given a list, guesses an ultimate period. Try:
#GuessPeriod1([2,3,seq((-1)^i+5,i=1..10)]);
GuessPeriod:=proc(Lis) local gu:
gu:=GuessPeriod1(Lis):

if gu[1][-1]=gu[2][-1] then
 gu:=[[op(1..nops(gu[1])-1,gu[1])],[gu[1][-1],op(1..nops(gu[2])-1,gu[2])] ]:
fi:
gu:

end:



#IsLinear(L,n): is the sequence L linear?
IsLinear:=proc(L,n) local i,mu:

if nops(convert(L,set))=1 then
  RETURN(true,L[1]):
fi:

mu:={seq(L[i+1]-L[i],i=1..nops(L)-1)}:

if nops(mu)=1 then
  RETURN(true,L[1]+mu[1]*(n-1)):
fi:
false,0:
end:





#FindPeriod(L,P): given a sequence L, finds the smallest period <=P or returns FAIL
#Try:
#FindPeriod([seq(1+3*(-1)^i,i=1..40)],10);

FindPeriod:=proc(L,P) local L1,p,n,i,j:
if nops(L)<60 then
 print(`Sequnence must be at least of length 60`):
 RETURN(FAIL):
fi:


L1:=[op(30..nops(L),L)]:

for p from 1 to min(P,nops(L1)/4)  do


 if {seq(IsLinear([seq(L1[p*i+j],i=1..(nops(L1)-j)/p)],n)[1],j=1..p)}={true} then
  RETURN(p):
 fi:
od:

FAIL:
end:


#SeqToQL(L,n,P): given a sequence L, converts it to a linear quasi-polynomial of quasi-period r<=P, say.
#The output is a list of initial values,INI, of length a multiple r, and a list of linear
#expressions (of length r), let's call it M, such that L[n] is equal to  M[i] if n mod r is i (i=1..r)
#and n>=nops(INI)
#Try:
#SeqToQL([seq(4+3*(-1)^i,i=1..200)],n,10);
SeqToQL:=proc(L,n,P) local L1,p,M,hat,L1a,ku,i,j,lu,L2:

p:=FindPeriod(L,P):

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

hat:=trunc(30/p)+1:

L1:=[op(hat*p+1..nops(L),L)]:

M:=[]:

for i from 1 to p do
 L1a:=[seq(L1[j*p+i],j=1..trunc((nops(L1)-i)/p) )]:
 ku:=IsLinear(L1a,n):
   if not ku[1] then
      RETURN(FAIL):
  else
  lu:=expand(subs(n=(n-i)/p,ku[2])):
   lu:=expand(subs(n=n-hat*p,lu)):
   M:=[op(M),lu]:
  fi:
od:

L2:=[]:

for i from 0 to trunc(nops(L)/p)-1 do
for j from 1 to p do
L2:=[op(L2),subs(n=i*p+j,M[j])]:
od:
od:

for i from nops(L2) by -1 to 1 while L2[i]=L[i] do od:

[[op(1..i,L)],M]:



end:


#Tovim2QL(K,n,P): all the pairs with parts<=K that lead to Quasi-linear functions with at mostperiod P
#expressed in terms of n.
#Try:
#Tovim2QL(7,n,30);
Tovim2QL:=proc(K,n,P) local N,K1,K2,L,gu,mu,vu:

N:=100+P*10:
mu:={}:

for K1 from 1 to K do
for K2 from 1 to K do
 L:=[K1,K2]:
 vu:=SeqQg(L,N):

if nops(vu)>60 then
 gu:=SeqToQL(vu,n,P):

  if gu<>FAIL then
    mu:=mu union {[L,gu]}:
  fi:
fi:
od:
od:
mu:
end:


#Tovim3QL(K,n,P): all the triples with parts<=K that lead to Quasi-linear functions with at mostperiod P
#expressed in terms of n.
#Try:
#Tovim3QL(4,n,30);
Tovim3QL:=proc(K,n,P) local N,K1,K2,K3,L,gu,mu,vu:

N:=100+P*10:
mu:={}:

for K1 from 1 to K do
for K2 from 1 to K do
for K3 from 1 to K do
 L:=[K1,K2,K3]:
 vu:=SeqQg(L,N):


if nops(vu)>60 and K3<>SeqQg([K1,K2],3)[3] then
 gu:=SeqToQL(vu,n,P):

  if gu<>FAIL then
    mu:=mu union {[L,gu]}:
  fi:
fi:
od:
od:
od:

mu:
end:

#LtoRmax(L): inputs a seqeunce L, and outputs the left-to-right maxima [place,value], try:
#LtoRmax([3,6,2,4,9,2]);
LtoRmax:=proc(L) local gu,ma,i:

if L=[] then
 RETURN([]):
fi:

ma:=L[1]:
gu:=[[1,ma]]:

for i from 2 to nops(L) do
 if L[i]>ma then
   ma:=L[i]:
  gu:=[op(gu),[i,ma]]:
 fi:
od:

gu:

end:


#LtoRmin(L): inputs a seqeunce L, and outputs the left-to-right minima [place,value], try:
#LtoRmin([3,6,2,4,9,2]);
LtoRmin:=proc(L) local gu,ma,i:

if L=[] then
 RETURN([]):
fi:

ma:=L[1]:
gu:=[[1,ma]]:

for i from 2 to nops(L) do
 if L[i]<ma then
   ma:=L[i]:
  gu:=[op(gu),[i,ma]]:
 fi:
od:

gu:

end:


#IsQL(f,t,M): inputs a rational generating function f, in the variable t, and a positive integer M,
#and outputs the smallest i<=M such (1-t^i)^2*f is a polynomial, together with the
#that polynomial or FAIL. Try:
#IsQL(1/(t^2-t+1),t,10);
IsQL:=proc(f,t,M) local i,lu:

for i from 1 to M do
 lu:=normal(f*(1-t^i)^2):
 if degree(denom(lu),t)=0 then
   RETURN([i,lu]):
 fi:
od:

FAIL:

end:



#Tg(L,n): Tanny's sequence  with initial conditions L, starting at n=0. Try:
##Tg([1,1,1],10);
Tg:=proc(L,n) option remember:

if n<0 then
 0:
elif n<=nops(L)-1 then
 L[n+1]:

elif
(Tg(L,n-1)<=-1 or Tg(L,n-2)<=-2 or Tg(L,n-1)=FAIL or Tg(L,n-2)=FAIL) then
 FAIL:

else
Tg(L,n-1-Tg(L,n-1))+Tg(L,n-2-Tg(L,n-2)):
fi:

end:

SeqTg:=proc(L,N) local n,lu,lu1:
lu:=[]:
for n from 0 to N do
lu1:=Tg(L,n):

if lu1=FAIL then
  RETURN(lu):
else
 lu:=[op(lu),lu1]:
fi:

od:

lu:

end:


#T1(n):Tanny's original sequence using the fast recurrence (2.2), (2.3) in  his article, Discrete Math 105(1992), 227-239
T1:=proc(n) local x,i:
option remember:
if n=0 or n=1 or n=2 then
 RETURN(1):
fi:

x:=trunc(log[2](n))-1:

i:=n-2^(x+1):

if i<2^x then
T1(2^x+i-1) + 2^(x-1):
else
i:=i-2^x:
T1(2^x+i)+2^x:
fi:

end: