#Please do not post homework
#Guy Adami, 2026-02-22, Assignment 8


#—————

with(combinat)

#—————

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:

#—————

#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]:
	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}:
			#S\T : S minus T 
		StillToDo:=StillToDo minus {op(C)}:
	od:
	S:
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 is the current world record
	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:


#—————

Foata:= proc(pi) local dec,cfc,n,i,P,phi:
	dec:= CycDec(pi):
	n:=nops(dec):
	cfc:=[]:
	for i from 1 to n do
		cfc:= [op(cfc),CFC(dec[i])]
	od:
	
	P:= sort(cfc):

	n:= nops(P):
	phi:=[]:
	for i from 1 to n do
		phi:= [op(phi),op(P[i])]:
	od:
	phi:
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:

#——————————————————————————————————————————————————

#Question 1

P := [seq(ithprime(i), i = 1 .. 14)]:

MODP := proc(P) local p, R, r; 
	R:= []:
	for p in P do 
		r := (xnP(p - 1, p) + p - 1)*xnP(p - 1, p) mod p:
		R := [op(R), r]:
	od; 
	R;
end:

#This list returns all 0s until the 14th prime=43


#—————

#Question 2

PI:= seq(permute(i),i=1..7):


LToR_EQ_CycDec := proc(PI) local pi, countfalse:
	countfalse:=0:
	for pi in PI do
		if nops(LToR(Foata(pi)))<> nops(CycDec(pi)) then
			countfalse:=countfalse + 1:
		fi:
	od:
	if countfalse <> 0 then
		RETURN(false);
	elif countfalse = 0 then
		RETURN(true);
	fi:
end:

#This function returned true, so this equality hold for all permutations of length <= 7

#—————

#Question 3

#Foata means putting a CycDec into canonical form will take any subcycle andalways put the greatest value at the start.
#By ordering it such that each subcycle of a permutation is ordered in ascending order of its first element
#means that for each sub cycle, the greatest value is always first, 
#in addition to having the following cycle begin with a greater value than the previous. 
#This will create a form of { [i, x<i …, i’>i, x’<i’…] }.
#As such, LtoR will only be incremented every time there is a new sub cycle by how it is ordered
#This means that the nops of LtoR and CycDec are equal

#—————

#Question 4

#Obtains the individual cycle components of a permutation assuming Foata notation
ObtainFoataCycles := proc(pi) local LR, i, CycPi; 
	LR := LtoR(pi); 
	CycPi := []; 
	for i to nops(LR) - 1 do  
		CycPi := [op(CycPi), pi[LR[i] .. LR[(i+1)] - 1]]; 
	od:
	CycPi := [op(CycPi), pi[LR[nops(LR)] .. nops(pi)]]:
	CycPi; 
end:

#Uses the cycles extracted from Foata notation and rearranges a list based on the cycles
AntiFoata := proc(pi) local i, j, CycPi, PI:

	CycPi := ObtainFoataCycles(pi):
	PI := pi:

	for i to nops(CycPi) do 
		for j to nops(CycPi[i]) do 
			if j = nops(CycPi[i]) then 
				PI[CycPi[i][j]] := CycPi[i][1]: 
			else 
				PI[CycPi[i][j]] := CycPi[i][j + 1]: 
			fi: 
		od: 
	od: 
	PI; 
end:


#—————

#Question 5


#xnk(n,k): The solution of the recurrence  xn(n)=(1+add(xn(i)^k,i=0..n-1))/n:

xnk:=proc(n) local i: option remember:
	if n=0 then
 		1:
 	else
 		(1+add(xn(i)^k,i=0..n-1))/n:
	fi:
end:


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

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


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

xnkP:=proc(n,k,p) option remember:
	if n = 0 then
		1:
	elif n=1 then
		2:
	else
		(xnkP((n-1),k,p)+n-1)*xnkP((n-1),k,p)/n mod p:
	fi:
end:


#Using a modifed MODP function, I was able to check if the result was zero or not which would correspond to values being prime

MODP_k := proc(P, k) local p, R, r; 
	R := []; 
	for p in P do 
		r := (xnkP(p - 1, k, p) + p - 1)*xnkP(p - 1, k, p) mod p; 
		R := [op(R), r]; 
	od:
	R;
end: