#   LibSequences:  MAPLE UTILITY FUNCTIONS FOR SEQUENCES 
#
#   F. Vivaldi.  Version 1,  Apr 2016
#
#  ===================================================== CONTENTS
#
#  -------------UTILITY FUNCTIONS FOR INTEGERS
# 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
#
#  -------------UTILITY FUNCTIONS FOR 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}}

#  --------------UTILITY FUNCTIONS FOR SEQUENCES
# Affine(L,a,b) Affine cyclic map x->a*x+b of the elements of L
# 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)   Distribution function of a frequency list
# Freq(L)       Frequency list of the elements of a sequence L.
# Index(x,L)	Index (position) of first occurrence of x in sequence L
# IndexW(w,L)	Index (position) of first occurrence of sequence w in sequence L
# PlotSeq(L)    Convert a sequence L into the plot structure [...[k,L[k]]...]
# PlotSeq(L,dx) ditto, with scaled abscissa [...[dx*k,L[k]]...]
# 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).
# Steps(L)      Step-function, from a plotting sequence
# SumSeq(L)     Partial sums of the elements of a list
#               L=[x1,x2,...] or L=[[k1,xk1],[k2,xk2],...]
# PSPS(L)       short-hand: plot(Steps(PlotSeq(L)))
# PSPS(L,options)   ditto, with plot options 

#  --------------UTILITY FUNCTIONS FOR SHIFT SPACES
# Complexity(L)	Complexity of sequence L: [C(1),C(2),...], where C(k) is 
# 		the # of substrings of length k. It stops when C(n+1)=C(n).
# Complexity(L,nmax)  ditto, with maximum length of substring.
# Language(L,n)	Set of distinct words of length n in sequence L.
# Wmap(F,L,x)   The image of x under the maps F driven by the word L
# Worbit(F,L,x) The orbit of x under the maps F driven by the word L

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

#  Affine__________________________________________________________General
#  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:

#  Cesaro__________________________________________________________General
#  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:

#  Complexity______________________________________________________General
#  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:

#  DeCo____________________________________________________________General
#  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_________________________________________________________General
# >Distrib(L);
# >Distrib(L,flag);
#  OUTPUT:  Distribution function of a frequency list L=[...,[x,#x],...]
#             flag=1  [default]: the measure of [x,#x] is proportional to x*#x
#             flag=2: the measure of [x,#x] is proportional to #x
# FUNCTIONS: none

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

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

end:

#  Farey___________________________________________________________General
#  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(a)..floor(-n*b))]:
 fi:
 return sort(%)
end:

#  Freq____________________________________________________________General
#    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:


#  H_______________________________________________________________General
#  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___________________________________________________________General
# Index(l,L)            index of first occurence of term l of sequence L
Index:=proc(l,L)
local i,d:
 d:=nops(L):
 if not member(l,L) then
   return NULL
 else
   for i while l<>L[i] to d do od:
   return i
 fi:
end:

# IndexW__________________________________________________________General
# IndexW(w,W)          index of first occurence of word w in word W
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:

# Int2Dig,Dig2Int,Int2Fibo________________________________________General
#   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________________________________________________________General
# 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______________________________________________________General
# isMinimal(F,W,i) Returns [i,j] where j=F(i) if the orbit of F driven by W
#            starting at i is minimal (there are no i,j between the first and the
#            last points. Otherwise it returns [].
isMinimal:=proc(F::procedure,W::list,i)
local j,k,nW,Orb,X:
  nW:=nops(W):
  X:=array(0..nW):
  X[0]:=i:
  for k to nW do
     X[k]:=F(X[k-1],W[k]):
     if X[k]=i and k<nW then
        return []
     fi:
  od:
  Orb:=[seq(X[k],k=0..nW)]:
  j:=Orb[-1]:
  if nops(Orb)=2 or not member(j,Orb[2..-2]) then
    return [i,j]
  else
    return []
  fi:
end:

# Wmap, Worbit___________________________________________________General
# 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:

# PlotSeq_________________________________________________________General
#    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:

# PWA____________________________________________________________General
# 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: 

#  Steps___________________________________________________________General
#  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:

#  SumSeq__________________________________________________________General
#  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:

# PSPS____________________________________________________________General
# 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: