# OK to post homework
# Aurora Hiveley, 2/19/26, Assignment 8

Help:=proc(): print(`AntiFoata(pi)`):   end:

with(combinat):

### Problem 1
# Show that for all primes p less than 43, if you compute (xnP(p-1,p)+p-1)*xnP(p-1,p) mod p
# you get 0, not enabling you to deduce that it stops producing integers, but for p=43 it breaks down, 
# concluding that xn(43) is NOT an integer (even though it is too big to compute directly)

# S := {}:
# p := 2:
# while p < 43 do 
#     S := S union {(xnP(p-1,p)+p-1)*xnP(p-1,p) mod p}:
#     p := nextprime(p):
# od:
# S:  # get S = {0}, as desired!!


### Problem 2
# verify that for all permutations of size ≤ 7, nops(LtoR(Foata(pi))=nops(CycDec(pi)))

# {seq( seq( evalb( nops(LtoR(Foata(pi)))=nops(CycDec(pi)) ), pi in permute(n) ), n=1..7)};
# returns {true} -- yay!!


### Problem 3
# this means that for all integers n and k ≤ n, the number of permutations with k cycles equals 
# the number of permutations with k Left-To-Right maxima, since Foata is a bijection. 



### Problem 4
# AntiFoata(pi): inverse of Foata(pi)

## fixed from hw6 
Foata := proc(pi) local C,c:
C := sort(convert( sort({seq( CFC(c), c in CycDec(pi) )}), list), (a,b) -> a[1] < b[1]):
[seq(op(c), c in C)]:
end:


AntiFoata := proc(pi) local C,M,i:
M := LtoR(pi):
C := [seq( pi[M[i]..M[i+1]-1], i=1..nops(M)-1), pi[M[-1]..-1]]: # set of cycles
CycToPer(convert(C,set)):
end:

# for every permutation pi of size ≤ 7, AntiFoata(Foata(pi))=pi and Foata(AntiFoata(pi))=pi

# {seq( seq( evalb(AntiFoata(Foata(pi))=pi), pi in permute(n) ), n=1..7)};
# returns {true}

# {seq( seq( evalb(Foata(AntiFoata(pi))=pi), pi in permute(n) ), n=1..7)};
# returns {true}



### Problem 5
# xnk(n,k), xnkC(n,k), and xnkP(n,k,p): give the first value (and those mod p) of the sequence defined by the recursion
# xnk(n,k)= (1+xnk(0,k)^k+...+xnk(n-1,k)^k)/n

#xnk(n,k): The solution of the recurrence  xn(n)=(1+add(xn(i)^2,i=0..n-1))/n:
xnk:=proc(n,k) local i: option remember:
if n=0 then
 1:
 else
 (1+add(xnk(i,k)^k,i=0..n-1))/n:
fi:
end:

#xnkC(n,k): clever version of xn(n) only having to remember what happened yesterday

# xnk(n,k)= (1+xnk(0,k)^k+...+xnk(n-1,k)^k)/n
#         = (1+xnk(0,k)^k+...+xnk(n-2,k)^k)/n + (xnk(n-1,k)^k)/n
#         = (n-1)/n* (1+xnk(0,k)^k+...+xnk(n-2,k)^k)/(n-1) + (xnk(n-1,k)^k)/n
#         = (n-1)/n* xnk(n-1,k) + (xnk(n-1,k)^k)/n

xnkC:=proc(n,k) option remember:
if n=1 then
2:
else
(xnkC(n-1,k)+n-1)*xnkC(n-1,k)/n:
fi:
end:

#xnkP(n,k,p): For a prime p, a fast way to compute xnk(n,k) mod p for n<p, as long as it is an integer
xnkP:=proc(n,k,p) option remember:
if n=1 then
2:
else
(xnkP(n-1,k,p)+n-1)*xnkP(n-1,k,p)*n^(-1) mod p:
fi:
end:


# [Note that xnk(n,2) equals xn(n)]
# tested: {seq( evalb(xnk(n,2)=xn(n)), n=1..10)}; --> returns {true}

# verify that for k=3, and k=4, the first 15 terms are integers.
# {seq( type(xnk(n,3), integer), n=1..15)}; --> returns {true}
# {seq( type(xnk(n,4), integer), n=1..15)}; --> returns {true}



### Problem 6
# [Optional challenge, 5 dollars, no peeking, I know the answer]
# Find the smallest n such that xnk(n,3) is NOT an integer (similar to 43 for k=2)

# p := 2:
# while (xnkP(p-1,3,p)+p-1)*xnkP(p-1,3,p) mod p = 0 and p < 300 do   # safeguard to prevent infinite loop     
#     p := nextprime(p):
# od:

# p; --> returns 43 again??
# although oeis entry A005166 claims its 89


### Problem 7
# [Optional challenge, 6 dollars, no peeking, I don't know the answer]
# Find the smallest n such that xnk(n,4) is NOT an integer

# p := 2:
# while (xnkP(p-1,4,p)+p-1)*xnkP(p-1,4,p) mod p = 0 and p < 300 do   # safeguard to prevent infinite loop     
#     p := nextprime(p):
# od:

# p; --> once again, returns 43
# oeis entry A005166 claims its 97



### copied from C8.txt 


HelpOld:=proc(): print(`LtoR(pi), ExtractCycle(pi,i), CycDec(pi)`): end:

Help8:=proc(): print(`xn(n), xnC(n), xnP(n,p), Foata(pi) `): end:

#xn(n): The solution of the recurrence  xn(n)=(1+add(xn(i)^2,i=0..n-1))/n:
xn:=proc(n) local i: option remember:
if n=0 then
 1:
 else
 (1+add(xn(i)^2,i=0..n-1))/n:
fi:
end:

#xnC(n): clever version of xn(n) only having to remember what happened yesterday
xnC:=proc(n) option remember:
if n=1 then
2:
else
(xnC(n-1)+n-1)*xnC(n-1)/n:
fi:
end:

#xnP(n,p): For a prime p, a fast way to compute xn(n) mod p for n<p, as long as it is an integer

xnP:=proc(n,p) option remember:
if n=1 then
2:
else
(xnP(n-1,p)+n-1)*xnP(n-1,p)*n^(-1) mod p:
fi:
end:

#LtoR(pi): The list places that are larger than anything to the left pi=[2,1,4,3]
LtoR:=proc(pi): local n, L, i, ma: n:=nops(pi): L:=[1]: ma:=pi[1]: for i from 2 to n do if pi[i]>ma then L:=[op(L), i ]; ma:=pi[i]; end if; end do; L; end proc:

#ExtractCycle(pi,i): The cycle corresponding to the member i of the permutation pi For example ExtractCycle([2,1,4,3],1)=[1,2]
ExtractCycle:=proc(pi,i): local C, ng: C:=[i]: ng:=pi[i]: while ng<>i do C:=[op(C),ng]: ng:=pi[ng]: end do; C; end proc:

CycDec:=proc(pi): local n, i, StillToDo, S, C, ng: n:=nops(pi): StillToDo:={seq(i,i=1..n)}: S:={}: while StillToDo<>{} do ng:=StillToDo[1]:
C:=ExtractCycle(pi,ng): S:=S union {C}: StillToDo:=StillToDo minus {op(C)}: end do; S; end proc:

#CFC(C): inputs a list of numbers coming from a cycle and outputs the equivalent cycle where the largest entry is the first. CFC([4,6,1])= [6,4,1]
CFC:=proc(C): local k: k:=max[index](C): [op(k..nops(C),C),op(1..k-1,C)]; end proc:

#CycToPer(C): The reverse of CycDec(pi). Given a permutation in cycle structure, outputs it in 1-line notation. For example CycToPer({[1,2],[3,4]})=[2,1,4,3]
CycToPer:=proc(C): local n,i,j,T: n:=add(nops(C[i]),i=1..nops(C)): for i from 1 to nops(C) do for j from 1 to nops(C[i])-1 do T[C[i][j]]:=C[i][j+1]:
end do; T[C[i][-1]]:=C[i][1]: end do; [seq(T[i],i=1..n)]; end proc: