# OK to post homework
# Aurora Hiveley, 4/28/26, Assignment 26

Help:=proc(): print(`necklaceColors(n)`):   end:

### Problem 1 
# Check that procedure CubeGp() is correct. If it is not, correct it.

# now updated! X was the same as Y, now X is distinct 

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

# f := CIP(CubeGp(),x);
#          6         2     2         2              3         2
# f := x[1]  + 3 x[1]  x[2]  + 6 x[1]  x[4] + 6 x[2]  + 8 x[3] 

# seq(subs({seq(x[i]=n, i=1..4)},f)/nops(CubeGp()), n=0..8);
#            0, 1, 10, 57, 240, 800, 2226, 5390, 11712

# in OEIS as entry A047780 



### Problem 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.

necklaceColors := proc(n) local S,G,f,i,c:
# generate the symmetry group, i.e. the cyclic group of order n  
S := { seq([seq(i..n),seq(1..i-1)], i=1..n) }: 
G := GenGp(S): 
f := CIP(G,x):
subs({seq(x[i]=c, i=1..n)}, f):
end:





### copied from C26.txt 


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

with(combinat):

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:=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/ only rotations
CubeGp:=proc() local X,Y,Z:
#[U,L,F,R,B,D]=[1,2,3,4,5,6]
# X:=[5,2,1,4,6,3]: --> was same as Y, so let's change it:

X:=[4,1,3,6,5,2]:
Y:=[5,2,1,4,6,3]:
Z:=[1,5,2,3,4,6]:
GenGp({X,Y,Z}):
end:


### copied from C3.txt 

#C3.txt, Jan. 29, 2026
Help3:=proc(): print(`FP(pi), Der(n), d(n) , ExtractCycle(pi,i) , CycDec(pi) `): end:

#FP(pi): The number of fixed points of per. pi
FP:=proc(pi) local n, i,co:
n:=nops(pi):

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

FP([2,1,4,3,5,6]);
 

# Der(n): The set of derangements (i.e. permutations w/o fixed points) of {1, ..., 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:

#a(n):=nops(permute(n))
#
#Solutions SOl. 
#a(n)-n*a(n-1)=0, a(0)=1 HOMOG. linear (first-order) recurrence
#
#d(n):=|Der(n)|
#

d:=proc(n) option remember: if n=1 then 0 else n*d(n-1)+(-1)^n 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: