#C24.txt, April 21, 2025
Help24:=proc(): print(` FibMod(INI,m), RF(n) , FindPer(f,i),AllF(m,n) , GFperiod(n,x) `):
print(`GFperiodPer(n,x), EstimateGFperiod(n,x,K)`):
end:
with(combinat):
#FibMod(INI,m):the trajectory of the mapping [a,b]->[b,a+b mod m] starting with INI
FibMod:=proc(INI,m) local pt,L,i:
pt:=INI:
L:=[]:
while not member(pt,L) do
 L:=[op(L),pt]:
 pt:=[pt[2],pt[1]+pt[2] mod m]:
od:
for i from 1 to nops(L) while L[i]<>pt do od:
[[op(1..i-1,L)],[op(i..nops(L),L)]]:
end:


#A random function from {1,...,n} to {1,...,n}
RF:=proc(n) local i,ra:
ra:=rand(1..n):
[seq(ra(),i=1..n)]:
end:

#FindPer(f,i): inputs a function from {1,..,n} to {1,...,n} where n:=nops(f):
#and outputs the pair [L1,L2] where L1 is the beginning segment and L2 is the cycle
FindPer:=proc(f,i) local L,pt,i1:
pt:=i:
L:=[]:
while not member(pt,L) do
 L:=[op(L),pt]:
 pt:=f[pt]:
od:
for i1 from 1 to nops(L) while L[i1]<>pt do od:
[[op(1..i1-1,L)],[op(i1..nops(L),L)]]:
end:

#AllF(m,n): all lists of length n with entries from {1,..,m}
AllF:=proc(m,n) local S,s,i:
if n=0 then
   RETURN({[]}):
fi:
S:=AllF(m,n-1):
{seq(seq([op(s),i],s in S),i=1..m)}:
end:

#GFperiod(n,x):the polynomial whose coeff. of x^j is the prob. that if you
#take a random function from [1,n] to [1,n] and a random starting point
#the length of the ultimate cycle is j
GFperiod:=proc(n,x) local S,i,s,f,av,sd:
S:=AllF(n,n):
f:=add(add(x^nops(FindPer(s,i)[2]),i=1..n),s in S)/(nops(S)*n):
av:=subs(x=1,diff(f,x)):
sd:=sqrt(subs(x=1,diff(x*diff(f,x),x))-av^2):
[f,[av,sd]]:

end:

#GFperiodPer(n,x):the polynomial whose coeff. of x^j is the prob. that if you
#take a random permutation from [1,n] to [1,n] and a random starting point
#the length of the ultimate cycle is j, followed by the average and standard-deviation
GFperiodPer:=proc(n,x) local S,i,s,f,av,sd:
S:=permute(n):
f:=add(add(x^nops(FindPer(s,i)[2]),i=1..n),s in S)/(nops(S)*n):
av:=subs(x=1,diff(f,x)):
sd:=sqrt(subs(x=1,diff(x*diff(f,x),x))-av^2):
[f,[av,sd]]:
end:

#EstimateGFperiod(n,x,K)
#samples K random functions from [1,n] to [1,n] and a random starting point. and finds the prob. gen. function for this run
#followed by the sample average and standard-deviation. Try:
#EstimateGFperiod(20,x,100);
EstimateGFperiod:=proc(n,x,K) local ra,i,j,f,av,sd:
ra:=rand(1..n):
f:=evalf(add(x^nops(FindPer(RF(n),ra())[2]),j=1..K)/K):
av:=subs(x=1,diff(f,x)):
sd:=sqrt(subs(x=1,diff(x*diff(f,x),x))-av^2):
[f,[av,sd]]:
end: