#OK to post homework
#Austin DeCicco, 5/3/26, Assignment 26

#BEGIN REFORMATTED C3
############################################################################################################################################################

#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; end if; end do; co; end proc:

# 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}; end if; end do; DE; end proc:

#d(n):=|Der(n)|
d:=proc(n) option remember: if n=1 then 0; else n*d(n-1)+(-1)^n; end if; 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(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)}: end do; S; end proc:

############################################################################################################################################################
#END REFORMATTED C3

#BEGIN REFORMATTED C26
############################################################################################################################################################

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): end do; [seq(coeff(f,x,i),i=1..N)]; end proc:

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): end do; [seq(coeff(f,x,i1),i1=1..N)]; end proc:

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

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

Mul:=proc(pi,sig): local i: [seq(sig[pi[i]],i=1..nops(pi))]; end proc:

#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)}: end do; G2; end proc:

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

############################################################################################################################################################
#END REFORMATTED C26

#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?
 
#Corrected CubeGp()

Cube:=proc(): local G, c, k:
G:=CubeGp():
expand(subs({seq(x[k]=c, k=1..6)},CIP(G, x))/nops(G));
end proc:

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

Necklace:=proc(n): local a, b, G, c, i, k:
a:=[seq(i, i=1..n)]:
b:=[seq(n+1 - i, i=1..n)]:
G:=GenGp({a,b}):
expand(subs({seq(x[k]=c, k=1..n)},CIP(G, x))/nops(G));
end proc: