# OK to post homework
# Aurora Hiveley, 4/22/25, Assignment 24

Help := proc(): print(`CompareFibMod(INI,M)`): end:

with(combinat):


### project update 

# group met to discuss the algorithm implementation outlined in the paper. it's quite dense so we each had to read it several times. 
# i found some supplemental papers that seemed to help fill in some of the details. implementation of smaller parts of the algorithm
# (helper procedures, mainly) has started, and pablo was kind enough to write up a skeleton of our overall procedure in pseudocode 
# and with data types identified so that we are all on the same page about how our graphs, spanning trees, vertex labels, etc. are encoded. 
# i also started working on a BFS procedure that will help us build an initial spanning forest to start the first iteration of the algorithm.


# CompareFibMod(INI,M): for m from 2 to M compares FibMod(INI,m) to EstimateGFperiod(m^2,x,10000)[2][1] 

CompareFibMod := proc(INI,M) local x,f,a,b,m:
f := 0:

for m from 2 to M do 
    a := nops(FibMod(INI,m)[2]):
    b := EstimateGFperiod(m^2,x,10000)[2][1]: 

    f := f + x[evalb(a<b)]:
od:

f:
end:

# Is the former ever larger than the latter? when INI=[0,1] and M=300?
# CompareFibMod([0,1],300):
# 300 took over 4 hours to run so i cancelled the execution 
# CompareFibMod([0,1],100):
# output: 87*x[false] + 12*x[true] 
# so it is possible for either one to be larger, although typically the estimation is smaller 








### copied from C24.txt 

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