######################################################################
##SlaterValez.txt: Save this file as  SlaterValez.txt                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read SlaterValez.txt<Enter>                        #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################

read `DataSlaterValez.txt`: 
#Created:  May 31, 2019
 
print(`Created:  May 31, 2019`):
print(` This is SlaterValez.txt `):
print(`It is one of the packages that accompany the article `):
print(` Guessing the Elusive Patterns in the Slater-Valez sequence (aka OEIS A081145) `):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):
print(`dedicated to Aviezri S. Fraenkel (b. June 7, 1929) on his 90th birthday`):

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/  .`):

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 plotting procedures type ezraP();, 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):

ezraP:=proc()

if args=NULL then
 print(` The plotting procedures are: Plota, Plotd  `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: AlDel,AlDelF,  mex , SSeqa, SSeqd `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: a , d, Estc, Mila, Seqa, SeqaPC, Seqd, SeqdPC,  StSt `):
 print(` `):



elif nops([args])=1 and op(1,[args])=a then
print(`a(n): the n-th term of the Slater-Valez sequence, A081145. Try:`):
print(`a(10);`):

elif nops([args])=1 and op(1,[args])=AlDel then
print(`AlDel(c): the equation satisfied by heuristic slope alpha if the prob. of going from state 2 to state 1 is c. Try:`):
print(` AlDel(c); `):

elif nops([args])=1 and op(1,[args])=AlDelF then
print(`AlDelF(c): the heuristic slopes of the three lines, if the prob. of going from state 2 to state 1 is c. `):
print(`It returns the three slopes followed by their frequency. Try:`):
print(`AlDelF(0.075); `):

elif nops([args])=1 and op(1,[args])=d then
print(`d(n): abs(a(n)-a(n-1)). Try:`):
print(`d(10);`):

elif nops([args])=1 and op(1,[args])=Estc then
print(`Estc(N1,N2): estimates the probability of the transition 2->1 in the Slater-Valez sequence. Try:`):
print(`Estc(1000,2000):`):

elif nops([args])=1 and op(1,[args])=mex then
print(`mex(S): inputs a set of positive integers, outputs the smallest  pos. integer that is NOT in S. Try:`):
print(`mex({1,2,3,7});`):

elif nops([args])=1 and op(1,[args])=Mila then
print(`Mila(N): the word, of length N,  in the alphabet {1,2,3} describing the states of OEIS sequence A081145, the Slater-Valez sequence`):
print(`an integer is in state 1 if a(n)/n<=1, in state 2 if 1<a(n)/n<=2 and in state 3 if a(n)/n>2 . Try:`):
print(`Mila(1000);`):

elif nops([args])=1 and op(1,[args])=Plota then
print(`Plota(N): plots all the points [n,a(n)] for n<=N. Try:`):
print(`Plota(100);`):

elif nops([args])=1 and op(1,[args])=Plotd then
print(`Plotd(N): plots all the points [n,d(n)] for n<=N. Try:`):
print(`Plotd(100);`):

elif nops([args])=1 and op(1,[args])=Seqa then
print(`Seqa(N): The list consisting of the first N terms of the Slater-Valez sequence, A081145. Try:`):
print(`Seqa(100);`):

elif nops([args])=1 and op(1,[args])=SeqaPC then
print(`SeqaPC(N): a pre-computed, hence much faster version of Seqa(N) (q.v.) for N<=150000. Requires the data file`):
print(` DataSlaterValez.txt , avaialble from `):
print(``):
 print(`http://sites.math.rutgers.edu/~zeilberg/tokhniot/DataSlaterValez.txt  .`):
print(``):
print(`Try:  `):
print(`SeqaPC(10000);`):

elif nops([args])=1 and op(1,[args])=Seqd then
print(`Seqd(N): The list consisting of the first N terms d(n):=abs(a(n)-a(n-1)), starting at n=2. Try:`):
print(`Seqd(100);`):

elif nops([args])=1 and op(1,[args])=SeqdPC then
print(`SeqdPC(N): Like Seqd(N), but faster using pre-computed values. Try:`):
print(`Requires the same data set as for SeqaPC(N) (q.v.). Try: `):
print(`SeqdPC(100);`):

elif nops([args])=1 and op(1,[args])=SSeqa then
print(`SSeqa(N): The three subsequences of Seqa(N), according to their type. Try:`):
print(`SSeqa(300);`):

elif nops([args])=1 and op(1,[args])=SSeqd then
print(`SSeqd(N): The two subsequences according to their type. Try:`):
print(`SSeqd(300);`):

elif nops([args])=1 and op(1,[args])=StSt then
print(`StSt(c): The steady-state distribution of the states if 1->2->3->1 but prob. of 2->1 is c.`):
print(` Try: `):
print(`StSt(c); `):
print(`StSt(0.07516); `):

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






Digits:=30:

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

#a(n): the n-th term of the Slater-Valez sequence, A081145. Try:
#a(10);
a:=proc(n) local x,i1:
option remember:
if n=1 then
 RETURN(1):
fi:

mex({seq(a(i1),i1=1..n-1),seq(a(n-1)+d(i1),i1=2..n-1),seq(a(n-1)-d(i1),i1=2..n-1)  } ):
end:

d:=proc(n)
option remember:
if n=1 then
 RETURN(1):
fi:
abs(a(n)-a(n-1)):
end:

Seqa:=proc(N) local n:
[seq(a(n),n=1..N)]:
end:

Seqd:=proc(N) local n:
[seq(d(n),n=2..N)]:
end:


#SeqdPC(N): The first N terms of the sequence d(n). It uses pre-computer values. Try:
#SeqdPC(1000);
SeqdPC:=proc(N) local n,gu:
gu:=SeqaPC(N):
[seq(abs(gu[n]-gu[n-1]),n=2..N)]:
end:


 Plota:=proc(N) local gu,n:

gu:=SeqaPC(N):

gu:={seq([n,a(n)],n=1..N)}:

plot(gu,style=point):

end:

 Plotd:=proc(N) local gu,n:
gu:=Seqa(N):
gu:={seq([n,abs(gu[n]-gu[n-1])],n=2..N)}:

plot(gu,style=point):

end:


#Mila(N): the word in the alphabet {1,2,3} describing the states of OEIS sequence A081145, the Slater-Valez sequence
Mila:=proc(N)  local gu,lu,n:
gu:=SeqaPC(N):
lu:=[]:


for n from 1 to N do
  if gu[n]<=n then
  lu:=[op(lu),1]:
  elif gu[n]<=2*n then
  lu:=[op(lu),2]:
   else
  lu:=[op(lu),3]:
 fi:
od:
lu:
end:


#SSeqa(N): The three subsequences according to their type. Try:
#SSeqa(300);
SSeqa:=proc(N) local n,gu,gu1,gu2,gu3:

gu:=SeqaPC(N):

gu1:=[]:gu2:=[]:gu3:=[]:

for n from 1 to N do
  if gu[n]<=n then
    gu1:=[op(gu1),[n,gu[n]] ]:
  elif gu[n]<=2*n then
    gu2:=[op(gu2),[n,gu[n]] ]:
   else
    gu3:=[op(gu3),[n, gu[n]] ]:
 fi:
od:

[gu1,gu2,gu3]:
end:




#SSeqd(N): The two subsequences according to their type. Try:
#SSeqd(300);
SSeqd:=proc(N) local n,gu,gu1,gu2,lu:

gu:=SeqaPC(N):

gu1:=[]:gu2:=[]:

for n from 2 to nops(gu) do
 lu:=abs(gu[n]-gu[n-1]):
  if  lu<=n then
    gu1:=[op(gu1),[n,lu]]:
  else
    gu2:=[op(gu2),[n,lu]]:
 fi:
od:

[gu1,gu2]:
end:




#StSt(c): The steady-state distribution of the states if 1->2->3->1 but prob. of 2->1 is c.
#Try:
#StSt(0.07516);
StSt:=proc(c) local x,eq,var:
var:={x[1],x[2],x[3]}:
eq:={x[1]=c*x[2]+x[3],x[2]=x[1],x[3]=(1-c)*x[2],x[1]+x[2]+x[3]=1}:
var:=solve(eq,var):
subs(var,[x[1],x[2],x[3]]):
end:


#AlDel(c): the equation satisfied by heuristic slope alpha if the prob. of going from state 2 to state 1 is c. Try: AlDel(c)
AlDel:=proc(c) local del,al:

del:=normal(2/(3-c)+(1-c)/(3-c)/2):

numer(1/(3-c)*1/(al)+ 1/(3-c)*1/(al+del)+(1-c)/(3-c)*1/(al+2*del)-1):


end:


#AlDelF(c): the heuristic slopes if the prob. of going from state 2 to state 1 is c. 
#It returns the three slopes followed by their frequency. Try:
#AlDelF(0.0.75)`
AlDelF:=proc(c) local del,al:

del:=normal(2/(3-c)+(1-c)/(3-c)/2):



al:=max(fsolve(numer(1/(3-c)*1/(al)+ 1/(3-c)*1/(al+del)+(1-c)/(3-c)*1/(al+2*del)-1),al)):

[[al,al+del,al+2*del],[1/(3-c),1/(3-c),(1-c)/(3-c)]]:


end:


#Estc(N1,N2): estimates the probability of the transition 2->1 in the Slater-Valez sequence. Try:
#Est(1000,2000):
Estc:=proc(N1,N2) local gu,lu,lu1,i:
gu:=Mila(N2):

gu:=[op(N1..N2,gu)]:

lu:=0:
lu1:=0:

for i from 1 to  nops(gu)-1 do
 if gu[i]=2 then
   lu:=lu+1:
   if gu[i+1]=1 then
     lu1:=lu1+1:
   fi:
 fi:
od:

evalf(lu1/lu):

end: