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


print(`First posted: Aug. 4, 2010: tested on Maple 13 `):
print(`Version of Aug. 4, 2010  `):
print():
print(`This is DALIA. It is a Maple package`):
print(`to investigate the recurrences in Aviezri S. Fraenkel's and Dalia Krieger's article  `):
print(`The Structure of Complementary sets of integers: a 3-shift theorem`):
print(`Internat. J. Pure and Appl. Math. 10 (2004) 1--49 `):
print(` available on line from: http://www.wisdom.weizmann.ac.il/~fraenkel/Papers/tsopaper1.pdf `):

print():
print(`The most current version of this Maple package is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/DALIA .`):
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(`AnS, AnSFromFK,Apn, F,Findpq,  Irregulars, mex,SnS, `):
print(` tton,   RecSeq, TovW, Wn, WnS `):
else
ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then

print(` DALIA: A Maple package for  studying the Fraenkel-Krieger phenomenon`):
print(`as described in their article`):
print(`Aveizri S. Fraenkel and Dalia Krieger, `):
print(`The structure of complementary sets of integers: A 3-shift theorem`):
print(`Inter. J. of Pure and Applied Math. 10 (2004) 1-49, available on-line`):
print(`http://www.wisdom.weizmann.ac.il/~fraenkel/Papers/tsopaper1.pdf`):
print(``):
print(`The MAIN procedures are`):
print(`An, AnFK, CheckFK, FK, FKv `):
print(` `):

elif nargs=1 and args[1]=Apn then
print(` Apn(c,n): The n-th term of the Fraenkel-Krieger recurrence`):
print(`for n0=1 and X={} given by the nice formula `):
print(`For example, try:`):
print( ` Apn(2,618); `):

elif nargs=1 and args[1]=An then
print(` An(n0,X,c,n): The n-th term of the Fraenkel-Krieger recurrence`):
print(`parametrized by n0,X, and c. `):
print(` This is the slow version, using the definition `):
print(`For example, try:`):
print( ` An(6,{1,5},2,618); `):

elif nargs=1 and args[1]=AnS then
print(` AnS(n0,X,c,N): The first N terms of An(n0,X,c,n)`):
print(`For example, try:`):
print( ` AnS(6,{1,5},2,300); `):

elif nargs=1 and args[1]=AnFK then
print(` AnFK(n0,X,c,n): The same as An(n0,X,c,n) `):
print(`BUT much faster, using the amazing Fraenkel-Krieger algorithm`):
print(`Once the initial investment of finding FK(n0,c,X) has been spent`):
print(`it is poly time (w.r.t. to bit-size) rather than exp. time `):
print(`For example, try:`):
print( ` AnFK(6,{1,5},2,618); `):


elif nargs=1 and args[1]=CheckFK then
print(` CheckFK(n0,X,c,N): empirically checks that the Fraenkel-Krieger`):
print(`algorithm is correct for the given parameters and for the`):
print(` first N terms`):
print(`For example, try:`):
print(`CheckFK(6,{1,5},2,1000);`):

elif nargs=1 and args[1]=F then
print(`F(w,c): applies the morphism in Definition 5`):
print(`For example, try:`):
print(`F([1,2,1],3);`):


elif nargs=1 and args[1]=Findpq then
print(`Findpq(c,n0,X): finds the places p and q `):
print(`that are synchronizing and for which all the pairs`):
print(`are adjacent `):
print(`It is promised in the proof of Theorem 1 `):
print(`For example, try:`):
print(`Findpq(2,6,{1,5});`):


elif nargs=1 and args[1]=FK then
print(` FK(n0,c,X): implements the Fraenkel-Krieger algorithm`):
print(`described in their beautiful paper.`):
print(`The inputs are: `):
print(`c: a positive integer `):
print(`n0: a positive integer `):
print(`X: a finite set (of course it is finite, there are no other kinds)`):
print(` f positive integers `):
print(``):
print(`The output is a list of length 4 that completely describes`):
print(`the sequence described by the Fraenkel-Krieger recurrence`):
print(`A(n)=mex({A(i),i=1..n-1} union {A(i)+c*i,i=1..n-1} union X )`):
print(`for n>=n0 `):
print(``):

print(`The first component is the "noisy beginning" `):
print(``):
print(`the second component is the regular shift, gamma `):
print(``):
print(`the third component is the list of sequenes for the places`):
print(`of the irregular shifts, given in terms of their initial conditons`):
print(`that start with the LOWER irregular shift, gamma-1 .`):
print(`As proved by Aviezri Fraenkel and Dalia Krieger, `):
print(`it satisfies a recurrence n[j]=c*n[j-1]+n[j-2] , so only`):
print(`[n[0],n[1]] has to be given `):
print(``):
print(`the fourth component is the list of sequenes for the places`):
print(`of the irregular shifts, given in terms of their initial conditons`):
print(`that start with the HIGHER irregular shift, gamma+1 .`):
print(`As proved by Aviezri Fraenkel and Dalia Krieger, `):
print(`it satisfies a recurrence n[j]=c*n[j-1]+n[j-2] , so only`):
print(`[n[0],n[1]] has to be given `):
print(`For example, try:`):
print(`FK(6,2,{1,5});`):


elif nargs=1 and args[1]=FKv then
print(` FKv(n0,c,X): A verbose rendition of`):
print(` FK(n0,c,X)`):
print(`For example, try:`):
print(`FKv(6,,2{1,5});`):


elif nargs=1 and args[1]=Irregulars then
print(`Irregulars(n0,X,c,N): the low-irregular indices`):
print(`and the high ones <=N`):
print(`For example, try:`):
print(`Irregulars(6,{1,5},2,1000);`):


elif nargs=1 and args[1]=mex1 then
print(` mex1(X): the smallest positive integer not in X`):
print(`For example, try:`):
print( ` mex1({1,2,5}; `):


elif nargs=1 and args[1]=RecSeq then
print(`RecSeq(m0,m1,c,N): all the terms of the sequence`):
print(`obeying the recurrence n[j]=c*n[j-1]+n[j-2] that are`):
print(`less than N`);
print(`For example, try:`):
print( ` RecSeq(1,1,1,1000); `):


elif nargs=1 and args[1]=Sn then
print(` Sn(n0,X,c,n): AnS(1,{},c,n) - AnS(n0,X,c,n) `):
print(`For example, try:`):
print( ` Sn(6,{1,5},2,300); `):

elif nargs=1 and args[1]=SnS then
print(` SnS(n0,X,c,N): The first N terms of Sn(n0,X,c,n)`):
print(`For example, try:`):
print( ` SnS(6,{1,5},2,300); `):

elif nargs=1 and args[1]=TovW then
print(`TovW(v1,v2): Given a pair of words in {1,2}`):
print(`decides whether [1,2] is followed immediately by a [2,1]`);
print(`and vice versa. For example, try:`):
print(`TovW([1,2],[2,1]);`):


elif nargs=1 and args[1]=tton then
print(`tton(c,n0,X,t): Given a Fraenkel-Krieger recurrence in terms`):
print(`of the parameters c, n0, X, and an index t>=max(n0,X)`):
print(` finds the matching n (as in Lemma 2). `):
print(` For example, try: `):
print(` tton(2,6,{1,4},10); `):


elif nargs=1 and args[1]=Wn then
print(` Wn(n0,X,c,n): An(n0,X,c,n+1)-An(n0,X,c,n) `):
print(`For example, try:`):
print( ` Wn(6,{1,5},2,618); `):

elif nargs=1 and args[1]=WnS then
print(` WnS(n0,X,c,N): The first N terms of Wn(n0,X,c,n)`):
print(`For example, try:`):
print( ` WnS(6,{1,5},2,300); `):

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

fi:


end:





Digits:=100:

mex1:=proc(S) local i:
for i from  1 while member(i,S) do od:
i:
end:

An:=proc(n0,X,c,n) local i:
option remember:
mex1({seq(An(n0,X,c,i),i=n0..n-1),seq(An(n0,X,c,i)+c*i,i=n0..n-1)}
union X):
end:


Wn:=proc(n0,X,c,n)
option remember:
An(n0,X,c,n+1)-An(n0,X,c,n):
end:
AnS:=proc(n0,X,c,N) local n:
option remember:
[seq(An(n0,X,c,n),n=1..N)]:
end:

WnS:=proc(n0,X,c,N) local i,gu:
gu:=AnS(n0,X,c,N):
[seq(gu[i+1]-gu[i],i=1..nops(gu)-1)]:
end:


SnS:=proc(n0,X,c,N) local gu1,gu2,i:
gu1:=AnS(1,{},c,N): 
gu2:=AnS(n0,X,c,N): 
[seq(gu1[i]-gu2[i],i=1..N)]:
end:

F:=proc(w,c) local i,w1:
w1:=[]:
for i from 1 to nops(w) do
 if w[i]=1 then
    w1:=[op(w1),1$(c-1),2]:
 elif w[i]=2 then
    w1:=[op(w1),1$c,2]:
  else
   ERROR(`Something is wrong`):
 fi:
od:
w1:
end:


Lemma2:=proc(n0,X,c,N) local n,gu,u,i,t,n1,u1,gu1:
for n from n0 while An(n0,X,c,n)<=max(op(X)) do od:
t:=n:


for n from t while An(n0,X,c,n)<>An(n0,X,c,t)+c*t+1 do od:
 u:=[seq(An(n0,X,c,n1),n1=t..n)]:
 u:=[seq(u[i+1]-u[i],i=1..nops(u)-1)]:

u1:=u:
gu:=[]:
for i from 1 while nops(gu)<N do
gu:=[op(gu),op(u1)]:
u1:=F(u1,c): 
od:

gu1:=W(n0,X,c,N):

gu:=[op(1..nops(gu1),gu)]:
evalb(gu1=gu):
end:


#TovW(v1,v2): Given a pair of words in {1,2}
#decides whether [1,2] is followed immediately by a [2,1]
#and vice versa. For example, try:
#TovW([1,2],[2,1]);
TovW:=proc(v1,v2) local v1a,v2a:

if nops(v1)<>nops(v2) then
 RETURN(FAIL):
fi:

v1a:=v1:
v2a:=v2:

while nops(v1a)>1 do

if v1a[1]=v2a[1] then
 v1a:=[op(2..nops(v1a),v1a)]:
 v2a:=[op(2..nops(v2a),v2a)]:

 elif v1a[1]=1 and v2a[1]=2 then
 
  if not (v1a[2]=2 and v2a[2]=1) then
    RETURN(false):
   else
    v1a:=[op(3..nops(v1a),v1a)]: 
    v2a:=[op(3..nops(v2a),v2a)]:
  fi:

 elif v1a[1]=2 and v2a[1]=1 then

  if not (v1a[2]=1 and v2a[2]=2) then
    RETURN(false):
   else
    v1a:=[op(3..nops(v1a),v1a)]: 
    v2a:=[op(3..nops(v2a),v2a)]:
  fi:

else
   ERROR(`Something is wrong`):
fi:

od:

true:

end:


#tton(c,n0,X,t): Given a F-K recurrence in terms
#of the parameters c, n0, X, and an index t>=max(n0,X)
#finds the matching n (as in Lemma 2).
#For example, try:
#tton(2,100,3,{1,4},10);
tton:=proc(c,n0,X,t) local n:
option remember:

for n from t+1 while An(n0,X,c,n) <>An(n0,X,c,t)+c*t+1 do od:
n:


end:


Findpq:=proc(c,n0,X) local n,t,v1,v2,i,m0:

for m0 from n0 while An(n0,X,c,m0)<=max(op(X)) do od:

for  t from m0 do
 if tton(c,n0,X,t)=tton(c,1,{},t) then
     n:=tton(c,n0,X,t):
     v1:=[seq(Wn(n0,X,c,i),i=t..n-1)]:
     v2:=[seq(Wn(1,{},c,i),i=t..n-1)]:

       if TovW(v1,v2) then
         RETURN([t,n]):
       fi:
  fi:
od:

end:




#FK(n0,c,X): finds the noisy beginning, of the
#An sequence, the regular shift, gammar, and the starting
#places for sequenes of irregular shifts, starting with
#gamma-1 followed by those starting with gamma+1
#For example, try:
#FK(6,2,{1,5});
FK:=proc(n0,c,X) local p,q,pq,v1,v2,u1,u2,u1F,u2F,gu1,gu2,dalia,dalia1,
d1, luLow,luHigh,i, luLowA,luHighA:
option remember:
pq:=Findpq(c,n0,X):
p:=pq[1]: q:=pq[2]:


v1:=WnS(n0,X,c,q):   
v2:=WnS(1,{},c,q):   
u1:=[op(p..q-1,v1)]:
u2:=[op(p..q-1,v2)]:
v1:=[op(1..p-1,v1)]:
v2:=[op(1..p-1,v2)]:

u1F:=F(u1,c):
u2F:=F(u2,c):



gu1:=[op(v1),op(u1),op(u1F)]:
gu2:=[op(v2),op(u2),op(u2F)]:


if WnS(n0,X,c,nops(gu1)+1)<>gu1 then
 RETURN(FAIL):
fi:


if WnS(1,{},c,nops(gu2)+1)<>gu2 then
 RETURN(FAIL):
fi:


dalia:=SnS(n0,X,c,nops(gu1)+10):
dalia1:={op(p..nops(dalia),dalia)}:



if not ( nops(dalia1)=3 or nops(dalia1)=1 ) then
 print(`The Fraenkel-Krieger theorem is wrong!`):

 RETURN(FAIL):
fi:


d1:=min(op(dalia1)):


if not (dalia1={d1} or dalia1={d1,d1+1,d1+2}) then
 print(`The Fraenkel-Krieger theorem is wrong!`):
 RETURN(FAIL):
fi:


if nops(dalia1)=3 then
 d1:=d1+1:
fi:

luLow:=[]:
luHigh:=[]:

for i from p to q-1 do
 if dalia[i]=d1-1 then
   luLow:=[op(luLow),i]:
 elif dalia[i]=d1+1 then
    luHigh:=[op(luHigh),i]:
  fi:
od:


luLowA:=[]:
luHighA:=[]:

for i from q to nops(gu1) do
 if dalia[i]=d1-1 then
   luLowA:=[op(luLowA),i]:
 elif dalia[i]=d1+1 then
    luHighA:=[op(luHighA),i]:
  fi:
od:



if nops(luLow)<>nops(luHighA) then
 print(`The Fraenkel-Krieger theorem is wrong!`):
 RETURN(FAIL):
fi:


if nops(luHigh)<>nops(luLowA) then
 print(`The Fraenkel-Krieger theorem is wrong!`):
 RETURN(FAIL):
fi:



[AnS(n0,X,c,p-1),d1,
[seq([luLow[i],luHighA[i]],i=1..nops(luLow))],
[seq([luHigh[i],luLowA[i]],i=1..nops(luHigh))]]:

end:



#RecSeq(m0,m1,c,N): all the terms of the sequence
#obeying the recurrence n[j]=c*n[j-1]+n[j-2] that are
#less than N
RecSeq:=proc(m0,m1,c,N) local gu:
gu:=[m0,m1]:
while gu[nops(gu)]<=N do
 gu:=[op(gu),c*gu[nops(gu)]+gu[nops(gu)-1]]:
od:
gu:
end:






#AnSFromFK(n0,X,c,N): Like AnS but using the Fraenkel-Krieger algorithm
#just for checking
AnSFromFK:=proc(n0,X,c,N) local lu,gu,i,T,j,cu:
lu:=FK(n0,c,X):
if N<=nops(lu[1]) then
 RETURN([op(1..N,gu[1])]):
fi:

gu:=AnS(1,{},c,N):

gu:=[op(AnS(n0,X,c,nops(lu[1]))),seq(gu[i]-lu[2],i=nops(lu[1])+1..N)]:

for i from 1 to nops(gu) do
T[i]:=gu[i]:
od:

for i from 1 to nops(lu[3]) do
 cu:=lu[3][i]:
 cu:=RecSeq(op(cu),c,N):
  if cu[nops(cu)]>N then
   cu:=[op(1..nops(cu)-1,cu)]:
  fi:

  for j from 1 to nops(cu) do
   T[cu[j]]:=  T[cu[j]]+(-1)^(j-1):
  od:
od:


for i from 1 to nops(lu[4]) do
 cu:=lu[4][i]:
 cu:=RecSeq(op(cu),c,N):
  if cu[nops(cu)]>N then
   cu:=[op(1..nops(cu)-1,cu)]:
  fi:

  for j from 1 to nops(cu) do
   T[cu[j]]:=  T[cu[j]]+(-1)^(j):
  od:
od:

 
[seq(T[i],i=1..N)]:
end:




CheckFK:=proc(n0,X,c,N):
evalb(AnSFromFK(n0,X,c,N) =AnS(n0,X,c,N) ):
end:



#Apn(c,n): computing A_c'(n) via the formula
Apn:=proc(c,n) local alpha:
alpha:=evalf((2-c+sqrt(c^2+4))/2):
trunc(alpha*n):
end:




#Irregulars(n0,X,c,N): the low-irregular indices
#and the high ones <=N
Irregulars:=proc(n0,X,c,N) local Namukh, Gavoa,cu,i,j,lu:
lu:=FK(n0,c,X):

Namukh:={}:
Gavoa:={}:

for i from 1 to nops(lu[3]) do
 cu:=lu[3][i]:
 cu:=RecSeq(op(cu),c,N):
  if cu[nops(cu)]>N then
   cu:=[op(1..nops(cu)-1,cu)]:
  fi:

  for j from 1 to nops(cu) do
    if j mod 2 =1 then
     Namukh:=Namukh union {cu[j]}:
    else
      Gavoa:=Gavoa union {cu[j]}:
    fi:
   od:

od:


for i from 1 to nops(lu[4]) do
 cu:=lu[4][i]:
 cu:=RecSeq(op(cu),c,N):
  if cu[nops(cu)]>N then
   cu:=[op(1..nops(cu)-1,cu)]:
  fi:

  for j from 1 to nops(cu) do
    if j mod 2 =1 then
     Gavoa:=Gavoa union {cu[j]}:
    else
      Namukh:=Namukh union {cu[j]}:
    fi:
  od:

od:

Namukh,Gavoa:
end:

AnFK:=proc(n0,X,c,n) local lu,gu,mu:
option remember:

lu:=FK(n0,c,X):
gu:=Irregulars(n0,X,c,n+3):

if n<=nops(lu[1]) then
 RETURN(An(n0,X,c,n)):
fi:


mu:=Apn(c,n)-lu[2]:

if member(n,gu[1]) then
mu+1:
elif member(n,gu[2]) then
mu-1:
else
mu:
fi:


end:





#FKv(n0,c,X): Verbose form of FK
FKv:=proc(n0,c,X) local gu,alpha,j,i1,r,n,t0:

alpha:=(2-c+sqrt(c^2+4))/2:

gu:=FK(n0,c,X):
print(`The Fraenkel-Krieger Solution to a certain mex recurrence`):
print(``):
print(` By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Define the sequence A(n) for n>=n0, with n0=`, n0 ):
print(`in terms of the following recurrence`):
print(``):
print(`A(n)=mex_1({A(i);i=n0..n-1)} Union {A(i)+c*i;i=n0..n-1} Union X ) `):
print(`where c=`, c, `n0=`, n0, `(as above) and `):
print(`where the set X is`, X, `and mex_1(S) is the smallest pos. integer`):
print(`NOT in S, for example mex({1,2,3,5})=4, mex({3,4,9})=1. `):

print(``):
print(`We have the following very fast way of computing A(n),`):
print(`without using the recurrence (it is polynomial time vs.`):
print(`exponential time (in bit-size, i.e. in log(n),`):
print(` the number of digits of n))`):
print(``):
if nops(gu[1])>0 then
print(`The first`, nops(gu[1])-n0+1 , ` members of the sequence are `):
print(`starting with n=`, n0, ` all the way to n=`, nops(gu[1]) ):
print(op(n0..nops(gu[1]),gu[1])):
fi:

if gu[3]=[] and gu[4]=[] then

print(`Starting at n=`, nops(gu[1])+1):
print(`A(n) is given by the succinct expression`):
print(IntegerPart(alpha*n)+gu[2]):

fi:

if gu[3]<>[] then
print(`Starting at n=`, nops(gu[1])+1):
print(`A(n) is given by the succinct expression`):
print(IntegerPart(alpha*n)+gu[2]+r(n)):
print(`where r(n) is usually 0, and otherwise it is 1 or -1`):

print(`as follows: `):

print(`consider the following`, nops(gu[3]), `sequences `):
print(`each satisfying the recurrence `):
print(n[j]=c*n[j-1]+n[j-2]):
print(`with the following initial conditions`):

for i1 from 1 to nops(gu[3]) do
 print(n[0]=gu[3][i1][1], n[1]=gu[3][i1][2]):
od:

print(`for each `, n[j], `j>=0 , with j even, r(n):=1 `):
print(`for each `, n[j], `j>=1 , with j odd, r(n):=-1 `):

fi:

if gu[4]<>[] then
print(`Also consider the following`, nops(gu[4]), `sequences `):
print(`each satisfying the recurrence `):
print(n[j]=c*n[j-1]+n[j-2]):
print(`with the following initial conditions`):

for i1 from 1 to nops(gu[4]) do
 print(n[0]=gu[4][i1][1], n[1]=gu[4][i1][2]):
od:

print(`for each `, n[j], `j>=1 , with j odd, r(n):=1 `):
print(`for each `, n[j], `j>=0 , with j even, r(n):=-1 `):

fi:

print(`Proof: it follows from the beautiful proof of Dalia Krieger`):
print(`and Aveizri Fraenkel, in their paper: `):
print(`The structure of complementary sets of integers: A 3-shift theorem`):
print(`Inter. J. of Pure and Applied Math. 10 (2004) 1-49, available on-line`):
print(`http://www.wisdom.weizmann.ac.il/~fraenkel/Papers/tsopaper1.pdf  .`):

print(`To sum up the "noisy" beginning has`, nops(gu[1])-n0+1, `terms `):
print(`The regular-shift is`, gu[2]):
print(`The irregular shifts belong to`, nops(gu[3])+nops(gu[4]), `recursive sequences. `):
print(``):
print(`Just to demonstrate how fast this is, the millionth-term,`):
t0:=time():
print(`A(1000000) equals `):
print(AnFK(n0,X,c,1000000)):
print(`and this took`, time()-t0, `seconds (once we computed the F-K scheme)`):
print(`Note that we didn't have to compute all the previous`):
print(` almost million terms! `):
print(``):
print(`Thanks for reading this computer-generated paper!`):
print(`generated by a humanly-generated Maple program`):
print(` (written by Doron Zeilberger) based on`):
print(`a humanly-generated brilliant algorithm invented (or discovered(?))`):
print(`by Dalia Krieger and Aviezri Fraenkel. `):
end: