# OK to post homework
# Lucy Martinez, 04-29-2026, Assignment 26

with(combinat):


# Question 1:
# Check that procedure CubeGp() is correct. If it is not, correct it.
# Using the (possibly corrected version), Find an expression for the number of ways
#  to color a cube with c colors (up to rotational symmetry). Is it in the OEIS?
# ANSWER:
# [seq(CubeColorPoly(c),c=1..10)]
# returns
#  [1, 10, 57, 240, 800, 2226, 5390, 11712, 23355, 43450]
# YES! It is in the OEIS:  
#  It is A047780 with the following description
#  Number of inequivalent ways to color faces of a cube using at most n colors. 
CubeColorPoly:=proc(c) local i,S,pi:
S:=CubeGp():
1/nops(S)*(add(c^nops(CycDec(pi)),pi in S)):
end:



# Question 2:
# Write a procedure that inputs a positive integer n (not necessarilty prime)
#  and outputs the polynomial in c, that tells you the number of ways to color a necklace (w/o clasp)
#  of n beads with c colors.
#We first need to check for each member in the group (in this case the set is {0,1,2,...,n-1}
# how many cycles does each member have? 
# For example: Let n=6 then we look at the different ways to shift the beads so in this example,
#  we can shift the beads by 2 to the right so 
#  2->4->0->2 which means the orbit {0,2,4} has a cycle of length 3
#  Now 1->3->5->1 which means the orbit {1,3,5} has a cycle of length 3
#  Here, the necklace beads are grouped into:
#  Cycle 1: {0,2,4}
#  Cycle 2: {1,3,5}
#  which means each bead in each cycle will have to get the same color
# Notice that moving forward by 4 (so bead 0 goes to position 5) produces the following cycles
#  Cycle 1: {0,4,2}
#  Cycle 2: {1,5,3}
# So we get the "same" rotations but these are different permutations!
# In conclusion, the number of cycles for each element in G=Z_n is gcd(n,k) where k=0,1,...,n-1
# We really are doing : i+k (mod n) and trying to check when do we get back to the start
#  for each k=0,1,...,n-1
NeckColorPoly:=proc(n,c) local k:
add(c^gcd(n,k),k=0..n-1)/n:
end:


################################From previous classes:
#April 23, 2026 C24.txt

Help26:=proc():
 print(`BTseq(N), BTseqE(N), CycDec(pi), WtPi(pi,x)`):
 print(`CIP(G,x), Mul(pi,sig), GenGp(S), CubeGp()`):
end:



BTseq:=proc(N) local x,f,i:
f:=x:
for i from 1 to N+1 do
f:=expand(x+f^2):
f:=taylor(f,x=0,N+2):
f:=add(coeff(f,x,i)*x^i,i=1..N):
od:
[seq(coeff(f,x,i),i=1..N)]:
end:


BTseqE:=proc(N) local x,f,i,i1:
f:=x:
for i from 1 to N+1 do
f:=expand(x+(f^2+subs(x=x^2,f))/2):
f:=taylor(f,x=0,N+2):
f:=add(coeff(f,x,i1)*x^i1,i1=1..N):
od:
[seq(coeff(f,x,i1),i1=1..N)]:
end:


#WtPi(pi,x): The weight of pi
WtPi:=proc(pi,x) local L,i:
L:=CycDec(pi):
mul(x[nops(L[i])],i=1..nops(L)):
end:



#CIP(G,x): The cycle index polynomial of G
CIP:=proc(G,x) local g:
add(WtPi(g,x),g in G):
end:


#Mul(pi,sig): The permutation pi times the permutation sig from left to right
Mul:=proc(pi,sig) local i:
[seq(sig[pi[i]],i=1..nops(pi))]:
end:


#GenGp(S): The group generated by S
GenGp:=proc(S) local G1,G2,s,g1,g2:
G1:=S:
G2:=G1 union {seq(seq(Mul(s,g1),g1 in G1),s in S)}:
while G1<>G2 do
G1:=G2:
G2:=G2 union {seq(seq(Mul(s,g2),g2 in G2),s in S)}:
od:
G2:
end:


#CubeGp(): The cube symmetry group w/o rotations
CubeGp:=proc() local X,Y,Z:
#[U,L,R,R,B,D]=[1,2,3,4,5,6]
X:=[5,2,1,4,6,3]:
Y:=[5,2,1,4,6,3]:
Z:=[1,5,2,3,4,6]:
GenGp({X,Y,Z}):
end:


################################From previous classes:
# C3.txt, Jan. 29, 2026

Help3:=proc() :
 print(`FP(pi), WtE(n,x), Der(n), d(n), ExtractCycle(pi,i), CycDec(pi)`):
end:


# FP(pi): inputs a permutation pi of {1,...,n}, where n:=nops(pi) 
#  and outputs the number of places of fixed points of pi
#  For example NuFP([2,1,3,5,4])=1
FP:=proc(pi) local i,n,count:
n:=nops(pi):
count:=0:

if not type(pi,list) then
 print(`The input should be a list`):
 return(FAIL):
fi:

if {seq(type(pi[i],integer),i=1..n)}<>{true} then
 print(`The elements in the list should be integers`):
 return(FAIL):
fi:

for i from 1 to n do 
 if pi[i] = i then 
  count ++:
 fi:
od:
count:
end:


# WtE(n,x): inputs a non-negative integer and a variable x 
#  and outputs the sum of x^NuFP(pi) summed over all
#  permutations of length n.
WtE := proc(n,x) local pi:
add(x^NuFP(pi), pi in permute(n)):
end:


# Der(n): The set of derangements (i.e. permutations without fixed pts)
#  of {1,2,...,n}
Der:=proc(n) local S,pi,DE:
S:=permute(n):
DE:={}:
for pi in S do
 if FP(pi)=0 then
  DE:=DE union {pi}:
 fi:
od:
DE:
end:


#Linear first order recurrence (homogeneous)
#if a(n):=nops(permute(n))
# then a(n)=n*a(n-1)
# a(n)-n*a(n-1)=0 with a(0)=1

# d(n):=|Der(n)|
# [seq(L[n]-n*L[n-1],n=2..nop(L))]=[1, -1, 1, -1, 1, -1]
# This means Der(n)-n*Der(n-1)-(-1)^n=0
# so Der(n)=n*Der(n-1)+(-1)^n

#d(n): the number of derangements of length n
d:=proc(n) option remember:
if n=1 then
 0:
else
 n*d(n-1)+(-1)^n:
fi:
end:

#limit of d(n)/n! as n goes to infinity is 1/e


# Cycle structure
# [3,1,4,2] ONE-LINE-NOTATION

# Two-Line Notation
# 1 2 3 4
# 3 1 4 2

# 3142=(1342)
# above: 1 goes to 3, 3 goes to 4, 4 goes to 2, 2 goes to 1
# Cycle decomposition of (3,1,4,2) is (1342)


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

# How many cyclic permutations are there of length n?
# (n-1)!
# So the probability of getting a cyclic perm is (n-1)!/n!=1/n


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

##############################
# C2.txt, Jan. 26, 2026

Help2:=proc() :
 print(`redu(L), SubSeqs3(L,k), Contain3(pi,sig), AvoidPer1(n,sig)`):
end:

#permutations: arranging objects with orders
with(combinat):

#if a(n):=number of permutations of {1,2,...,n}
# conjecture: a(n)/a(n-1)=n
# first-order recurrence
# so a(n)=n*a(n-1) where a(0)=1
# 
# A different recurrence: Fibonacci numbers
# F(n)=F(n-1)+F(n-2), second-order linear recurrence
# 
# Back to permutations:
# n!=1*2*3*...*(n-1)*n

#Per(n): the set of permutations of length n
# Per(n)=Union of P(n-1) with n insertd at the first place
# T:=Per(n)->{1,2,...,n} X Per(n-1)
# T(pi)=[i,pi']
#  where i=location of n and
#  pi' is the permutation of {1,2,...,n-1} obtained by deleting n
# Implementation in Maple:
#  B:=proc(pi) : inputs a permutation pi of {1,2,...,n}, outputs [i,pi']

# Pattern avoidance
# 123-avoiding permutation
# i.e. you never have a subsequence of lenght 3 that is increasing


#redu(L): inputs a list of distinct numbers and outputs
# its reduction according to their order
# For example: redu([5,9,1])=[2,3,1]
# redu([Pi,e])=[2,1]
redu:=proc(L) local n, L1,T,i:
n:=nops(L):
L1:=sort(L):
for i from 1 to n do
 T[L1[i]]:=i 
 #the first number in the sorted version of L1 will be called 1
od:
[seq(T[L[i]],i=1..n)]:
end:


#SubSeqs3(L): The set of subsequences of the list L of length k
# For example: SubSeqs3([1,6,2,4])={[1,6,2],[1,6,4],[1,2,4],[6,2,4]}
SubSeqs3:=proc(L) local n,i1,i2,i3,S:
n:=nops(L):
S:={}:
for i1 from 1 to n do
 for i2 from i1+1 to n do
  for i3 from i2+1 to n do
   S:= S union {[L[i1],L[i2],L[i3]]}:
  od:
 od:
od:
S:
end:


#Contain3(pi,sig): inputs a permutation pi and outputs if pi
# contains the pattern sig
# Contain3([2,1,3,4],[1,2,3]) returns true
# Contain3([4,3,2,1],[1,2,3]) returns false
Contain3:=proc(pi,sig) local n,S,s:
n:=nops(pi):
if nops(sig)<>3 then
 return(FAIL):
fi:
S:=SubSeqs3(pi):
for s in S do
 if redu(s)=sig then
  return(true):
 fi:
od:
false:
end:


#AvoidPer1(n,sig): inputs a pos. integer n and pattern of length 3
# outpus the subset of permute(n) that avoid the pattern sig
AvoidPer1:=proc(n,sig) local A,pi,G:
A:=permute(n):
G:={}:
for pi in A do
 if not Contain3(pi,sig) then
  G:=G union {pi}:
 fi:
od:
G:
end: