# Ok to post homework
# Lucy Martinez, 04-24-2025, Assignment 24


######################
#Problem 1: Briefly describe the progress on your team's project
# My team has been reading and trying to understand the 
# paper of the algorithm we want to implement. We had 3 group
# meetings this past week and have been discussing the definitions
# from the paper.




#Problem 2: Write a procedure CompareFibMod(INI,M)
# that for m from 2 to M compares FibMod(INI,m)
# to EstimateGFperiod(m^2,x,10000)[2][1].
# Is the former ever larger than the latter? when INI=[0,1] and M=300?
# I ran this for M=100 since M=300 was taking too long
# The following output means 87 times
# the former was larger than the latter and 12 times it was not
# Answer: 87*x^true + 12*x^false

CompareFibMod:=proc(INI,M) local m,C,f1,f2,f:
f:=0:
for m from 2 to M do
 f1:=nops(FibMod(INI,m)[2]):
 f2:=EstimateGFperiod(m^2,x,10000)[2][1]:
 if f1>f2 then
  f:=f+x^true:
 else
  f:=f+x^false:
 fi:
od:
f:
end:


####################################from previous classes:
#C24.txt: April 21, 2025

Help24:=proc():
print(`FibMod(INI,m), RF(n), FindPer(f,i), AllF(m,n)`):
print(`GFperiod(n,x), GFperiodPer(n,x), EstimateGFperiod(n,x,K)`):
end:

with(combinat):

#Vertices {[a,b], a<=a,b<=m-1 [a,b]->[b,a+b mod m]}

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


#RF(n): 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 pt,L,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 probability
# 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: