# OK to post homework
# Lucy Martinez, 02-19-2026, Assignment 8

with(combinat):


# Question 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) 
# Proof:
# Recall that n*xn(n) = xn(n-1)(n-1+xn(n-1))
#
# Assume that (p-1)*xnP(p,p) is an integer for p<=43 where p is a prime.
# Then, p*xnP(p,p) mod p= (p-1 + xnP(p-1,p))*xnP(p-1,p) mod p
# But p*xnP(p,p) mod p = 0 
# therefore
# (p-1 + xnP(p-1,p))*xnP(p-1,p) mod p = 0.





# Question 4:
# Write a procedure AntiFoata(pi) that is the inverse of Foata(pi)
#  that for every permutation pi of size <= 7
#  AntiFoata(Foata(pi))=pi and Foata(AntiFoata(pi))=pi
AntiFoata:=proc(pi) local L,C,i,n,len:
n:=nops(pi):
L:=LtoR(pi):
len:=nops(L):
C:=[]:
for i from 1 to nops(L)-1 do
 C:=[op(C), [op(L[i]..L[i+1]-1,pi)]]:
od:
C:=[op(C),[op(L[len]..n,pi)]]:
CycToPer(C):
end:


# Question 5:
# Write procedures xnk(n,k), xnkC(n,k), and xnkP(n,k,p)
#  that 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
#  [Note that xnk(n,2) equals xn(n)]
# Verify that for k=3, and k=4, the first 15 terms are integers
#
# ANSWER: for k=3, I ran L:={seq(denom(xnkC(n,3)),n=1..15)}; which returns {1}
#  which means that the denominator is always 1 so the values are integers
# 
# For k=4, L:={seq(denom(xnkC(n,4)),n=1..15)}; returns {1}
#  which means that the denominator is always 1 so the values are integers
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:


#Notice that n*xnk(n,k)=1+xnk(0,k)^k+...+xnk(n-1,k)^k
# (n-1)*xnk(n-1,k)=1+xnk(0,k)^k+...+xnk(n-2,k)^k
# SO, n*xnk(n,k)-(n-1)*xnk(n-1,k)=xnk(n-1,k)^k
# this implies that
# n*xnk(n,k)=(n-1)*xnk(n-1,k)+xnk(n-1,k)^k=xnk(n-1,k)(n-1+xnk(n-1,k)^(k-1))
# xnk(n,k)=xnk(n-1,k)(n-1+xnk(n-1,k)^(k-1))/n
#xnkC(n,k): clever version of xnk(n,k)
xnkC:=proc(n,k) local i:
option remember:
if n=1 then
 2:
else
 xnkC(n-1,k)*(n-1+xnkC(n-1,k)^(k-1))*n^(-1):
fi:
end:


#xnkP(n,k,p): clever version of xnkC(n,k) mod p
xnkP:=proc(n,k,p) local i:
option remember:
if n=1 then
 2:
else
 xnkC(n-1,k)*(n-1+xnkC(n-1,k)^(k-1))*n^(-1) mod p:
fi:
end:



##################################
#### From my homework from before:
Foata:=proc(pi) local cycle, CanForm,C, sig,i,L,s,newsig:
cycle:=CycDec(pi):
CanForm:={}:
for C in cycle do
 CanForm:=CanForm union {CFC(C)}:
od:
L:=[]:
for C in CanForm do
 L:=[op(L),C[1]]:
od:
L:=sort(L):

sig:=[]:

for i from 1 to nops(L) do
 for C in CanForm do
  if L[i]=C[1] then
   sig:=[op(sig), C]:
  fi:
 od:
od:

newsig:=[]:
for s in sig do
 newsig:=[op(newsig),op(s)]:
od:
newsig:
end:


###Other classes:
#added after class: 
#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]:
  od:
   T[C[i][-1]]:=C[i][1]:
  od:

[seq(T[i],i=1..n)]:
end:


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


# ExtractCycle(pi,i): The cyle 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]:
od:
C:
end:


#CycDec(pi): The full cyclic decomposition of pi
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)}:
od:
S:
end:


#LtoR(pi): The list of places that are larger than anything to the left
# For example: if pi=[2,1,4,3]
#  then LtoR(pi)={1,3}
LtoR:=proc(pi) local i,L,ma,n:
n:=nops(pi):
L:=[1]:
#ma is the current world record (the max)
ma:=pi[1]:
for i from 2 to n do
 if pi[i]>ma then
  L:=[op(L),i]:
  ma:=pi[i]:
 fi:
od:
L:
end:

##################################
#### From before:

# C8.txt, Feb. 16, 2026

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

#Gobel's sequence, tricks you to conjecture that
# it is always an integer

#x_n=(1+x_0^2+...+x_(n-1)^2)/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:

#Notice that n*xn(n)=1+x_0^2+...+x_(n-1)^2
# (n-1)*xn(n-1)=1+x_0^2+...+x_(n-2)^2
# n*xn(n)-(n-1)*xn(n-1)=x_(n-1)^2
# so
# n*xn(n)=(n-1)*xn(n-1)+x_(n-1)^2=xn(n-1)(n-1+xn(n-1))
# xn(n)=xn(n-1)(n-1+xn(n-1))/n
#xnC(n): clever version of xn(n)
xnC:=proc(n) local i:
option remember:
if n=1 then
 2:
else
 xnC(n-1)*(n-1+xnC(n-1))*n^(-1):
fi:
end:

#The sequence above is still slow. It will take a long time for n=40
# We have to prove that there is a place p, such that it is
# integers until p-1 but stops being an integer 
#xnP(n,p): clever version mod p
xnP:=proc(n,p) local i:
option remember:
if n=1 then
 2:
else
 xnP(n-1,p)*(n-1+xnP(n-1,p))*n^(-1) mod p:
fi:
end:

#Assume that xn(43) is an integer n*xn(n)=xn(n-1)(n-1+xn(n-1))
# then (xn(42)+42)xn(42) mod 43 =43*xn(43)  mod 43 = 0
# This means (xn(42)+42)xn(42)=0
# (33+42)*42 mod 43;