#   LibDyn:  MAPLE TOOLS FOR DISCRETE DYNAMICS
#
#   F. Vivaldi.  Version 1.2,  Feb 2016
#
# ------------------------------- Local libraries
if lavoro then
  read(`/home/network/fvivaldi/research/lib/LibSequences`):
else
  read(`/Users/fv/Documents/lib/LibSequences`):
fi:

#  ===================================================== NOTATION
#
#  L,S,LS	a list, a set, a list or set
#  z		z=[x,y] a point
#  Xs		the boundary points of a partition, left to right.
#  Cs		a list of integers in the range 1,...,n+1, where n=#Xs
#		representing an orbit code.
#  Fs		a list of n+1 functions, defining a map
#  Js		a list of n+1 entries for the corresponding jacobians
#
#  ===================================================== MAPS (at the end of this document)
#  F(z)		Piecewise affine map with 3 pieces
#                  function list: Flist=[f1,f2,f3], 
#                  parameters: aF[1,2,3], bF[1,2,3] (linear and constant coeffs)
#		               xF[1,2]     (location of boundaries of partition)
#		   if F is continuous, specify the additional parameter
#			       yF=F(xF[1])
#                  and then set:
#                  bF:=[yF-aF[1]*xF[1],yF-aF[2]*xF[1],yF-aF[2]*xF[1]+aF[2]*xF[2]-aF[3]*xF[2]]:
#  Fdsp(z)	Dissipative version. The Jacobian is the global parameter "dsp".
#  G(z)		Continuous piecewise affine map with 2 pieces
#                  function list: Glist, 
#                  parameters: aG[1,2], bG[1,2], xG
#		   For G continuous, specify the parameter yG=G(xG), and then set:
#		   bG:=[yG-aG[1]*xG[1],yG-aG[2]*xG[1]]:
#
#  ===================================================== CONTENTS
#  Code(z,Xs)	Code of the first co-ordinate of z wrt the partition boundaries Xs
#  isOK(z,Fs,Cs,Xs)	Is the periodiic point z ok?
#  Jac(Js,Cs)	Jacobian of the Jacobian list Js and code Cs
#  Tra(Js,Cs)	Sequence of traces of the Jacobian list Js and code Cs
#  Map(z,Fs,Cs) Map generated by the function list Fs and the code Cs
#  PeriodicPoint(Fs,Cs)	The periodic point for the function list Fs and periodic code Cs
#  
#  Orbit(F,z,length)  	Orbit of map F through z, of given length; 
#  			if the orbit is periodic, it returns one period
#  OrbitData(Fdata,z,length)  
#  OrbitData(Fdata,z,length,flag)  
#               list of pvalues and codes of points of the
#		orbit of map F through z, of given length; 
#               The flag (default=0) determines the argument of pvalue
#               flag=0,1,2  argument=z,z[1],z[2]
#  		if the orbit is periodic, it returns one period.
#		The map data are Fdata=[F,P,Xs], where P=list of primes, and Xs is a partition.
#  GcdData(F,z,length)  
#  GcdData(F,z,length,delay)  
#  GcdData(F,z,length,delay,flag)  
#               various data lists about gcd(x_t,x_{t-n}).
#		F: map, z: initial cond., length: #output sequence (#iterates=length+delay), delay: n, 
#		Let x_t=a_t/b_t, 
#		flag		output
#		1		gcd(a_t,a_{t-n})  [Default]
#		11		ditto, in factored form
#		12		ditto, frequency of primes
#		13		[t,{newPs}] where t is the time of first appearance of the new primes newPs
#		2		gcd(a_t,a_{t-n})/gcd(b_t,b_{t-n})
#		21		ditto, in factored form
#		22		ditto, frequency of primes
#		3		same as 1, but for y_t
#		31		ditto, in factored form
#		32		ditto, frequency of primes
#  AutoCorr(F,z,length,maxtau) 
#  		Autocorrelation function of 1/gcd(x_t,x_{t+tau})
#
#  -------------UTILITY FUNCTIONS FOR ARITHMETIC
#  Artin(f,p)	Decomposition of f modulo p  (the list of sorted inertia 
#  		degrees of the irreducible factors of f modulo p).
#  h(x)		Logarithmic height of x (floating-point)
#  PoBoH(n)	Points in Q^2 of height bounded by n (height=|n| if n is negative)
#  pvalue(r,p)	p-adic value of rational number r (p=infinity is allowed).
#  punit(r,p)	unit part of a rational r

#  --------------UTILITY FUNCTIONS FOR SEQUENCES
#  ZeroGaps(L)  Sequence of gaps between consecutive occurrences
#  		of the value 0 in the sequence L.
#  Zsig(L)	Zsigmondy set of an integer sequence

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

# Artin___________________________________________________________NumThy
# Artin - Artin map: cycle decomposition of f modulo p into irreducbles
# >Artin(f,p);
#   INPUT:  f:polynomial - a polynomial over Z
#           p:prime - the characteristic 
#   OUTPUT: A list of lists containing sorted inertia degrees of the
#           irreducible factors of f(x) modulo p.
#   FUNCTIONS: none
#
 alias(artin=Artin):

 Artin:=proc(f::polynom(integer),p::prime)
 local x, d, fmodp, g, L:

 x:=indets(f):
 if not nops(xnow)=1 then
   ERROR(f,`is not a univariate polynomial`)
 fi:
 d:=degree(f,x):
 fmodp:=Factor(f) mod p:
 if type(fmodp,`+`) or type(fmodp,symbol) then
   [d]
 elif type(fmodp,`*`) then
     L:=NULL:
     for g in op(fmodp) do
     if type(g,`+`) or type(g,symbol) then
           L:=L,degree(g,x)
     else
           L:=L,-degree(op(1,g))
     fi
   od:
   sort([L]):
 elif type(fmodp,`^`) then
   sort([-degree(op(1,fmodp),x)])
 fi
end:

#  Code____________________________________________________________General
#  Code(z)   Code of the first co-ordinate of z with respect to the partition L
Code:=proc(z::list,L::list)
local  x,n,N:
   x:=z[1]:
   N:=nops(L):
   for n to N do
     if x<L[n] then return n fi
   od:
   return N+1
end:

#  AutoCorr___________________________________________________________General
#  Autocorrelation function
#
AutoCorr:=proc(F,z0,length,maxtau)
local gcdNum, gcdDen,AC,Outtype,Z,
      tau, t, xy, i:
# -------------------------------- functions
  gcdNum:=(z,w,i)->igcd(numer(z[i]),numer(w[i])):
  gcdDen:=(z,w,i)->igcd(denom(z[i]),denom(w[i])):
# -------------------------------- setting I/O
  if nargs=4 then
     xy:=1:
  else
     xy:=args[5]
  fi:
  Outtype:=(x,y)->gcdNum(x,y,xy):
# -------------------------------- main loop
  AC:=array(1..maxtau):
  Z:=array(0..length):
  Z[0]:=z0:
  for t to length do
    Z[t]:=F(Z[t-1])
  od:
  for tau to maxtau do
     AC[tau]:=0:
     for t to length-tau do
        AC[tau]:=AC[tau]+1/Outtype(Z[t],Z[t+tau])
     od:
     AC[tau]:=AC[tau]/(length-tau):
  od:
  return [seq(AC[t],t=1..maxtau)]
end:

#  h_______________________________________________________________General
#  h(x)		Logarithmic height of x (floating-point)
#  
h:=proc(x)
global H:
  return evalf(log(H(x)))
end:

#  isOK____________________________________________________________
#  is the periodic point OK?
isOK:=proc(z0,Fs,Cs,Xs)
global Code:
local c, n, z:
  z:=z0:
  for c in Cs do
    if Code(z,Xs) <> c then return false fi:
    z:=[Fs[c](z[1])-z[2],z[1]]
  od:
  return evalb(Code(z,Xs)=Cs[1])
end:

#  Jac_____________________________________________________________
#  The Jacobian of the list Js and code Cs
Jac:=proc(Js,Cs)
local c:
  evalm([[1,0],[0,1]]):
  for c in Cs do
     evalm([[Js[c],-1],[1,0]]&*%)
  od:
  return [[%[1,1],%[1,2]],[%[2,1],%[2,2]]]
end:

#  Map_____________________________________________________________
#  The map generated by the function list Fs and the code Cs
Map:=proc(z,Fs,Cs)
local  c:
  z:
  for c in Cs do
    [Fs[c](%[1])-%[2],%[1]]
  od:
  return %
end:

#  Mat_____________________________________________________________
#  The matrix of the function list Fs and the (periodic) code Cs
Mat:=proc(Fs,Cs)
global Map:
local  Aff, Lin, x,y, z:
   Aff:=Map([x,y],Fs,Cs);
   solve({Aff[1]=x,Aff[2]=y},{x,y}):
   z:=subs(%,[x,y]):
   Lin:=subs(x=(x+z[1]),y=(y+z[2]),Aff-z);
   [coeff(Lin[1],x,1),coeff(Lin[1],y,1)],
   [coeff(Lin[2],x,1),coeff(Lin[2],y,1)]:
   return [%]
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:

#  OrbitData___________________________________________________
#       Data (p-values, code) of orbit of F through z0, with given length
OrbitData:=proc(Fdata::list,z0,length::posint)
global Code, pvalue:
local n, nit, nP, t, z, F, P, L,
       pvalueSeq, CodeSeq, flag, zh:
  if nargs=3 then flag:=0 else flag:=args[4] fi:
  F,P,L:=op(Fdata):
  nP:=nops(P):
  pvalueSeq:=array(1..length,1..nP):
  CodeSeq:=array(1..length):
  z:=z0:
  CodeSeq[1]:=Code(z,L):
  if flag=0 then zh:=z else zh:=z[flag] fi:
  for n to nP do
     pvalueSeq[1,n]:=pvalue(zh,P[n])
  od:
  z:=F(z):
  for t from 2 while z<>z0 to length do
     CodeSeq[t]:=Code(z,L):
     if flag=0 then zh:=z else zh:=z[flag] fi:
     for n to nP do
        pvalueSeq[t,n]:=pvalue(zh,P[n])
     od:
    z:=F(z):
  od:
  nit:=min(length,t-1):
  [seq([seq(pvalueSeq[t,n],t=1..nit)],n=1..nP)],[seq(CodeSeq[t],t=1..nit)]:
  return [%]
end:

#  PeriodicPoint_______________________________________________
PeriodicPoint:=proc(Fs,Cs)
global Map:
local  x,y:
   Map([x,y],Fs,Cs);
   solve({%[1]=x,%[2]=y},{x,y}):
   return subs(%,[x,y])
end:

#  PoBoH___________________________________________________________General
#    points of bounded hieght
PoBoH:=proc(n::{posint,negint})
global Farey:
local Recip,FareySet,FareySetn,FareySetn1,BH,BHn,BHn1:
  if abs(n)=1 then
     [0,0],[1,0],[0,1]
  else
    Recip:=x->if x=0 then x else 1/x fi:
    if n > 0 then
       FareySet:={op(Farey(abs(n)))}:
       BH:=FareySet union map(Recip,FareySet):
       seq(seq([x,y],x=BH),y=BH):
    else
       FareySetn:={op(Farey(n))}:
       FareySetn1:={op(Farey(abs(n)-1))}:
       BHn:=FareySetn union map(Recip,FareySetn):
       BHn1:=FareySetn1 union map(Recip,FareySetn1):
       seq(seq([x,y],x=BHn),y=BHn),seq(seq([x,y],x=BHn),y=BHn1),seq(seq([x,y],x=BHn1),y=BHn):
    fi:
  fi:
  return {%} union {op(map(z->[z[1],-z[2]],[%]))} union {op(map(z->[-z[1],z[2]],[%]))} union {op(map(z->[-z[1],-z[2]],[%]))}
end:


#  pvalue__________________________________________________________General
#        p-adic value
pvalue:=proc(r::{integer,fraction,list},p::{integer,extended_numeric})
local  n, d, nd, k, pk, segno:
  if type(r,list) then
     if p=infinity then
        return max(pvalue(r[1],p),pvalue(r[2],p))
     else
        return min(pvalue(r[1],p),pvalue(r[2],p))
     fi:
  fi:
  if p=infinity then return log(max(abs(numer(r)),abs(denom(r)))) fi:
  if r=0 then return(infinity) fi:
  n:=numer(r):
  d:=denom(r):
  if igcd(n,p)=1 then
    nd:=d:
    segno:=-1
  else
    nd:=n:
    segno:=1
  fi:
  pk:=p:
  for k while igcd(nd,pk) = pk do pk:=pk*p od:
  return (k-1)*segno
end:

# GcdData______________________________________________________
#       gcd of pairs of points of orbit of F through z0, 
#       with given length and time-delay
GcdData:=proc(F,z0,length)
global Freq:
local All, gcdNum, gcdDen, gcdNumDen, 
      delay, k, n, nit, t, z, w, i,
      pNew,pNow,pOld,pTimes,
      Outtype, Outf, OutSeq:
# -------------------------------- functions
  gcdNum:=(z,w,i)->igcd(numer(z[i]),numer(w[i])):
  gcdDen:=(z,w,i)->igcd(denom(z[i]),denom(w[i])):
  gcdNumDen:=(z,w,i)->gcdNum(z,w,i)/gcdDen(z,w,i):
# -------------------------------- setting the output
  if nargs>3 then
    delay:=args[4]
  else
    delay:=1
  fi:
  if nargs>4 then
    Outtype:=args[5]
  else
    Outtype:=1
  fi:
  if Outtype=1 then
    Outf:=(x,y)->gcdNum(x,y,1):
  elif Outtype=11 or Outtype=12 then 
    Outf:=(x,y)->ifactors(gcdNum(x,y,1))[2]:
  elif Outtype=13 then 
    Outf:=proc(x,y) #--- list of prime divisors of gcd
      ifactors(gcdNum(x,y,1))[2]:
      map(z->z[1],%):
      return {op(%)}
    end:
  elif Outtype=2 then
    Outf:=(x,y)->gcdNumDen(x,y,1):
  elif Outtype=21 or Outtype=22 then 
    Outf:=(x,y)->ifactors(gcdNumDen(x,y,1))[2]:
  elif Outtype=3 then
    Outf:=(x,y)->gcdNum(x,y,2):
  elif Outtype=31 or Outtype=32 then 
    Outf:=(x,y)->ifactors(gcdNum(x,y,2))[2]:
  fi:
  OutSeq:=array(1..length):
# -------------------------------- main loop
  z,w:=z0,z0:
  for k to delay do z:=F(z) od:
  OutSeq[1]:=Outf(z,w):
  z,w:=F(z),F(w):
  for t from 2 while z<>z0 to length do
    OutSeq[t]:=Outf(z,w):
    z,w:=F(z),F(w):
  od:
  nit:=min(length,t-1):
  All:=[seq(OutSeq[t],t=1..nit)]:
# -------------------------------- output 
  if Outtype=1 or Outtype=3 then
    return All
  elif Outtype=11 or Outtype=31 then
    return All
  elif Outtype=12 or Outtype=32 then
    map(z->z[1],map(op,All));
    return Freq(%)
  elif Outtype=2 then
    return All
  elif Outtype=21 then
    return All
  elif Outtype=22 then
    map(z->z[1],map(op,All));
    return Freq(%)
  elif Outtype=13 then
    pTimes:=NULL:
    pOld:={}:
    for t to nit do
      pNow:=All[t]:
      if not (pNow subset pOld) then
        pNew:=pNow minus pOld:
        pOld:=pOld union pNow:
        pTimes:=pTimes,[t,pNew]:
      fi:
    od:
    return [pTimes]
  fi:
end:

# punit___________________________________________________________NumThy 
#      p-adic unit part of a rational
punit:=proc(r::{integer,fraction,list,set},p::integer)
global pvalue:
 if type(r,list) then
    return [punit(r[1],p),punit(r[2],p)]
 elif type(r,set) then
    return {punit(r[1],p),punit(r[2],p)}
 elif r=0 then
    return(0)
 else
    return(r/p^pvalue(r,p))
 fi:
end:

#  ZeroGaps________________________________________________________General
#  ZeroGaps(L)   seq of gaps between consecutive occurrences of 0
#  		   in the sequence L
#    FUNCTIONS: none
ZeroGaps:=proc(L::list)
local kis0:
 kis0:=(k,L)->if L[k]=0 then k else NULL fi:
 [seq(kis0(k,L),k=1..nops(L))]:
 [seq(%[k+1]-%[k],k=1..nops(%)-1)]:
 return %
end:


#  Tra_____________________________________________________________
#  Sequence of traces of the Jacobian of the list Js and code Cs
Tra:=proc(Js,Cs)
local c, n, N, Tr, M:
  N:=nops(Cs):
  Tr:=array(1..N):
  M:=evalm([[1,0],[0,1]]):
  for n to N do
     c:=Cs[n]:
     M:=evalm([[Js[c],-1],[1,0]]&*M):
     Tr[n]:=trace(M):
  od:
  return [seq(Tr[n],n=1..N)]
end:

#------------------ Zsigmondy set of an integer sequence
Zsig:=proc(L)
local j,AllPrimes,NewPrimes,PrimesNow,Z:
 AllPrimes:={}:Z:=NULL:
 for j to nops(L) do
    PrimesNow:={op(map(x->x[1],ifactors(L[j])[2]))}:
    NewPrimes:=select(x-> not member(x,AllPrimes),PrimesNow):
    if NewPrimes={} then Z:=Z,j fi:
    AllPrimes:=AllPrimes union NewPrimes:
 od:
 return [Z]
end:

#============================================================  MAPS

#  G______________________________________________________________
#  The continuous piecewise-affine map with 2 pieces
#  Parameters: aG (two slopes), 
#              xG (discontinuity),
#              xG (height at discontinuity: yG=G(xG))
#   then set:  bG:=[yG-aG[1]*xG,yG-aG[2]*xG]
G:=proc(z)
global Glist, xG:
local x,y:
   x,y:=op(z):
   if x<xG then 
      Glist[1](x) 
   else 
      Glist[2](x) 
   fi:
   return [%-y,x]
end:

#  ----------------------------- auxiliary functions for G 
Glist:=[proc(x) global aG,bG: return aG[1]*x+bG[1] end,
        proc(x) global aG,bG: return aG[2]*x+bG[2] end]:

#  F______________________________________________________________
#  The continuous piecewise-affine map with 3 pieces
#  Parameters: aF (3 slopes), 
#              xF (2 discontinuities), 
#              yF (height at leftmost discontinuity: F(xF[1])
#    then set: bF:=[yF-aF[1]*xF[1],yF-aF[2]*xF[1],yF-aF[2]*xF[1]+aF[2]*xF[2]-aF[3]*xF[2]]:
F:=proc(z)
global Flist, xF:
local x,y:
   x,y:=op(z):
   if x<xF[1] then 
      Flist[1](x) 
   elif x<xF[2] then 
      Flist[2](x) 
   else 
      Flist[3](x) 
   fi:
   return [%-y,x]
end:

Fdsp:=proc(z)
global Flist, xF, dsp:
local x,y:
   x,y:=op(z):
   if x<xF[1] then
      Flist[1](x)
   elif x<xF[2] then
      Flist[2](x)
   else
      Flist[3](x)
   fi:
   return [%-y,dsp*x]
end:


#  ----------------------------- auxiliary functions for F 
Flist:=[proc(x) global aF,bF: return aF[1]*x+bF[1] end,
        proc(x) global aF,bF: return aF[2]*x+bF[2] end,
        proc(x) global aF,bF: return aF[3]*x+bF[3] end]: