#   LibSequences:  MAPLE UTILITY FUNCTIONS FOR SEQUENCES 
#
#   F. Vivaldi  Mar 2020
#
#  ===================================================== CONTENTS
#
# -------------- LIBRARIES
#
with(NumberTheory):
#
# -------------- CONVERSIONS
# cf2x(L)         Number in ]0,1] having L=[a1,...,an] as continued fraction
# x2cf(x,n)       n terms of the continued fraction of x
# x2cf(x)         all terms of the continued fraction of a rational x
# x2cvg(x,n)      n convergents of x
# x2cvg(x)        all convergents of a rational x
# Int2Dig(n,b)    List of digits of non-negative integer n to the base b
# Int2Dig(n,b,nd) Ditto, the first nd digits
# Dig2Int(L,b)    Integer having L as list of digits to the base b
# Int2Fibo(n)     Fibonacci expansion of a positive integer
#
# --------------RATIONAL NUMBERS
# Farey(n)     	Farey fractions in closed unit interval, and denom <=n
# Farey(-n)    	Ditto, with denom =n
# Farey(n,a,b) 	Ditto, in closed interval [a,b]
# H(x)		Height of x (x::{integer,fraction,list{integer,fraction}}
# Mdt(r,s)      mediant of r and s

# --------------SEQUENCES
# Affine(L,a,b) Affine cyclic map x->a*x+b of the elements of L
# LShift(L)     left shift on a list
# RShift(L)     right shift on a list
# isMinimal(L)  has the periodic list L minimal period?
# Cesaro(L)     Cesaro sum of a sequence
# DeCo(L)	Delay co-ordinates: the sequence of pairs of consecutive 
#  		terms of the sequence L.
# Distrib(FL)   Cumulative distribution function of a frequency list
# Freq(L)       Frequency list of the elements of a sequence L.
# FreqMerge(L)  A frequency list from a list of frequency items,
#               possibly unordered/duplicated.
# Hamming(L1,L2) Hamming distance of two sequences.
# SumSeq(L)     Partial sums of the elements of a list
#               L=[x1,x2,...] or L=[[k1,xk1],[k2,xk2],...]
# PWA(L)	Represents the piecewise affine sequence L as a list of turning points
#               [...[k,L[k]],...] 
#               If L is not PWA, then the output is the same as PlotSeq(L).
#
# --------------DYNAMICAL SYSTEMS
# Compose(F,G)  G(F(x)) (F & G are lists)
# Orbit(f,z,length)  	Orbit of map f through z, of given length; 
#  			if the orbit is periodic, it returns one period
# IOrbit(f,ic)  Orbit of invertible map f (slower than Orbit)
# IOrbits(f,L)  All periodic orbits of invertible map f over set L
# Wmap(F,W,x)   The image of x under the maps F(x,w) driven by the word W
# Worbit(F,W,x) The orbit of x under the maps F(x,w) driven by the word W
# FunOrbit(Fs,W) Functional orbit of maps Fs, driven by word W (Fs are lists)

# --------------SYMBOL SEQUENCES
# Blocks(w,W)   Decomposition of W into w-blocks (greedy algorithm)
# Blocks(w,W,any)  ditto, safe algorithm
# Complexity(W)	Complexity of word W: [C(1),C(2),...], where C(n) is 
# 		the # of substrings of length n. It stops when C(n+1)=C(n).
# Complexity(W,nmax)  ditto, with maximum length of substring.
# Index(i,L)	Index (position) of first occurrence of i in sequence L
# Indices(i,L)    List of indices of all occurrences of i in sequence L.
# IndexW(w,L)	Index (position) of first occurrence of sequence w in sequence L
# IndicesW(w,L)	  List of indices of all occurrences of seq w in seq L
# IndicesW(w,L,anything)	  ditto, but without overlapping ws
# Language(W,n)	Set of distinct words of length n in word W.
# MatrixWord(W)	Left matrix product for word W=[w_1,...,w_n]:
# 		M[w_n]M[w_{n-1}]...M[w_1] where M[x]=[[x,-1],[1,0]]
# ReCode(w,W,x) Replace all (disjoint) occurrences of w in W with x
#
# --------------STURMIAN SEQUENCES
# Sturmian(theta,gamma,n) Sturmian n-word with parameters theta,gamma
# Sturmian(p/q,gamma)     Sturmian q-word for parameter p/q,gamma
# Sturmian(p/q)           Sturmian q-word for parameter p/q,p/q
# SturmianR(theta,k)      Recursive Sturmian word (k-step recursion)
#                         for special word (gamma=theta)
# AllSturmian(n,r)        All Sturmian words of lenght n at the rational r
# AllSturmian(n,a,b)      Ditto, in interval [a,b]
#               
# --------------PLOT FACILITIES
# PlotSeq(L)    Convert a sequence L into the plot structure [...[k,L[k]]...]
# PlotSeq(L,dx) ditto, with scaled abscissa [...[dx*k,L[k]]...]
# PSPS(L)       short-hand: plot(Steps(PlotSeq(L)))
# PSPS(L,options)   ditto, with plot options 
# Steps(L)      Step-function, from a plotting sequence
# Thin(L,n)     thin a sequence: L->[L[n],L[2*n],L[3*n],...]

#===============================================================================

#  Affine_________________________________________________________________
#  Affine(L,a,b)         affine cyclic map x->a*x+b of the elements of L
#  FUNCTIONS: none
Affine:=proc(L,a,b)
local i, n:
 n:=nops(L):
 [seq(L[(i*a+b-1 mod n)+1],i=1..n)]:
 return %:
end:

#  Blocks_________________________________________________________________
#  Blocks(w,W)   decomposition of W into w-blocks: greedy algorithm
#  Blocks(w,W,any)  ditto, safe algorithm

Blocks:=proc(w,W)
global Index,IndicesW:
local greedy,B,b,i,ib,head,mayb,n,
      Iw,allIw,nIw:

   greedy:=evalb(nargs=2):
   allIw:=IndicesW(w,W):
   n,nIw:=nops(w),nops(allIw):
   if nIw=0 then return [nops(W)] fi:
   if greedy then
      Iw:=allIw:
   else
      select(i->allIw[i+1]-allIw[i]>=n,[$1..nIw-1]):
      Iw:=[seq(allIw[i],i=%),allIw[-1]]:
      nIw:=nops(Iw):
   fi:
   ib,b,mayb:=1,Iw[1],Iw[1]:
   if b=1 then head:=NULL else head:=b-1 fi:
   B:=head,0:

   while mayb<Iw[-1] do
     if Iw[ib+1]=b+n then
        B:=B,0:
        ib:=ib+1:
        b:=b+n:
        mayb:=b:
     else
        if member(b+n,Iw) then
          B:=B,0:
          b:=b+n:
          ib:=Index(b,Iw):
          mayb:=b:
        else
          ib:=ib+1:
          mayb:=Iw[ib]:
          while ib<nIw and mayb-n<b do
            ib:=ib+1:
            mayb:=Iw[ib]:
          od:
          if mayb-n>b then
             B:=B,mayb-n-b,0:
             b:=mayb:
          fi:
        fi:
     fi:
   od:

   nops(W)-b-n+1:
   if %=0 then
      return [B]
   else
      return [B,%]
   fi:
end:

#  Cesaro_________________________________________________________________
#  Cesaro(L)             Cesaro sum of a sequence
Cesaro:=proc(L::list)
global SumSeq:
local  i, isList, N:
  N:=nops(L):
  if N=0 then
    return L
  else
    isList:=type(L[1],list):
  fi:
  SumSeq(L):
  if isList then
    return [seq([%[i][1],%[i][2]/i],i=1..N)]
  else
    return [seq(%[i]/i,i=1..N)]
  fi:
end:

#  cf2x___________________________________________________________________
#  cf2x(L)    number in ]0,1] having L=[a1,...,an] as continued fraction
#  x2cf(x,n)  n terms of the continued fraction expansion of x
#  x2cf(x)    all terms of continued fraction expansion of a rational x
#  x2cvg(x,n) n convergents of x
#  x2cvg(x)   all convergents of a rational x
cf2x:=proc(L::list)
local i,n:
  n:=nops(L):
  1/L[n]:
  for i to n-1 do
     %+L[n-i]:
     1/%:
  od:
  return simplify(%)
end:

x2cf:=proc(x)
  ContinuedFraction(x):
  if nargs=2 then
     return Term(%,1..args[2])
  else
     return Term(%,all)
  fi:
end:

x2cvg:=proc(x)
  ContinuedFraction(x):
  if nargs=2 then
     return Convergent(%,1..args[2])
  else
     return Convergent(%,all)
  fi:
end:


#  Complexity_____________________________________________________________
#  Complexity(L)          Complexity of a sequence
#  Complexity(L,nmax)  
Complexity:=proc(L)
local C, j, n, k, N, Nmax, isIncreasing:
   N:=nops(L):
   if nargs=1 then Nmax:=N else Nmax:=args[2] fi:
   C:=array(1..Nmax):
   C[1]:=nops({op(L)}):
   isIncreasing:=true:
   for n from 2 to Nmax while isIncreasing do 
      {seq(L[k..(k+n-1)],k=1..N-n+1)}:
      C[n]:=nops(%):
      if C[n]=C[n-1] then isIncreasing:=false fi:
   od:
   return [seq(C[j],j=1..min(Nmax,n-1))]
end:

#  Compose________________________________________________________________
#  Compose(F,G)         G(F(x))
Compose:=proc(F,G)
local i:
  return [seq(G[F[i]],i=1..nops(F))]
end:


#  DeCo___________________________________________________________________
#  DeCo(L)             Delay co-ordinates of a sequence
DeCo:=proc(L::list)
local i,N:
  N:=nops(L):
  if N < 2 then 
     return [] 
  else
     return [seq([L[i],L[i+1]],i=1..N-1)]
 fi
end:

#  Distrib________________________________________________________________
# >Distrib(L);
# >Distrib(L,flag);
#  OUTPUT: Cumulative distribution function of a discrete measure,
#             given as a frequency list L=[...,[x,#x],...]
#             flag=1: [default] the measure of [x,#x] is (proportional to) #x
#             flag=2: the measure of [x,#x] is (proportional to) x*#x
# FUNCTIONS: none

Distrib:=proc(L::list)
local  i, flag, n, N, A, sortL, S:

  if nargs=1 then
     flag:=1:
  else
     flag:=args[2]
  fi:
  if flag=1 then
     N:=add(x[2],x=L):
  else
     N:=add(x[1]*x[2],x=L):
  fi:
  n:=nops(L):
  sortL:=sort(L,(x,y)->x[1]<y[1]):
  A:=array(1..n):
  S:=0:
  for i to n do
    sortL[i]:
    if flag=1 then
       S:=S+%[2]:
    else
       S:=S+%[1]*%[2]:
    fi:
    A[i]:=S/N:
  od:
  [[sortL[1][1],0],seq([sortL[i][1],A[i]],i=1..n)]:

end:

#  Farey__________________________________________________________________
#  Farey(n)   Farey fractions in closed unit interval with denominator =< n
#  Farey(-n)  Farey fractions in closed unit interval with denominator = n
#  Farey(n,a,b)   Ditto, in closed interval [a,b]
Farey:=proc(n::{posint,negint})
local Coprime, x,y,a,b:
 Coprime:=(x,y)->if igcd(x,y)=1 then x/y else NULL fi:
 if nargs=3 then
   a,b:=args[2],args[3]
 else
   a,b:=0,1:
 fi:
 if abs(n)=1 then
    [a,b]
 elif n>0 then
    [seq(seq(Coprime(x,y),x=ceil(y*a)..floor(y*b)),y=1..n)]:
 else
    [seq(Coprime(x,-n),x=ceil(-n*a)..floor(-n*b))]:
 fi:
 return sort(%)
end:

#  Freq___________________________________________________________________
#    frequency of occurrence of the elements of a list
#   >Freq(L);
#    INPUT:  L:list
#    OUTPUT: a list of lists
#    FUNCTIONS: none
#    SYNOPSIS:
#      Compute the frequency of occurrence of the elements of a list.
#      The number of operands is equal to the number of distinct
#      elements of the input, which are sorted.
#      Each operand consists of the given element, with its frequency.
#
Freq:=proc(L::list)
local  sortOpsL, x:
  sortOpsL:=sort([op({op(L)})]):
  [seq([x, nops(select((l,x)->evalb(l=x), L, x))], x=sortOpsL)]
end:

#  FreqMerge______________________________________________________________
#    A frequency list from a list of frequency items
#   >FreqMerge(L);
#    INPUT:  L:list of frequency items [x,n]
#    OUTPUT: a frequency list (see Freq)
#    FUNCTIONS: none
#    SYNOPSIS:
#      Turn a pseudo-frequency list with (possibly) unordered/duplicated 
#      terms, into a frequency list (n need not be an integer).
#
FreqMerge:=proc(L::list)
local i,k,n,nL,sortL,mf,xnow,fnow:
   nL:=nops(L):
   n:=nops({op(map(z->z[1],L))}):mf:=array(1..n):
   sortL:=sort(L,(x,y)->x[1]<y[1]):
   i:=1:fnow:=0:
   for k to nL do
     xnow:=sortL[k][1]:
     fnow:=fnow+sortL[k][2]:
     if k=nL or sortL[k+1][1]<>xnow then
       mf[i]:=[xnow,fnow]:
       fnow:=0:
       i:=i+1:
     fi:
   od:
   return [seq(mf[j],j=1..n)]
end:

#  FunOrbit_______________________________________________________________
#  FunOrbit(Fs,W) functional orbits of maps (lists) Fs driven by W
FunOrbit:=proc(Fs,W)
global Compose:
local n,i,fff,f:
  f[0],f[1]:=Fs[1],Fs[2]:
  n:=nops(W):
  fff:=array(1..n):
  fff[1]:=f[W[1]]:
  for i from 2 to n do fff[i]:=Compose(fff[i-1],f[W[i]]) od:
  return [seq(fff[i],i=1..n)]:
end:


#  H______________________________________________________________________
#  H(x)		Height of x (x::{integer,fraction,list{integer,fraction}}
# 
H:=proc(x)
 if type(x,integer) then
   return max(abs(x),1)
 else
   return max(map(H,[op(x)]))
 fi:
end:

# Index__________________________________________________________________
# Index(l,L)            index of first occurence of term l of sequence L
# Indices(l,L)          indices of all occurrences
Index:=proc(l,L)
local i,n:
 n:=nops(L):
 if not member(l,L) then
   return NULL
 else
   for i while l<>L[i] to n do od:
   return i
 fi:
end:

Indices:=proc(j,L)
global Index:
local i,n,inow,Is:
  n:=nops(L):
  i:=1:
  Is:=NULL:
  while i<n do 
     inow:=Index(j,L[i..n]):
     if inow=NULL then
        return [Is]
     else
        i+inow:
        Is:=Is,%-1:
        i:=%%:
     fi:
  od:
  return [Is]:
end:

# IndexW_________________________________________________________________
# IndexW(w,W)          index of first occurence of word w in word W
# IndicesW(w,W)        index of all occurences of word w in word W
# IndicesW(w,W,any)      ditto, without overlapping ws
#
IndexW:=proc(w,W)
global Index:
local i,nw,nW:
 nw,nW:=nops(w),nops(W):
 if nw>nW then 
    return NULL 
 elif nw=1 then
    return Index(w[1],W)
 else
    for i to nW-nw+1 do
       if W[i..i+nw-1]=w then
          return i
       fi:
    od:
    return NULL
 fi:
end:

IndicesW:=proc(w,W)
global IndexW:
local i,inow,nw,nW,idxrel,Is,noverlap:
  noverlap:=evalb(nargs>2):
  nw,nW:=nops(w),nops(W):
  i:=1:
  Is:=NULL:
  while i<=nW-nw+1 do 
     idxrel:=IndexW(w,W[i..nW]):
     if idxrel=NULL then
        return [Is]
     else
        inow:=i+idxrel-1:
        Is:=Is,inow:
        if noverlap then
           i:=inow+nw:
        else
           i:=inow+1:
        fi:
     fi:
  od:
  return [Is]:
end:

# Hamming________________________________________________________________
#   Hamming distance of two sequences of equal length
#  >Hamming(L1,L2)
Hamming:=proc(L1,L2)
local i, h, n:
  n:=nops(L1):
  if nops(L2)<>n then ERROR("Sequences have unequal length") fi:
  h:=0:
  for i to n do
    if L1[i]<>L2[i] then h:=h+1 fi
  od:
  return h
end:

# Int2Dig,Dig2Int,Int2Fibo_______________________________________________
#   convert from b-adic digits to integer, and vice-versa
#  >Int2Dig(i,b);
#  >Int2Dig(i,b,nd);
#  >Dig2Int(L,b);
#  >Int2Fibo(n);

Int2Dig:=proc(x::nonnegint,b::posint,nd::nonnegint)
local n, i:
 if x=0 then
    if nargs = 2 then return [0] else return([0$nd]) fi:
 fi:
 if nargs = 2 then
    n:=max(0,floor(evalf(log[b](abs(x)))) + 1)
 else
    n:=nd
 fi:
 [seq(irem(iquo(x,b^(i-1)),b),i=1..n)]
end:

Dig2Int:=proc(L::list(nonnegint),b::posint)
local i:
 add(op(i,L)*b^(i-1),i=1..nops(L)):
end:

Int2Fibo:=proc(n::posint)
local i, idig, F, FD, x, nFibo:
 F:=[1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,6765,10946,
     17711,28657,46368,75025,121393,196418,317811,514229,832040,1346269,
     2178309,3524578,5702887,9227465,14930352,24157817,39088169,63245986,102334155,165580141]:
 if n>F[-2]+F[-1] then ERROR("not enough Fibonacci numbers") fi:
 x:=n:
 for i while F[i]<=x do od:
 nFibo:=i-1:FD:=array(1..nFibo,[0$nFibo]):
 idig:=nFibo:#FD[idig]:=1:
 while x>0 do
   FD[idig]:=1:
   x:=x-F[idig]:
   for i while F[i]<x do od:
   idig:=i:
 od:
 return [seq(FD[i],i=1..nFibo)]
end:

# Language_______________________________________________________________
# Language(L,n)       set of all distinct n-words in sequence L
Language:=proc(L,n)
local N, k:
  N:=nops(L):
  {seq(L[k..(k+n-1)],k=1..N-n+1)}:
  return %
end:

# isMinimal______________________________________________________________
#          check whether of not a periodic list L has minimal period
isMinimal:=proc(L::list)
local n, t, Lnow:
  n:=nops(L):
  Lnow:=[op(L[2..n]),L[1]]:
  for t while L <> Lnow do
    Lnow:=[op(Lnow[2..n]),Lnow[1]]
  od:
  evalb(t=n)
end:

#  MatrixWord__________________________________________________________
#  MatrixWord(W) left product of matrices M[x]=[[x,-1],[1,0]] associated with word W
MatrixWord:=proc(W)
local A,n:
  A:=Matrix([[x,-1],[1,0]]):
  Matrix([[1,0],[0,1]]):
  for n to nops(W) do
     Multiply(subs(x=W[n],A),%)
  od:
  return map(factor,evalm(%))
end:

#  Mdt_________________________________________________________________
#  Mdt(r,s)  mediant of r and r
Mdt:=proc(r::{integer,fraction},s::{integer,fraction})
 return (num(r)+num(s))/(denom(r)+denom(s))
end:

#  Orbit_______________________________________________________________
#       orbit of F through ic, with given length
Orbit:=proc(f::procedure,z0,length::posint)
local  nit, t, z, Z:
  Z:=array(1..length):
  z:=z0:
  Z[1]:=z:
  z:=f(z):
  for t from 2 while z<>z0 to length do
     Z[t]:=z:
     z:=f(z):
  od:
  nit:=min(length,t-1):
  return([seq(Z[t],t=1..nit)])
end:

# IOrbit___________________________________________________________________
#   orbit of an invertible map f with initial condition ic
#  >IOrbit(f,ic);
#   OUTPUT: a list [cycle]

IOrbit:=proc(f::procedure,ic)
local T, X, t:

 t:=1:
 T:=table([t=ic]):
 X:=f(ic):
 while X <> ic do
   t:=t+1:
   T[t]:=X:
   X:=f(X)
 od:
 return([seq(T[i],i=1..t)]):
end:

# IOrbits_______________________________________________________________
#   cycle decomposition of invertible map f over L
#  >IOrbits(f,L);
#   OUTPUT: a list of lists [[cycle_1],...,[cycle_n]]

IOrbits:=proc(f::procedure,L::list)
local ic, norbits, L_now, orbit_now, T:
global IOrbit:

 T:=table():
 norbits:=0:
 L_now:={op(L)}:
 while nops(L_now) > 0 do
   ic:=op(1,L_now):
   orbit_now:=IOrbit(f,ic):
   L_now:=L_now minus {op(orbit_now)}:
   norbits:=norbits+1:
   T[norbits]:=orbit_now:
 od:
 map(x->op(x),[entries(T)]):

end:

# PlotSeq________________________________________________________________
#    the plot structure of a list
PlotSeq:=proc(L::list)
local dx:
  if nargs=1 then dx:=1 else dx:=args[2] fi:
  return([seq([dx*k,L[k]],k=1..nops(L))])
end:

# PSPS___________________________________________________________________
# PSPS(L)   short-hand for plot(Steps(PlotSeq(L)))
PSPS:=proc(L)
global PlotSeq,Steps:
local opzioni:
  if nargs=1 then
    opzioni:=NULL
  else
    opzioni:=args[2]
  fi:
  return plot(Steps(PlotSeq(L)),opzioni)
end:

# PWA___________________________________________________________________
# PWA(L) represent a piecewise affine sequence as a sequence of turning points
PWA:=proc(L::list)
local n, i, old, new, V:
  n:=nops(L):
  if n=1 then
     ERROR("input list must have at least two elements"):
  fi:
  old:=infinity:V:=NULL:
  for i to n-1 do
     new:=L[i+1]-L[i]:
     if new <> old then
        V:=V,[i,L[i]]:
        old:=new
     fi:
  od:
  V:=V,[n,L[n]]:
  return [V]
end: 

#  ReCode_________________________________________________________________
#  ReCode(w,W,x)  replace every disjoint occurrence of w in W with x
ReCode:=proc(w,W,x)
global IndicesW:
local i,k,n,nws,N,Nout,Iw,L,inside,outside:
   n,N:=nops(w),nops(W):
   Iw:=IndicesW(w,W,1):
   nws:=nops(%):
   Nout:=N+nws*(1-n):
   inside:=map(z->(seq(k,k=z..z+n-1)),Iw):
   outside:=select(z->not member(z,inside),[$1..N]):
   L:=array(1..Nout):
   i,k:=0,1:
   while k <= N do
     i:=i+1:
     if member(k,outside) then
        L[i]:=W[k]:
        k:=k+1:
     else
        L[i]:=x:
        k:=k+n:
     fi:
   od:
   return [seq(L[i],i=1..Nout)]:
end:

#  Shift________________________________________________________________
#  LShift(L) - left shift of a list                         
#  RShift(L) - right shift of a list                         
#
LShift := proc(L::list)
local i, n:                                                                      
  n := nops(L):
  [seq(L[i],i=2..n), L[1]]
end:

RShift := proc(L::list)
local i, n:
  n := nops(L):
  [L[n], seq(L[i],i=1..n-1)]
end:


#  Steps__________________________________________________________________
#  Steps(L)             step-function, from a plotting sequence
#    FUNCTIONS: none
Steps:=proc(L)
local  N, S, k:
  N:=nops(L):
  S:=array(1..2*N-1):
  S[1]:=L[1]:
  for k from 1 to N-1 do
    S[2*k]:=[L[k+1][1],L[k][2]]:
    S[2*k+1]:=L[k+1]:
  od:
  return [seq(S[k],k=1..2*N-1)]:
end:

# Sturmian_______________________________________________________________
# Sturmian(theta,gamma,n)  Sturmian n-word for parameters theta, gamma
# Sturmian(p/q,gamma)      Sturmian q-word for parameter p/q,gamma
# Sturmian(p/q)            Sturmian q-word for parameter p/q,p/q
# SturmianR(theta,n)       Recursive Sturmian word (n-step recursion)
#                          for special word (parameters theta,theta)
# AllSturmian(n,r)         All Sturmian words of lenght n at the rational r
# AllSturmian(n,a,b)       Ditto, in interval [a,b]

Sturmian:=proc(theta)
local gamma,j,L,n,partition,x:
  if nargs=1 then
     gamma:=theta:
     n:=denom(theta):
  elif nargs=2 then
     gamma:=args[2]:
     n:=denom(theta):
  else 
     gamma:=args[2]:
     n:=args[3]:
  fi:
  x:=gamma:
  L:=array(1..n):
  partition:=1-theta:
  for j to n do
      if x<partition then
          L[j]:=0
      else
          L[j]:=1
      fi:
      x:=frac(x+theta)
  od:
  return [seq(L[j],j=1..n)]:
end:

SturmianR:=proc(theta,n)
local i,cf,cf2x:
  cf:=ContinuedFraction(theta):
  if n=0 then
     return [0]
  elif n=1 then
     return [seq(0,k=1..Term(cf,1)-1),1]
  fi:
  cf2x:=[Term(cf,1)-1,op(Term(cf,2..n))]:
  1:0:
  for i to n do
     seq(op([%]),k=1..cf2x[i]),%%
  od:
  return [%]
end:

AllSturmianQ:=proc(n::posint,r)
global Affine,Sturmian:
local q,k,S,Sbig:
  if n=1 then return {[0],[1]} fi:
  S:=Sturmian(r):
  q:=denom(r):
  Sbig:=[seq(op(S),k=1..ceil(n/q))]:
  {seq(Affine(Sbig,1,k)[1..n],k=0..q-1)}:
  return %
end:

AllSturmian:=proc(n::posint)
global Farey,AllSturmianQ:
local a,b,allF,low,high:
  if n=1 then 
     return {[0],[1]} 
  elif nargs=2 or args[2]=args[3] then 
     return AllSturmianQ(n,args[2]) 
  fi:
  a,b:=args[2],args[3]:
  if a>b then
     ERROR("second argument cannot be greater than third"):
  fi:
  allF:=Farey(n,0,1/2):
  low:=max(op(select(z->z<=a,allF))):
  high:=min(op(select(z->z>=b,allF))):
  select(z->z>=low and z<=high,allF):
  {seq(op(AllSturmianQ(n,r)),r=%)}:
  return %
end:

#  SumSeq_________________________________________________________________
#  SumSeq(L)             sequence of partial sums
SumSeq:=proc(L::list)
local  i, isList, N, S, LS:
  N:=nops(L):
  if N=0 then
    return L
  else
    isList:=type(L[1],list):
  fi:
  LS:=array(1..N):
  S:=0:
  for i to N do
    if isList then
       S:=S+L[i][2]:
       LS[i]:=[L[i][1],S]
    else
       S:=S+L[i]:
       LS[i]:=S
    fi
  od:
  return [seq(LS[i],i=1..N)]
end:

#  Thin___________________________________________________________________
#  Thin(L,n)              extract the subsequence L_{ni} from L_i
Thin:=proc(L::list,n::posint)
  return [seq(L[n*i],i=1..iquo(nops(L),n))]
end:

# Wmap, Worbit__________________________________________________________
# Wmap:=proc(F,W,ic)       The image of ic under the maps F driven by the word W
# Worbit:=proc(F,W,ic)     The orbit of ic under the maps F driven by the word W
Wmap:=proc(F::procedure,W::list,ic)
local i, x:
  x:=ic:
  for i to nops(W) do
     x:=F(x,W[i])
  od:
  return x
end:

Worbit:=proc(F::procedure,W::list,ic)
local i, nW, X:
  nW:=nops(W):
  X:=array(0..nW):
  X[0]:=ic:
  for i to nW do
     X[i]:=F(X[i-1],W[i])
  od:
  return [seq(X[i],i=0..nW)]
end: